- Email: [email protected]
Lithospheric necking: a dynamic model for rift morphology
Earth and Planetary Science Letters, 77 (1986) 373-383
373
Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
[61
Lithospheric necking: a dynamic model for rift morphology M.T. Z u b e r * and E.M. P a r m e n t i e r Department of Geological Sciences, Brown Uniuersity, Providence, RI 02912 (U.S.A.) Received August 26, 1985; revised version received January 3, 1986 Rifting is examined in terms of the growth of a necking instability in a lithosphere consisting of a strong plastic or viscous surface layer of uniform strength overlying a weaker viscous substrate in which strength is either uniform or decreases exponentially with depth. As the lithosphere extends, deformation localizes about a small imposed initial perturbation in the strong layer thickness. For a narrow perturbation, the resulting surface topography consists of a central depression and uplifted flanks; the layer thins beneath the central depression. The width of the rift zone is related to the dominant wavelength of the necking instability, which in turn is controlled by the layer thickness and the mechanical properties of the lithosphere. For an initial thickness perturbation with a width less than the dominant wavelength, deformation concentrates into a zone comparable to the dominant wavelength. If the initial perturbation is wider than the dominant wavelength, then the width of the zone of deformation is controlled by the width of the initial perturbation; deformation concentrates in the region of enhanced thinning and develops periodically at the dominant wavelength. A surface layer with limiting plastic (stress exponent n = ~ ) behavior produces a rift-like structure with a width typical of continental rifts for a strong layer thickness consistent with various estimates of the m a x i m u m depth of brittle deformation in the continental lithosphere. The width of the rift is essentially independent of the layer/substrate strength ratio. For a power law viscous surface layer (n = 3), the dominant wavelength varies with layer/substrate strength ratio to the one-third power and is always larger than for a plastic surface layer of the same thickness. The great widths of rift zones on Venus may be explained by unstable extension of a strong viscous surface layer.
1. Introduction
Rift zones are areas of localized lithospheric extension characterized by a central depression, uplifted flanks, and thinning of the underlying crust. These features have been identified on many of the planets and their satellites; on the earth, rifts are found both on the continents and in ocean basins, and represent the initial stage of continental breakup and seafloor spreading. Topographic profiles illustrating the general morphology of rifts on continents and on the surface of Venus are shown in Fig. 1. High heat flow, broad regional uplift, and local magmatism are often associated with rifts, which suggests that thermal as well as mechanical effects play a role in determining their morphology. In this study, we investigate the mechanical aspects of rifting by developing a simple model in which the dynamical consequences of flow in an extending plastic or * Now at: Geodynamics Branch, N A S A / G o d d a r d Space Flight Center. 0012-821X/86/$03.50
© 1986 Elsevier Science Publishers B.V.
viscous lithosphere are evaluated. For simple rheological stratifications we calculate the pattern of near-surface deformation which arises due to horizontal extension and compare the results to major morphological features of rifts, such as those shown in Fig. 1. The width and morphology of rift zones was first explained by Vening Meinesz [1], who suggested that rifts form in an extending elastic-brittle layer which fails by normal faulting. Flexure of the layer occurs in response to motion on a normal fault and a second normal fault forms where the bending stresses in the layer are a maximum, thus defining the width of the rift. Flanking highs form by isostatic upbending of the elastic layer in response to graben subsidence. The elastic properties of the layer thus determine the width and morphology of the rift. Artemjev and Artyushkov [2] qualitatively considered rifting due to necking in a ductile crust which is strong near the surface and weaker at depth. A perturbation in lower crustal thickness localizes during uniform extension and a corresponding stress concentration in the brittle
374 (a) GULF OF SUEZ
2ooo -7
1'
~
,~o ~
RHINE GRABEN
RIO GRANDE RIFT
I 0
2"r~" A~
.,o
I tO0 DISTANCE (km)
286". 37N
'F"
,Oo
I
_....~.
i
"3,
oE.
B~
.o
-/
CD,D~
4
/ /
....... DARK
G~
BRIGHT- LINEAR E D G E OF B R I G H T A R F A
- -
.. o I
I
C'
0
1
,
/
E~ FD,.
Hm,-
I 200
BRIGHT
o~"-
:,
I
I
,':
\
J~ dD,-
ot::~-
:E'
""~"~°"" I
--q:
Oe.-~
d . ' 4,00
(b)
"'1o .,I 0
T ',"°,~ .
. -?"
Fig. 1. (a) Topographic profiles across several continental rift zones. From Buck [15]. (b) Map of radar bright and dark linear features, interpreted as faults, and topographic profiles across rift zone in Beta Regio, Venus. Note the presence of a central trough bounded by uplifted flanks. The arrows represent the approximate locations of major bounding faults shown on map. Vertical lines mark edge of the map. The locations of Rhea and Theia Mons, interpreted as volcanic shields, are shown for reference. From Campbell et al. [251.
375
upper crust results in a narrow zone of normal faulting. Bott [3] developed a model for rifting by crustal stretching in response to external tension on the basis of the elastic model of Vening Meinesz. In this model, the elastic upper crust responds to tension by normal faulting and graben subsidence while the ductile lower crust undergoes thinning due to horizontal flow beneath the subsiding graben. In this study, we examine rifting as the growth of a necking instability. To nucleate a rift, a small thickness perturbation is imposed at the base of a strong layer which overlies a weaker substrate. For a range of rheological parameters we evaluate the conditions for which the initial disturbance will amplify as the lithosphere extends, and we determine the associated pattern of near-surface deformation. As shown later (see Fig. 5), the pattern of deformation consists of a central depression beneath which the strong layer thins due to necking, and flanking uplifts, all of which are characteristic of rift zones. Horizontal extension of the lithosphere, which drives the instability, could occur in response to remotely applied forces or horizontal forces due to regional doming on a scale much larger than the width of the rift zone. The model results are thus applicable to both passive and active rifting (e.g. [4]). The initial thickness perturbation could correspond to a preexisting structural weakness such as a tectonic suture or a thermal anomaly due to an igneous intrusion, diapiric upwelling, or convective transport of heat to the base of the lithosphere. All of these have been suggested in association with terrestrial rifts. A thickness perturbation due to diapirism or convective upwelling may be much wider than the rift, while that due to igneous intrusion or a pre-existing weakness may be narrower than the rift. As discussed later, differences in the width of the perturbation may have implications for the character of rifting. 2. Model formulation 2.1. One-dimensional extension
The simplest formulation for the growth of an instability in an extending layer treats the flow as locally one-dimensional with a uniaxial stress Oxx and strain rate c** [5]. Consider a layer of thick-
ness h ( x ) which has a viscosity much higher than its surroundings. To satisfy the condition of equilibrium in the layer, the horizontal force F where: (1)
F=axx h
must be independent of x. To satisfy the condition of incompressibility: Cxx = - h - I
Oh/Ot
(2)
where the minus sign indicates that the layer thins as it extends. The constitutive relationship between stress and strain rate is:
,xx =AoL
(3)
where n is the stress exponent of the layer and A is a constant. By combining (1) to (3), the rate of change of layer thickness can be written: Oh/Ot = - A F " h 1-"
(4)
For later comparison with our linearized two-dimensional formulation, h can be expressed: h = n ( t ) - 8 ( x , t)
(5)
where H ( t ) is the layer thickness for uniform thinning and 8 is a small (<< h) thickness variation. By substituting (5) into (4), expanding, and integrating with respect to time, the amplitude of the layer thickness variation may be written: ~(x, t) = 6 ( x , 0) e x p [ ( n - 1),xxt ]
(6)
where 8(x, 0) describes the shape of the initial perturbation and Cxx is independent of time. The stress exponent n defines the growth rate of a thickness perturbation. If n > 1 an initial disturbance will grow with time, while if n < 1 a disturbance will decay and the layer will thin uniformly. Equation (6) demonstrates that n > 1 is necessary for extensional instability and that the horizontal layer thickness variation retains the same shape as the perturbation grows in amplitude. 2.2. Linearized two-dimensional extension
In a two-dimensional formulation, the total flow can be expressed as the sum of a mean flow or basic state of uniform horizontal extension and a perturbing flow which arises due to instability. The lithosphere is represented as a strong layer of thickness h overlying a weaker substrate. Two strength stratifications, illustrated in Fig. 2, are considered. For creep deformation, the strain rate
376
d-Model
C-Model
viscous fluid:
-fh
+ R=Ga,d
oij = 21acij - p3ij
(8)
where # is the dynamic viscosity and p is the pressure. Substituting the stresses and strain rates from (7) into (8) and linearizing about the basic state gives:
6,, X = ( 2 ~ t / n ) g x x - ~ ;/h
G = 6xz = 2~gxz
(9)
Fig. 2. Strength stratification of strength jump (J) and continuous strength (C) models. The former is described in terms of the ratio of layer/substrate strength R = #~)/#~2), while the latter is described by the ratio a = ~/h of the viscosity decay depth in the substrate to the strong layer thickness. The layer and substrate are described by uniform densities (O> P2) and power law exponents (nl, n2).
where ~ is the viscosity evaluated at the stress or strain rate of the basic state. For a single Fourier harmonic, perturbing velocities in the vertical and horizontal directions:
(c) is proportional to stress (o) to a power n I = 3. Deformation by distributed faulting can be idealized using a perfectly plastic material with a stress exponent of n I = m. The substrate is assumed to deform by creep with n 2 = 3. In the strength j u m p (J) model, the layer and substrate both have a uniform strength with Oxx~'(1)> Uxx~'(2),where the superscripts 1 and 2 refer to the layer and substrate, respectively. In the continuous strength (C) model, the strength is uniform in the layer, continuous across the layer-substrate interface, and decays exponentially in the substrate with an e-folding depth ~'. Since we are interested primarily in dynamic, stress-supported topography and numerous previous studies have considered the effect of isostatically compensated crustal thickness variations, the layer and substrate are each represented by a uniform density (p) with I)1 = 02. In a layered medium, any small perturbation along an interface will cause deviatoric stresses proportional to the amplitude of the perturbation resulting in the growth of an instability. The total velocities (u, w), stresses (oo), and strain rates ( q j) are written:
where D = d / d z and k ( = 27r/)~) is the wave number, satisfy the incompressibility condition. Within any layer in which the viscosity varies exponentially with depth, the equations of equilibrium are satisfied if:
u=~+fi
w=~+~
% = o-,j + a,j
(7)
Eij = Eij + Eij
where a bar represents the mean flow and a tilde the additional or perturbing flow. For an isotropic
= Wcos kx =-k
(10)
1DWsin k x
D4W + 2~-lD3W+
-2k2~-'(2/,-
[~-2 _ 2 k 2 ( Z / n _ 1)] D 2 W 1 ) D W + k 2 [ k 2 + ~,-2] W = 0 (11)
Solutions of (11) are given in Fletcher and Hallet [6]. The flow in any layer with uniform strength or viscosity can be obtained from (11) by taking the limit as ~"+ ~ . The velocities u and w and the stress components %, and ozz must be continuous across each interface. At the surface, the vertical normal stress must equal the weight of the surface topography and the shear stress must vanish. Expressed in terms of the basic and perturbing stresses, the stress continuity conditions at each interface give:
5x~(x, d i ) = [~x~/-') - ~ i ' ] d 3 , / d x
G(x,
(12)
4) = (p,-, - 0,)ga,
where d i is the depth of the interface. These equations have been linearized for small distortions (3i) of the interface (cf. [7,8]). Vertical velocities arising from these perturbing stresses amplify initial deformation of the interfaces further enhancing the perturbing stresses resulting in instability. At a given time, the shape of the ith inter-
377
face can be represented by the superposition of Fourier harmonics:
6i(x, t ) = ~ A i ( k ,
t) sin
kx
(13)
k
where i = 1, 2 refer to the free surface and the layer-substrate interface, respectively, and 8 and A are the spatial and wave number domain representations of the interface shapes. In terms of the perturbing vertical velocity, the time rate of change of amplitude A is:
A, = W(k,
d~) - ~x~A,
(14)
where A = d A / d t . The system of equations consisting of (11), (14) and that obtained by substituting (13) into (12) has solutions of the form Acc exp[(q - 1)~xt ], where q is the growth rate factor which measures the degree of instability of the disturbance at a particular wave number. The value of unity subtracted from q represents the kinematic interface distortion from the second term on the right-hand side of (14). Complete ~olutions of this system of equations can be written: Al(k, t) = A,1 e x p [ ( q , -
1)exx ]
+A12 exp[(q2 - 1)exx]
(15)
A2(k, t) =A21 exp[(q 1 - 1)exx ] +A22 exp[(q2 - 1)exx] where ql and q2 are eigenvalues, (An, A12 ) and (A21 , A22) are the eigenvectors corresponding to ql and q2, respectively, and exx ( = i x x t) is the total horizontal extension at a time t. The length of the eigenvectors are determined from the prescribed initial conditions. The surface topography is assumed to be initially flat [81(x, 0 ) = 0], and the initial shape of the layer-substrate interface is given by a prescribed function. The Fourier transform of the initial interface shapes is substituted into (15), and the perturbed surface topography is found from the inverse Fourier transform of Al(k, t). Patterns of deformation in the layer and substrate can be determined by calculating the displacements due to the perturbing flow at points on a coordinate grid. The number of growth rate factors at each wave number (in this case two) equals the number of interfaces. In general, only one value of q is
greater than unity and thus contributes to instability. While both values of q are required to completely describe the evolution of interface shapes, solutions of (15) which retain only the positive value of q yield a good approximation to the perturbed interface shape. Note the correspondence of (15) to the one-dimensional growth rate equation (6) with the substitution of q for n. As subsequent results for the two-dimensional problem show, q varies directly with the stress exponent in the layer, but unlike the one-dimensional growth rate, q is a function wave number. The J and C models are each expressed in terms of two dimensionless parameters. In both models, S = ( P a - Po)gh/[a~lx)/2] where P0 is the density above the layer, determines the relative importance of surface topography which stabilizes the deformation and stresses due to surface distortion which help drive the instability. In the J model, R = ~1)/7,(2) the layer/substrate strength ratio, conx/Uxx, trols the relative magnitude of stresses due to distortion of the surface and the layer-substrate interface. In the C model, the layer-substrate interface does not help to drive the instability because the strength is continuous across it. The amount of viscous resistance to deformation in the substrate is determined by the ratio of the viscosity decay depth to the layer thickness a = ~/h. 3. Results
In the limit of large layer/substrate viscosity contrast, the one-dimensional growth rate agrees with the growth rate of the dominant wavelength for the two-dimensional formulation, as noted by Emerman and Turcotte [9]. However, Fig. 3a shows that the growth rate spectrum for one-dimensional extension is independent of wave number so this model fails to predict a dominant wavelength. Thus, to explain the regular development of geologic structures using one-dimensional theory, initial thickness perturbations along layer interfaces must be distributed periodically. In the absence of a process that defines a natural length scale, the occurrence of periodic initial disturbances seems unlikely. Fig. 3 also shows growth rate spectra for twodimensional models with limiting plastic (n 1 = 10 4) and power law viscous (nl = 3) surface layers plotted as a function of dimensionless wave number k' (=2~rh/)Q. The wave number (k~) at
378
a)
b}
r, ~..-3
1.5 L
one-d xO.5
/~=m
I0
a) d-Model ......... ~ .
/
"I | , ~
-- nl= -- nt=3 ......... xlh
ql.O
~
1
0
_
2
b)
-I0 0.5
o
!
d model
\ k'
-2C
r~:~/[
-3o
[ l
C model
Fig. 3. (a) Growth rate spectra as a function of wave number non-dirnensionalized by the strong layer thickness k ' ( = kh). The growth rate for one-dimensional flow in a layer with stress exponent n = 3 is simply equal to n over the entire range of k'. Thus the one-dimensional model fails to predict a dominant wave number. For two-dimensional flow, the solid and dashed lines refer to cases with limiting plastic (n] = 104) and power law viscous (n 1 = 3 ) surface layers, respectively, for the J model. For both cases S = 1.2, R = 50 and n 2 = 3. The peaks in the growth rate spectra define the dominant wave numbers. (b) Two-dimensional growth rate spectra for the C model for plastic and power law viscous surface layers with S = 1.2, a = 5 and n 2 = 3. In this model the growth rate for a viscous surface layer is always less than unity and therefore the layer is stable with respect to necking.
which the growth rate factor is maximized ( q d ) defines the dominant wavelength )kd ( = 2 ' n ' / k d ) . For the J model (Fig. 3a), the magnitude of the growth rate for a lithosphere with a plastic surface layer is approximately two orders of magnitude greater than that for a power law viscous layer with the same value of R. This indicates that a lithosphere in which the near-surface deforms plastically is much more unstable in extension than a lithosphere which deforms viscously throughout. The growth rate spectrum for the C model is shown in Fig. 3b for both plastic and viscous surface layers. For n I = 3 the growth rate is everywhere less than unity, which indicates that the lithosphere is stable with respect to necking. An initial disturbance in strong layer thickness will decay with time, and deformation of the medium will be manifest as uniform thinning. This holds for any reasonable range of lithosphere mechanical properties. Topographic profiles for the two-dimensional growth rate spectra in Fig. 3 and an initial thickness perturbation o f the form ~2(X, O)(X e (~/d)2 where d << h are shown in Fig. 4. The topography in each case consists of a central depression and flanking uplifts which are dynamically supported
I
Fig. 4. (a) Rift morphology for limiting plastic (n 1 = 104, solid line) and power law viscous (nl = 3, dashed line) surface layers for the strength j u m p (J) model. For both cases S = 1.2, R = 50 and n 2 = 3. For comparison, the surface topography in each case is normalized by the depth of the central depression. (b) Rift morphology for the C model with plastic surface layer and with S = 1 . 2 , c t = 5 and n 2 = 3 .
and arise in response to the perturbing flow induced as the layer necks. For the J model, a plastic surface layer (n] = 10 4) results in deformation localized into a region overlying the initial thickness perturbation. With a power law viscous surface layer (hi = 3) deformation is more broadly distributed. For the C model, the relative amount of uplift of the flanks is much greater than in the J model, and small surface depressions occur outside the flanking highs. A variety of initial perturbation shapes including Gaussian, boxcar, and triangular functions have been examined. The topography is independent of the shape of the initial perturbation. The amplitude of the deformation is dependent on the amplitude of the initial disturbance, the amount of horizontal extension exx, and the magnitude of the growth rate factor. In the linearized theory, the amplitude of deformation for a given horizontal extension is exactly proportional to the amplitude of the initial disturbance. Since for actual rift zones the amplitude of this disturbance is unknown, only relative amplitudes, which depend on the growth rate spectrum and horizontal extension, are shown in Fig. 4. Fig. 4a shows that the flanking uplifts are narrower and have a greater amplitude relative to the depth of the central depression for the J model with large n t. For the same initial disturbance amplitude, the absolute magnitude of the topography for the plastic surface layer is greater than that for the viscous layer, reflecting the larger growth rate for the former case as shown in Fig. 3a. As for the topography discussed above, the
379
IIllllllllllll~ z']']l']'~'l I',li
.. ~ l l l l l l l l l l l l
IIIIIIlIIl]l
Illlllllllllllllll~
_/x
"~4.UIIIIIIIIIIIIIIIII
VE-Zx
Fig. 5. Deformation of an extending plastic surface layer with an initial thickness perturbation at its base. The width of the disturbance is determined by the growth rate spectrum in Fig. 3a. The shape of the initial perturbation of width d (<< h) is shown schematically by the dashed line. The central depression and uplifted flanks produced due to the unstable growth of the initial perturbation are characteristic topographic features of rift zones (cf. Fig. 1).
pattern of strong layer deformation shown in Fig. 5 reflects the shape of the growth rate spectrum. While the absolute amplitude is arbitrary, the relative amplitudes are everywhere properly represented. Hence if the topography at the surface or another interface is known, the amplitude of deformation at any point in the layer or substrate is determined. For example, Fig. 5 shows that the amount of upwarping at the base of the strong layer is approximately four times the depth of the central depression. Likewise, the amount of relative uplift of the flanks is about half the depth of the central depression. The width of the rift zone is controlled by the dominant wavelength. For the J model, the dominant wave numbers shown in Fig. 3 correspond to wavelength to layer thickness ratios of Xd/h = 4.0 and 10.6, respectively, for plastic and viscous surface layers. In the J model, the width of the rift from flank-to-flank is almost exactly equal to the dominant wavelength. In the C model for n t 10 4, the dominant wave number corresponds to ~d/h = 3.7, which is very similar to the result in the J
I0°
K~
102
R
i0~
I0
V3
I/6
V9
VI2
I/f5
Va
Fig. 7. (a) Dominant growth rate factor qd vs• R for the J model• (b) qd VS. 1/a for the C model• For r/1 = 3 , q is always less than unity and the lithosphere extends uniformly. S = 0 and n 2 = 3 for both cases.
model with a plastic surface layer. However, the flank-to-flank width of the rift zone for the C model is about twice the dominant wavelength. In this case the width of the central depression is approximately equal to the dominant wavelength. The relationships between rift zone width and the dominant wavelength demonstrated by the above examples hold for the complete range of cases we have examined for the J and C models. Fig. 6 shows the variation of k~ with R for the J model. For power law viscous near-surface behavior k d varies a s R -]/3. This relationship was first determined for flexural buckling of a layered medium [10] and was later shown to be valid for boudinage and folding of a layered viscous medium in the limit of large viscosity contrast [7,8,11]. In the plastic limit k~ is essentially independent of R. For this case the width of the rift will be
=
1.8 '~'- -~v ........
-' -' -" -' ~,6o'8°
15 ~0 ~k' 1.2t q ~ ' X , ~
nl=~o
%[
ioo
- - n,'lO'~ - - n,-3
120
% ,oo
~
0.9[
80
0.6
SO 40
•
~ - --Tkd. . . . 0
IOI
102
R
IOa
Fig. 6. D o m i n a n t wave number k~ as a function of strength contrast R for the J model with S = 0 and n 2 = 3.
2
4
20
. . . . . . . . . .
6
8
I0
S
12
14
Fig. 8. Dominant wave number k~l and dominant growth rate factor qd as a function of S for n] =104 (solid) and n 1 = 3 (dashed) with R = 50 and n 2 = 3.
380 primarily a function of the strong layer thickness. The dominant wavelength decreases with increasing stress exponent in the layer (n 1), which is consistent with results on folding and boudinage of a non-Newtonian layer embedded in a viscous medium [11]. The variation of the growth rate at the dominant wave number with R for the J model and with 1/a for the C model is shown in Fig. 7. Note that qd increases with R and l / a , indicating enhanced instability for large strength contrasts and small e-folding depths. Increased instability is a result of the relative decrease in viscous resistance in the substrate for both cases. In the plastic layer case for the J model the magnitude of qa is generally proportional to R, indicating that the strength contrast at the layer-substrate interface provides a significant contribution to driving the instability. Fig. 8 summarizes the effects of the b u o y a n c y / strength parameter S for the J model and shows that k~ and qa decrease with increasing S for both large and small n~. As would be expected, an increase in density with depth across an interface stabilizes the perturbing flow. The decrease in dominant wave number with increasing S is consistent with results for folding of a layered viscoelastic medium under the influence of gravity [12]. Increasing S results in deepening of the rift and damping of the bounding flanks. 4. Discussion
as a mechanical response to localized extension. In the present model the flanks are solely a consequence of viscous or plastic flow. A thermal anomaly that is broader than the region of crustal thinning may also explain flanking uplifts. Buck [15], Keen [16], and Steckler [17] suggest that thermally-produced uplifted flanks form as a result of horizontal heat transfer due to small-scale convection induced by the rift temperature structure. Steckler [17] interpreted the amount of uplift in the Gulf of Suez to be indicative of lithospheric heating greatly in excess of that expected due to uniform lithosphere extension. Buck [15] showed that the magnitude and distribution of uplift associated with continental rifts may be explained by small-scale convection. The relative contributions of thermal and mechanical mechanisms should be reflected by the timing of uplift on flanks. If uplift is dynamic, then deepening of the central trough and uplift of the flanks should occur simultaneously. If it is a thermal effect, then uplift should lag the formation of the central depression. Flanking uplifts along the Gulf of Suez formed during the main phase of rift development [17], but the relative timing of the formation of the rift valley and the flanks is 'not yet clearly defined. Thermal uplift, unless it is frozen in by a thickening elastic lithosphere, should decay due to cooling after extension ceases. Uplift due to viscous or plastic flow, produced while the extending lithosphere is at yield, must also be frozen-in as stresses fall below the yield stress at the cessation of extension.
4.1. Flanking uplifts 4.2. Continental rifts Uniform stretching of the lithosphere resulting in crustal thinning and subsequent conductive cooling may account for the subsidence of rift basins (e.g. [13]) but does not explain the flanking highs characteristic of rift zones. Flanking uplifts in our necking model form in response to dynamic upwarming produced by the unstable flow. Since there are no density contrasts at depth, the flanks are supported entirely by stresses in the layer. As in the models of Vening Meinesz [1] and Bott [3] which treat the near-surface as an elastic layer, the flanking uplifts regionally compensate the rift depression. Finite element solutions [14] incorporating elastic as well as plastic and viscous behavior have also shown that flanking uplifts can occur
The width of a rift zone formed by necking is a function of the thickness of the strong layer of the lithosphere. The strength stratification of the lithosphere is determined by brittle behavior at shallow depths and creep at greater depths (cf. [18]). The brittle strength increases linearly with depth to a maximum value determined by a creep strength which decreases with depth as temperature increases. Neither the J nor C model is an exact representation of this strength stratification since in both models the linear increase of brittle strength with depth is approximated by a layer of uniform strength. In addition, in the C model the strength falls to zero at depth while in the J model the
381 strength is discontinuous at the base of the strong layer. For the purposes of the present study, the strong layer thickness will be taken to correspond to the depth of the brittle-ductile transition. In a model lithosphere with a plastic surface layer, the layer-substrate interface corresponds to a change in power law exponent, indicating the change in the mode of deformation from faulting to ductile flow. Strength decreases rapidly with depth below the brittle-ductile transition. If estimated from this strength stratification, the strong layer thickness would not be significantly greater than the depth to the brittle-ductile transition. As shown in Fig. 4, a rift formed in a lithosphere with a plastic surface layer has a flank-toflank width = 4h. If the strong upper crustal region of the continental lithosphere is best described by a plastic theology, then a continental rift with a typically observed width in the range 35-60 km [19] requires a strong layer thickness of 9-15 km. An extending continental lithosphere in which a quartz rheology approximates flow in the crust undergoes the transition from brittle to ductile behavior at a depth of about 15-20 km depending on the geothermal gradient and strain rate (cf. [18]). This is comparable to the depth to which earthquakes are observed in the continental lithosphere and to which faulting penetrates in continental rifts as shown in seismic reflection profiles. For a strong viscous surface layer with nl = 3, the dominant wavelength, and therefore the rift zone width, depends on the ratio of layer to substrate strength as shown in Fig. 6. For a rift zone 60 km wide, typical of the well-studied East African Rift or Rhinegraben, a strong layer thickness of 15 km requires an R only slightly greater than unity. The width of the rift is controlled by the dominant wavelength as long as the initial perturbation width d < ?~d. However, if d > ?~a the width of the rift is governed by the width of the initial disturbance. If an initial perturbation is much wider than the dominant wavelength, corresponding perhaps to broad-scale lithospheric thinning, extension localizes in the thinned region and deformation within that region develops periodically at the dominant wavelength [6,20]. The vast regional extent of the Basin and Range Province ( > 103 km) may reflect an initial strong layer thickness per-
turbation much wider than the spacing of individual basins and ranges ( -~ 30 km). The depth distribution of seismicity [21], high heat flow [22], and low Pn velocities [23] in this area suggest that an upper mantle thermal anomaly of large horizontal extent could be the cause of the initial perturbation. 4.3. Rifting on Venus
A number of features on Venus with the morphology of rift zones have been recognized from Pioneer Venus radar altimetry [24]. High-resolution earth-based radar images have recently been obtained for a major rift in Beta Regio [25]. As shown in the across-strike profiles in Fig. lb, the topography consists of a central depression and bounding highs which may represent uplifted flanks. The width of the rift varies along strike but averages about 150 km from flank to flank. Two volcanic shields occur along the rift, Rhea Mons in the north and Theia Mons in the south. While volcanic construction could contribute to the flanking highs, they are best developed where the central depression is deepest and in areas not associated with the volcanoes (profiles G and H in Fig. lb). The high surface temperature of Venus, approximately 700 K, suggests that ductile deformation will occur at shallower depths than on the earth. On the basis of estimates of rock strength, surface temperature, and the spacing of features of presumed tectonic origin in the banded terrain, Solomon and Head [26] suggest that the elastic lithosphere thickness on Venus is in the range 1-10 km, and therefore too thin to explain the width of the Beta Regio rift by elastic flexure, which would require an elastic-brittle layer thickness in excess of 60 km. Schaber [24] estimated a brittle layer thickness in the range 47-69 km on the basis of the widths of other presumed rift valleys on Venus. Radar bright lineaments in Beta Regio, interpreted as faults [25], suggest the presence of a brittle surface layer. On the basis of Fig. 4, the width of the Beta Regio rift would require a plastic layer thickness of about 40 km. If the near-surface of Venus is weaker than the earth, S will be greater than that assumed in Fig. 4. The limit of large S corresponds to a rift zone width = 6h and a plastic layer thickness of 25 km. This
382
is greater than the depth to the brittle-ductile transition on Venus estimated by Solomon and Head [26] and would require a geothermal gradient less than that which they assume (20 K k m - l ) . The 10-20 km spacing of bright lineaments where they are most well-developed is comparable to the width of radar bright bands in the banded terrain. If the radar bright features are due to faulting with a spacing comparable to the brittle-elastic layer thickness, then this thickness in Beta Regio should be similar to that in areas of banded terrain and too thin to explain the rift zone width. If the plastic region of near-surface faulting is thin, then lithospheric necking may be controlled by a thicker viscous layer deforming by creep. In this case the width of the rift zone is determined by the strong layer thickness and the layer/substrate strength contrast. On the basis of Fig. 6, 10 and 20 km thick strong layers require values of R = 100 and 10, respectively, to explain the width of the Beta Regio Rift. 5. Summary Unstable extension of a lithosphere consisting of a strong surface layer overlying a weaker viscous substrate results in a pattern of deformation that is consistent with the major morphological characteristics of rift zones. The surface topography consists of a central depression, beneath which the layer thins by necking, and flanking uplifts. The rift is nucleated by introducing a small amplitude thickness perturbation at the base of the layer. For an initial perturbation narrower than the dominant wavelength, deformation concentrates in a zone of width comparable to the dominant wavelength. For an initial thickness perturbation wider than the dominant wavelength, deformation develops periodically at the dominant wavelength in the region above the perturbation. The dominant wavelength is controlled by the layer thickness and by the growth rate spectrum of extensional instability, which is a function of the layer/substrate strength contrast, the density stratification, and the stress exponents describing flow in the layer and substrate. Extension of a strong surface layer that deforms by distributed faulting, idealized by plastic behavior (n 1 = oo), results in a rift zone width approximately four times the plastic layer thickness. If the base of the
strong layer corresponds to the maximum depth of brittle deformation, determined either by experimental flow laws or the observed depth of earthquakes, then the width predicted by this model is generally consistent with that of a typical continental rift zone. Dominant wavelength and therefore rift zone width for a plastic surface layer is essentially independent of the layer/substrate strength contrast. For a strong viscous surface layer (n~ = 3), dominant wavelength varies as the strength ratio to the one-third power. Extension of a viscous surface layer always produces a rift zone that is wider than that for a plastic surface layer. The great width of rift zones on Venus, relative to the spacing of other features of presumed tectonic origin, may be explained by ductile necking of a strong viscous layer. Acknowledgements This research was supported by NASA grant NSG-7605. We thank Ray Fletcher for helpful discussions and Philip England for a constructive review of the manuscript. References 1 F.A. Vening Meinesz, Les "graben" africains r6sultat de compression ou de tension dans la croBte terrestre?, K. Belg. Kol. Inst. Bull. 21, 539-552, 1950. 2 M.E. Artemjev and E.V. Artyushkov, Structure and isostasy of the Baikal Rift and the mechanism of rifting, J. Geophys. Res. 76, 1971. 3 M.H.P. Bott, Formation of sedimentary basins of graben type by extension of continental crust, Tectonophysics 36, 77-86, 1976. 4 A.H.C. SengSr and K. Burke, Relative timing of rifting and volcanism on earth and its tectonic implications, Geophys. Res. Lett. 5, 419-421, 1978. 5 E.W. Hart, Theory of the tensile test, Act. Metall. 15, 351-355, 1967. 6 R.C. Fletcher and B. Hallet, Unstable extension of the lithosphere: a mechanical model for Basin and Range structure, J. Geophys. Res. 88, 7457-7466, 1983. 7 R.C. Fletcher, Folding of a single viscous layer: exact infinitesimal amplitude solution, Tectonophysics 39, 593-606, 1977. 8 R.B. Smith, Unified theory of the onset of folding, boudinage, and mullion structure, Geol. Soc. Am. Bull. 86, 1601-1609, 1975. 9 S.H. Emerman and D.L. Turcotte, A back-of-the-envelope approach to boudinage mechanics, Tectonophysics 110, 333-338, 1984. 10 M.A. Biot, Folding instability of a layered viscoelastic
383
11
12
13 14
15
16
17
18
19
medium under compression, Proc. R. Soc. London, Ser. A 242, 444-454, 1957. R.B. Smith, Formation of folds, boudinage, and mullions in non-Newtonian materials, Geol. Soc. Am. Bull. 88, 312-320, 1977. M.A. Biot, The influence of gravity on the folding of a layered viscoelastic medium under compression, J. Franklin Inst. 267, 211-228, 1959. D. McKenzie, Some remarks on the development of sedimentary basins, Earth Planet. Sci. Lett. 40, 25-32, 1978. H.D. Lynch, Numerical Models of the Formation of Continental Rifts by Processes of Lithospheric Necking, pp. 164-216, Ph.D. Thesis, New Mexico State University, Las Cruces, N.M., 1983. W.R. Buck, Small-scale convection induced by passive rifting: the cause for uplift of rift shoulders, submitted to Earth Planet. Sci. Lett., 1985. C.E. Keen, The dynamics of rifting: Deformation of the lithosphere by active and passive driving mechanisms, Geophys. J. R. Astron. Soc. 80, 95-120, 1985. M.S. Steckler, Uplift and extension at the Gulf of Suez: Indications of induced mantle convection, Nature 317, 135-139, 1985. W.F. Brace and D.L. Kohlstedt, Limits on lithospheric stress imposed from laboratory experiments, J. Geophys. Res. 85, 6248-6252, 1980. I.B. Ramberg and P. Morgan, Physical characteristics and
20
21
22
23
24
25
26
evolutionary trends of continental rifts, Proc. 27th Int. Geol. Congr. 7, 165-216, 1984. M.T. Zuber, E.M. Parmentier and R.C. Fletcher, Extension of continental lithosphere: A model for two scales of Basin and Range deformation, J. Geophys. Res., in press, 1985. G.P. Eaton, The Basin and Range Province: origin and tectonic significance, Annu. Rev. Earth Planet. Sci. 10, 409-440, 1982. A.H. Lachenbruch and J.H. Sass, Models of an extending lithosphere and heat flow in the Basin and Range Province, in: R.B. Smith and G.P. Eaton, eds., Geol. Soc. Am. Mem. 152, 209-250, 1978. R.B. Smith, Seismicity, crustal structure, and intraplate tectonics of the interior of the western cordillera, in: R.B. Smith and G.P. Eaton, eds., Geol. Soc. Am. Mem. 152, 111-144, 1978. G.G. Schaber, Venus: limited extension and volcanism along zones of lithospheric weakness, Geophys. Res. Lett. 9, 499-502, 1982. D.B. Campbell, J.W. Head, J.K. Harmon and A.A. Hine, Venus: volcanism and rift formation in Beta Regio, Science 226, 167-170, 1984. S.C. Solomon and J.W. Head, Venus banded terrain: tectonic models for band formation and their relationship to lithospheric thermal structure, J. Geophys. Res. 89, 6885-6897, 1984.
373
Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
[61
Lithospheric necking: a dynamic model for rift morphology M.T. Z u b e r * and E.M. P a r m e n t i e r Department of Geological Sciences, Brown Uniuersity, Providence, RI 02912 (U.S.A.) Received August 26, 1985; revised version received January 3, 1986 Rifting is examined in terms of the growth of a necking instability in a lithosphere consisting of a strong plastic or viscous surface layer of uniform strength overlying a weaker viscous substrate in which strength is either uniform or decreases exponentially with depth. As the lithosphere extends, deformation localizes about a small imposed initial perturbation in the strong layer thickness. For a narrow perturbation, the resulting surface topography consists of a central depression and uplifted flanks; the layer thins beneath the central depression. The width of the rift zone is related to the dominant wavelength of the necking instability, which in turn is controlled by the layer thickness and the mechanical properties of the lithosphere. For an initial thickness perturbation with a width less than the dominant wavelength, deformation concentrates into a zone comparable to the dominant wavelength. If the initial perturbation is wider than the dominant wavelength, then the width of the zone of deformation is controlled by the width of the initial perturbation; deformation concentrates in the region of enhanced thinning and develops periodically at the dominant wavelength. A surface layer with limiting plastic (stress exponent n = ~ ) behavior produces a rift-like structure with a width typical of continental rifts for a strong layer thickness consistent with various estimates of the m a x i m u m depth of brittle deformation in the continental lithosphere. The width of the rift is essentially independent of the layer/substrate strength ratio. For a power law viscous surface layer (n = 3), the dominant wavelength varies with layer/substrate strength ratio to the one-third power and is always larger than for a plastic surface layer of the same thickness. The great widths of rift zones on Venus may be explained by unstable extension of a strong viscous surface layer.
1. Introduction
Rift zones are areas of localized lithospheric extension characterized by a central depression, uplifted flanks, and thinning of the underlying crust. These features have been identified on many of the planets and their satellites; on the earth, rifts are found both on the continents and in ocean basins, and represent the initial stage of continental breakup and seafloor spreading. Topographic profiles illustrating the general morphology of rifts on continents and on the surface of Venus are shown in Fig. 1. High heat flow, broad regional uplift, and local magmatism are often associated with rifts, which suggests that thermal as well as mechanical effects play a role in determining their morphology. In this study, we investigate the mechanical aspects of rifting by developing a simple model in which the dynamical consequences of flow in an extending plastic or * Now at: Geodynamics Branch, N A S A / G o d d a r d Space Flight Center. 0012-821X/86/$03.50
© 1986 Elsevier Science Publishers B.V.
viscous lithosphere are evaluated. For simple rheological stratifications we calculate the pattern of near-surface deformation which arises due to horizontal extension and compare the results to major morphological features of rifts, such as those shown in Fig. 1. The width and morphology of rift zones was first explained by Vening Meinesz [1], who suggested that rifts form in an extending elastic-brittle layer which fails by normal faulting. Flexure of the layer occurs in response to motion on a normal fault and a second normal fault forms where the bending stresses in the layer are a maximum, thus defining the width of the rift. Flanking highs form by isostatic upbending of the elastic layer in response to graben subsidence. The elastic properties of the layer thus determine the width and morphology of the rift. Artemjev and Artyushkov [2] qualitatively considered rifting due to necking in a ductile crust which is strong near the surface and weaker at depth. A perturbation in lower crustal thickness localizes during uniform extension and a corresponding stress concentration in the brittle
374 (a) GULF OF SUEZ
2ooo -7
1'
~
,~o ~
RHINE GRABEN
RIO GRANDE RIFT
I 0
2"r~" A~
.,o
I tO0 DISTANCE (km)
286". 37N
'F"
,Oo
I
_....~.
i
"3,
oE.
B~
.o
-/
CD,D~
4
/ /
....... DARK
G~
BRIGHT- LINEAR E D G E OF B R I G H T A R F A
- -
.. o I
I
C'
0
1
,
/
E~ FD,.
Hm,-
I 200
BRIGHT
o~"-
:,
I
I
,':
\
J~ dD,-
ot::~-
:E'
""~"~°"" I
--q:
Oe.-~
d . ' 4,00
(b)
"'1o .,I 0
T ',"°,~ .
. -?"
Fig. 1. (a) Topographic profiles across several continental rift zones. From Buck [15]. (b) Map of radar bright and dark linear features, interpreted as faults, and topographic profiles across rift zone in Beta Regio, Venus. Note the presence of a central trough bounded by uplifted flanks. The arrows represent the approximate locations of major bounding faults shown on map. Vertical lines mark edge of the map. The locations of Rhea and Theia Mons, interpreted as volcanic shields, are shown for reference. From Campbell et al. [251.
375
upper crust results in a narrow zone of normal faulting. Bott [3] developed a model for rifting by crustal stretching in response to external tension on the basis of the elastic model of Vening Meinesz. In this model, the elastic upper crust responds to tension by normal faulting and graben subsidence while the ductile lower crust undergoes thinning due to horizontal flow beneath the subsiding graben. In this study, we examine rifting as the growth of a necking instability. To nucleate a rift, a small thickness perturbation is imposed at the base of a strong layer which overlies a weaker substrate. For a range of rheological parameters we evaluate the conditions for which the initial disturbance will amplify as the lithosphere extends, and we determine the associated pattern of near-surface deformation. As shown later (see Fig. 5), the pattern of deformation consists of a central depression beneath which the strong layer thins due to necking, and flanking uplifts, all of which are characteristic of rift zones. Horizontal extension of the lithosphere, which drives the instability, could occur in response to remotely applied forces or horizontal forces due to regional doming on a scale much larger than the width of the rift zone. The model results are thus applicable to both passive and active rifting (e.g. [4]). The initial thickness perturbation could correspond to a preexisting structural weakness such as a tectonic suture or a thermal anomaly due to an igneous intrusion, diapiric upwelling, or convective transport of heat to the base of the lithosphere. All of these have been suggested in association with terrestrial rifts. A thickness perturbation due to diapirism or convective upwelling may be much wider than the rift, while that due to igneous intrusion or a pre-existing weakness may be narrower than the rift. As discussed later, differences in the width of the perturbation may have implications for the character of rifting. 2. Model formulation 2.1. One-dimensional extension
The simplest formulation for the growth of an instability in an extending layer treats the flow as locally one-dimensional with a uniaxial stress Oxx and strain rate c** [5]. Consider a layer of thick-
ness h ( x ) which has a viscosity much higher than its surroundings. To satisfy the condition of equilibrium in the layer, the horizontal force F where: (1)
F=axx h
must be independent of x. To satisfy the condition of incompressibility: Cxx = - h - I
Oh/Ot
(2)
where the minus sign indicates that the layer thins as it extends. The constitutive relationship between stress and strain rate is:
,xx =AoL
(3)
where n is the stress exponent of the layer and A is a constant. By combining (1) to (3), the rate of change of layer thickness can be written: Oh/Ot = - A F " h 1-"
(4)
For later comparison with our linearized two-dimensional formulation, h can be expressed: h = n ( t ) - 8 ( x , t)
(5)
where H ( t ) is the layer thickness for uniform thinning and 8 is a small (<< h) thickness variation. By substituting (5) into (4), expanding, and integrating with respect to time, the amplitude of the layer thickness variation may be written: ~(x, t) = 6 ( x , 0) e x p [ ( n - 1),xxt ]
(6)
where 8(x, 0) describes the shape of the initial perturbation and Cxx is independent of time. The stress exponent n defines the growth rate of a thickness perturbation. If n > 1 an initial disturbance will grow with time, while if n < 1 a disturbance will decay and the layer will thin uniformly. Equation (6) demonstrates that n > 1 is necessary for extensional instability and that the horizontal layer thickness variation retains the same shape as the perturbation grows in amplitude. 2.2. Linearized two-dimensional extension
In a two-dimensional formulation, the total flow can be expressed as the sum of a mean flow or basic state of uniform horizontal extension and a perturbing flow which arises due to instability. The lithosphere is represented as a strong layer of thickness h overlying a weaker substrate. Two strength stratifications, illustrated in Fig. 2, are considered. For creep deformation, the strain rate
376
d-Model
C-Model
viscous fluid:
-fh
+ R=Ga,d
oij = 21acij - p3ij
(8)
where # is the dynamic viscosity and p is the pressure. Substituting the stresses and strain rates from (7) into (8) and linearizing about the basic state gives:
6,, X = ( 2 ~ t / n ) g x x - ~ ;/h
G = 6xz = 2~gxz
(9)
Fig. 2. Strength stratification of strength jump (J) and continuous strength (C) models. The former is described in terms of the ratio of layer/substrate strength R = #~)/#~2), while the latter is described by the ratio a = ~/h of the viscosity decay depth in the substrate to the strong layer thickness. The layer and substrate are described by uniform densities (O> P2) and power law exponents (nl, n2).
where ~ is the viscosity evaluated at the stress or strain rate of the basic state. For a single Fourier harmonic, perturbing velocities in the vertical and horizontal directions:
(c) is proportional to stress (o) to a power n I = 3. Deformation by distributed faulting can be idealized using a perfectly plastic material with a stress exponent of n I = m. The substrate is assumed to deform by creep with n 2 = 3. In the strength j u m p (J) model, the layer and substrate both have a uniform strength with Oxx~'(1)> Uxx~'(2),where the superscripts 1 and 2 refer to the layer and substrate, respectively. In the continuous strength (C) model, the strength is uniform in the layer, continuous across the layer-substrate interface, and decays exponentially in the substrate with an e-folding depth ~'. Since we are interested primarily in dynamic, stress-supported topography and numerous previous studies have considered the effect of isostatically compensated crustal thickness variations, the layer and substrate are each represented by a uniform density (p) with I)1 = 02. In a layered medium, any small perturbation along an interface will cause deviatoric stresses proportional to the amplitude of the perturbation resulting in the growth of an instability. The total velocities (u, w), stresses (oo), and strain rates ( q j) are written:
where D = d / d z and k ( = 27r/)~) is the wave number, satisfy the incompressibility condition. Within any layer in which the viscosity varies exponentially with depth, the equations of equilibrium are satisfied if:
u=~+fi
w=~+~
% = o-,j + a,j
(7)
Eij = Eij + Eij
where a bar represents the mean flow and a tilde the additional or perturbing flow. For an isotropic
= Wcos kx =-k
(10)
1DWsin k x
D4W + 2~-lD3W+
-2k2~-'(2/,-
[~-2 _ 2 k 2 ( Z / n _ 1)] D 2 W 1 ) D W + k 2 [ k 2 + ~,-2] W = 0 (11)
Solutions of (11) are given in Fletcher and Hallet [6]. The flow in any layer with uniform strength or viscosity can be obtained from (11) by taking the limit as ~"+ ~ . The velocities u and w and the stress components %, and ozz must be continuous across each interface. At the surface, the vertical normal stress must equal the weight of the surface topography and the shear stress must vanish. Expressed in terms of the basic and perturbing stresses, the stress continuity conditions at each interface give:
5x~(x, d i ) = [~x~/-') - ~ i ' ] d 3 , / d x
G(x,
(12)
4) = (p,-, - 0,)ga,
where d i is the depth of the interface. These equations have been linearized for small distortions (3i) of the interface (cf. [7,8]). Vertical velocities arising from these perturbing stresses amplify initial deformation of the interfaces further enhancing the perturbing stresses resulting in instability. At a given time, the shape of the ith inter-
377
face can be represented by the superposition of Fourier harmonics:
6i(x, t ) = ~ A i ( k ,
t) sin
kx
(13)
k
where i = 1, 2 refer to the free surface and the layer-substrate interface, respectively, and 8 and A are the spatial and wave number domain representations of the interface shapes. In terms of the perturbing vertical velocity, the time rate of change of amplitude A is:
A, = W(k,
d~) - ~x~A,
(14)
where A = d A / d t . The system of equations consisting of (11), (14) and that obtained by substituting (13) into (12) has solutions of the form Acc exp[(q - 1)~xt ], where q is the growth rate factor which measures the degree of instability of the disturbance at a particular wave number. The value of unity subtracted from q represents the kinematic interface distortion from the second term on the right-hand side of (14). Complete ~olutions of this system of equations can be written: Al(k, t) = A,1 e x p [ ( q , -
1)exx ]
+A12 exp[(q2 - 1)exx]
(15)
A2(k, t) =A21 exp[(q 1 - 1)exx ] +A22 exp[(q2 - 1)exx] where ql and q2 are eigenvalues, (An, A12 ) and (A21 , A22) are the eigenvectors corresponding to ql and q2, respectively, and exx ( = i x x t) is the total horizontal extension at a time t. The length of the eigenvectors are determined from the prescribed initial conditions. The surface topography is assumed to be initially flat [81(x, 0 ) = 0], and the initial shape of the layer-substrate interface is given by a prescribed function. The Fourier transform of the initial interface shapes is substituted into (15), and the perturbed surface topography is found from the inverse Fourier transform of Al(k, t). Patterns of deformation in the layer and substrate can be determined by calculating the displacements due to the perturbing flow at points on a coordinate grid. The number of growth rate factors at each wave number (in this case two) equals the number of interfaces. In general, only one value of q is
greater than unity and thus contributes to instability. While both values of q are required to completely describe the evolution of interface shapes, solutions of (15) which retain only the positive value of q yield a good approximation to the perturbed interface shape. Note the correspondence of (15) to the one-dimensional growth rate equation (6) with the substitution of q for n. As subsequent results for the two-dimensional problem show, q varies directly with the stress exponent in the layer, but unlike the one-dimensional growth rate, q is a function wave number. The J and C models are each expressed in terms of two dimensionless parameters. In both models, S = ( P a - Po)gh/[a~lx)/2] where P0 is the density above the layer, determines the relative importance of surface topography which stabilizes the deformation and stresses due to surface distortion which help drive the instability. In the J model, R = ~1)/7,(2) the layer/substrate strength ratio, conx/Uxx, trols the relative magnitude of stresses due to distortion of the surface and the layer-substrate interface. In the C model, the layer-substrate interface does not help to drive the instability because the strength is continuous across it. The amount of viscous resistance to deformation in the substrate is determined by the ratio of the viscosity decay depth to the layer thickness a = ~/h. 3. Results
In the limit of large layer/substrate viscosity contrast, the one-dimensional growth rate agrees with the growth rate of the dominant wavelength for the two-dimensional formulation, as noted by Emerman and Turcotte [9]. However, Fig. 3a shows that the growth rate spectrum for one-dimensional extension is independent of wave number so this model fails to predict a dominant wavelength. Thus, to explain the regular development of geologic structures using one-dimensional theory, initial thickness perturbations along layer interfaces must be distributed periodically. In the absence of a process that defines a natural length scale, the occurrence of periodic initial disturbances seems unlikely. Fig. 3 also shows growth rate spectra for twodimensional models with limiting plastic (n 1 = 10 4) and power law viscous (nl = 3) surface layers plotted as a function of dimensionless wave number k' (=2~rh/)Q. The wave number (k~) at
378
a)
b}
r, ~..-3
1.5 L
one-d xO.5
/~=m
I0
a) d-Model ......... ~ .
/
"I | , ~
-- nl= -- nt=3 ......... xlh
ql.O
~
1
0
_
2
b)
-I0 0.5
o
!
d model
\ k'
-2C
r~:~/[
-3o
[ l
C model
Fig. 3. (a) Growth rate spectra as a function of wave number non-dirnensionalized by the strong layer thickness k ' ( = kh). The growth rate for one-dimensional flow in a layer with stress exponent n = 3 is simply equal to n over the entire range of k'. Thus the one-dimensional model fails to predict a dominant wave number. For two-dimensional flow, the solid and dashed lines refer to cases with limiting plastic (n] = 104) and power law viscous (n 1 = 3 ) surface layers, respectively, for the J model. For both cases S = 1.2, R = 50 and n 2 = 3. The peaks in the growth rate spectra define the dominant wave numbers. (b) Two-dimensional growth rate spectra for the C model for plastic and power law viscous surface layers with S = 1.2, a = 5 and n 2 = 3. In this model the growth rate for a viscous surface layer is always less than unity and therefore the layer is stable with respect to necking.
which the growth rate factor is maximized ( q d ) defines the dominant wavelength )kd ( = 2 ' n ' / k d ) . For the J model (Fig. 3a), the magnitude of the growth rate for a lithosphere with a plastic surface layer is approximately two orders of magnitude greater than that for a power law viscous layer with the same value of R. This indicates that a lithosphere in which the near-surface deforms plastically is much more unstable in extension than a lithosphere which deforms viscously throughout. The growth rate spectrum for the C model is shown in Fig. 3b for both plastic and viscous surface layers. For n I = 3 the growth rate is everywhere less than unity, which indicates that the lithosphere is stable with respect to necking. An initial disturbance in strong layer thickness will decay with time, and deformation of the medium will be manifest as uniform thinning. This holds for any reasonable range of lithosphere mechanical properties. Topographic profiles for the two-dimensional growth rate spectra in Fig. 3 and an initial thickness perturbation o f the form ~2(X, O)(X e (~/d)2 where d << h are shown in Fig. 4. The topography in each case consists of a central depression and flanking uplifts which are dynamically supported
I
Fig. 4. (a) Rift morphology for limiting plastic (n 1 = 104, solid line) and power law viscous (nl = 3, dashed line) surface layers for the strength j u m p (J) model. For both cases S = 1.2, R = 50 and n 2 = 3. For comparison, the surface topography in each case is normalized by the depth of the central depression. (b) Rift morphology for the C model with plastic surface layer and with S = 1 . 2 , c t = 5 and n 2 = 3 .
and arise in response to the perturbing flow induced as the layer necks. For the J model, a plastic surface layer (n] = 10 4) results in deformation localized into a region overlying the initial thickness perturbation. With a power law viscous surface layer (hi = 3) deformation is more broadly distributed. For the C model, the relative amount of uplift of the flanks is much greater than in the J model, and small surface depressions occur outside the flanking highs. A variety of initial perturbation shapes including Gaussian, boxcar, and triangular functions have been examined. The topography is independent of the shape of the initial perturbation. The amplitude of the deformation is dependent on the amplitude of the initial disturbance, the amount of horizontal extension exx, and the magnitude of the growth rate factor. In the linearized theory, the amplitude of deformation for a given horizontal extension is exactly proportional to the amplitude of the initial disturbance. Since for actual rift zones the amplitude of this disturbance is unknown, only relative amplitudes, which depend on the growth rate spectrum and horizontal extension, are shown in Fig. 4. Fig. 4a shows that the flanking uplifts are narrower and have a greater amplitude relative to the depth of the central depression for the J model with large n t. For the same initial disturbance amplitude, the absolute magnitude of the topography for the plastic surface layer is greater than that for the viscous layer, reflecting the larger growth rate for the former case as shown in Fig. 3a. As for the topography discussed above, the
379
IIllllllllllll~ z']']l']'~'l I',li
.. ~ l l l l l l l l l l l l
IIIIIIlIIl]l
Illlllllllllllllll~
_/x
"~4.UIIIIIIIIIIIIIIIII
VE-Zx
Fig. 5. Deformation of an extending plastic surface layer with an initial thickness perturbation at its base. The width of the disturbance is determined by the growth rate spectrum in Fig. 3a. The shape of the initial perturbation of width d (<< h) is shown schematically by the dashed line. The central depression and uplifted flanks produced due to the unstable growth of the initial perturbation are characteristic topographic features of rift zones (cf. Fig. 1).
pattern of strong layer deformation shown in Fig. 5 reflects the shape of the growth rate spectrum. While the absolute amplitude is arbitrary, the relative amplitudes are everywhere properly represented. Hence if the topography at the surface or another interface is known, the amplitude of deformation at any point in the layer or substrate is determined. For example, Fig. 5 shows that the amount of upwarping at the base of the strong layer is approximately four times the depth of the central depression. Likewise, the amount of relative uplift of the flanks is about half the depth of the central depression. The width of the rift zone is controlled by the dominant wavelength. For the J model, the dominant wave numbers shown in Fig. 3 correspond to wavelength to layer thickness ratios of Xd/h = 4.0 and 10.6, respectively, for plastic and viscous surface layers. In the J model, the width of the rift from flank-to-flank is almost exactly equal to the dominant wavelength. In the C model for n t 10 4, the dominant wave number corresponds to ~d/h = 3.7, which is very similar to the result in the J
I0°
K~
102
R
i0~
I0
V3
I/6
V9
VI2
I/f5
Va
Fig. 7. (a) Dominant growth rate factor qd vs• R for the J model• (b) qd VS. 1/a for the C model• For r/1 = 3 , q is always less than unity and the lithosphere extends uniformly. S = 0 and n 2 = 3 for both cases.
model with a plastic surface layer. However, the flank-to-flank width of the rift zone for the C model is about twice the dominant wavelength. In this case the width of the central depression is approximately equal to the dominant wavelength. The relationships between rift zone width and the dominant wavelength demonstrated by the above examples hold for the complete range of cases we have examined for the J and C models. Fig. 6 shows the variation of k~ with R for the J model. For power law viscous near-surface behavior k d varies a s R -]/3. This relationship was first determined for flexural buckling of a layered medium [10] and was later shown to be valid for boudinage and folding of a layered viscous medium in the limit of large viscosity contrast [7,8,11]. In the plastic limit k~ is essentially independent of R. For this case the width of the rift will be
=
1.8 '~'- -~v ........
-' -' -" -' ~,6o'8°
15 ~0 ~k' 1.2t q ~ ' X , ~
nl=~o
%[
ioo
- - n,'lO'~ - - n,-3
120
% ,oo
~
0.9[
80
0.6
SO 40
•
~ - --Tkd. . . . 0
IOI
102
R
IOa
Fig. 6. D o m i n a n t wave number k~ as a function of strength contrast R for the J model with S = 0 and n 2 = 3.
2
4
20
. . . . . . . . . .
6
8
I0
S
12
14
Fig. 8. Dominant wave number k~l and dominant growth rate factor qd as a function of S for n] =104 (solid) and n 1 = 3 (dashed) with R = 50 and n 2 = 3.
380 primarily a function of the strong layer thickness. The dominant wavelength decreases with increasing stress exponent in the layer (n 1), which is consistent with results on folding and boudinage of a non-Newtonian layer embedded in a viscous medium [11]. The variation of the growth rate at the dominant wave number with R for the J model and with 1/a for the C model is shown in Fig. 7. Note that qd increases with R and l / a , indicating enhanced instability for large strength contrasts and small e-folding depths. Increased instability is a result of the relative decrease in viscous resistance in the substrate for both cases. In the plastic layer case for the J model the magnitude of qa is generally proportional to R, indicating that the strength contrast at the layer-substrate interface provides a significant contribution to driving the instability. Fig. 8 summarizes the effects of the b u o y a n c y / strength parameter S for the J model and shows that k~ and qa decrease with increasing S for both large and small n~. As would be expected, an increase in density with depth across an interface stabilizes the perturbing flow. The decrease in dominant wave number with increasing S is consistent with results for folding of a layered viscoelastic medium under the influence of gravity [12]. Increasing S results in deepening of the rift and damping of the bounding flanks. 4. Discussion
as a mechanical response to localized extension. In the present model the flanks are solely a consequence of viscous or plastic flow. A thermal anomaly that is broader than the region of crustal thinning may also explain flanking uplifts. Buck [15], Keen [16], and Steckler [17] suggest that thermally-produced uplifted flanks form as a result of horizontal heat transfer due to small-scale convection induced by the rift temperature structure. Steckler [17] interpreted the amount of uplift in the Gulf of Suez to be indicative of lithospheric heating greatly in excess of that expected due to uniform lithosphere extension. Buck [15] showed that the magnitude and distribution of uplift associated with continental rifts may be explained by small-scale convection. The relative contributions of thermal and mechanical mechanisms should be reflected by the timing of uplift on flanks. If uplift is dynamic, then deepening of the central trough and uplift of the flanks should occur simultaneously. If it is a thermal effect, then uplift should lag the formation of the central depression. Flanking uplifts along the Gulf of Suez formed during the main phase of rift development [17], but the relative timing of the formation of the rift valley and the flanks is 'not yet clearly defined. Thermal uplift, unless it is frozen in by a thickening elastic lithosphere, should decay due to cooling after extension ceases. Uplift due to viscous or plastic flow, produced while the extending lithosphere is at yield, must also be frozen-in as stresses fall below the yield stress at the cessation of extension.
4.1. Flanking uplifts 4.2. Continental rifts Uniform stretching of the lithosphere resulting in crustal thinning and subsequent conductive cooling may account for the subsidence of rift basins (e.g. [13]) but does not explain the flanking highs characteristic of rift zones. Flanking uplifts in our necking model form in response to dynamic upwarming produced by the unstable flow. Since there are no density contrasts at depth, the flanks are supported entirely by stresses in the layer. As in the models of Vening Meinesz [1] and Bott [3] which treat the near-surface as an elastic layer, the flanking uplifts regionally compensate the rift depression. Finite element solutions [14] incorporating elastic as well as plastic and viscous behavior have also shown that flanking uplifts can occur
The width of a rift zone formed by necking is a function of the thickness of the strong layer of the lithosphere. The strength stratification of the lithosphere is determined by brittle behavior at shallow depths and creep at greater depths (cf. [18]). The brittle strength increases linearly with depth to a maximum value determined by a creep strength which decreases with depth as temperature increases. Neither the J nor C model is an exact representation of this strength stratification since in both models the linear increase of brittle strength with depth is approximated by a layer of uniform strength. In addition, in the C model the strength falls to zero at depth while in the J model the
381 strength is discontinuous at the base of the strong layer. For the purposes of the present study, the strong layer thickness will be taken to correspond to the depth of the brittle-ductile transition. In a model lithosphere with a plastic surface layer, the layer-substrate interface corresponds to a change in power law exponent, indicating the change in the mode of deformation from faulting to ductile flow. Strength decreases rapidly with depth below the brittle-ductile transition. If estimated from this strength stratification, the strong layer thickness would not be significantly greater than the depth to the brittle-ductile transition. As shown in Fig. 4, a rift formed in a lithosphere with a plastic surface layer has a flank-toflank width = 4h. If the strong upper crustal region of the continental lithosphere is best described by a plastic theology, then a continental rift with a typically observed width in the range 35-60 km [19] requires a strong layer thickness of 9-15 km. An extending continental lithosphere in which a quartz rheology approximates flow in the crust undergoes the transition from brittle to ductile behavior at a depth of about 15-20 km depending on the geothermal gradient and strain rate (cf. [18]). This is comparable to the depth to which earthquakes are observed in the continental lithosphere and to which faulting penetrates in continental rifts as shown in seismic reflection profiles. For a strong viscous surface layer with nl = 3, the dominant wavelength, and therefore the rift zone width, depends on the ratio of layer to substrate strength as shown in Fig. 6. For a rift zone 60 km wide, typical of the well-studied East African Rift or Rhinegraben, a strong layer thickness of 15 km requires an R only slightly greater than unity. The width of the rift is controlled by the dominant wavelength as long as the initial perturbation width d < ?~d. However, if d > ?~a the width of the rift is governed by the width of the initial disturbance. If an initial perturbation is much wider than the dominant wavelength, corresponding perhaps to broad-scale lithospheric thinning, extension localizes in the thinned region and deformation within that region develops periodically at the dominant wavelength [6,20]. The vast regional extent of the Basin and Range Province ( > 103 km) may reflect an initial strong layer thickness per-
turbation much wider than the spacing of individual basins and ranges ( -~ 30 km). The depth distribution of seismicity [21], high heat flow [22], and low Pn velocities [23] in this area suggest that an upper mantle thermal anomaly of large horizontal extent could be the cause of the initial perturbation. 4.3. Rifting on Venus
A number of features on Venus with the morphology of rift zones have been recognized from Pioneer Venus radar altimetry [24]. High-resolution earth-based radar images have recently been obtained for a major rift in Beta Regio [25]. As shown in the across-strike profiles in Fig. lb, the topography consists of a central depression and bounding highs which may represent uplifted flanks. The width of the rift varies along strike but averages about 150 km from flank to flank. Two volcanic shields occur along the rift, Rhea Mons in the north and Theia Mons in the south. While volcanic construction could contribute to the flanking highs, they are best developed where the central depression is deepest and in areas not associated with the volcanoes (profiles G and H in Fig. lb). The high surface temperature of Venus, approximately 700 K, suggests that ductile deformation will occur at shallower depths than on the earth. On the basis of estimates of rock strength, surface temperature, and the spacing of features of presumed tectonic origin in the banded terrain, Solomon and Head [26] suggest that the elastic lithosphere thickness on Venus is in the range 1-10 km, and therefore too thin to explain the width of the Beta Regio rift by elastic flexure, which would require an elastic-brittle layer thickness in excess of 60 km. Schaber [24] estimated a brittle layer thickness in the range 47-69 km on the basis of the widths of other presumed rift valleys on Venus. Radar bright lineaments in Beta Regio, interpreted as faults [25], suggest the presence of a brittle surface layer. On the basis of Fig. 4, the width of the Beta Regio rift would require a plastic layer thickness of about 40 km. If the near-surface of Venus is weaker than the earth, S will be greater than that assumed in Fig. 4. The limit of large S corresponds to a rift zone width = 6h and a plastic layer thickness of 25 km. This
382
is greater than the depth to the brittle-ductile transition on Venus estimated by Solomon and Head [26] and would require a geothermal gradient less than that which they assume (20 K k m - l ) . The 10-20 km spacing of bright lineaments where they are most well-developed is comparable to the width of radar bright bands in the banded terrain. If the radar bright features are due to faulting with a spacing comparable to the brittle-elastic layer thickness, then this thickness in Beta Regio should be similar to that in areas of banded terrain and too thin to explain the rift zone width. If the plastic region of near-surface faulting is thin, then lithospheric necking may be controlled by a thicker viscous layer deforming by creep. In this case the width of the rift zone is determined by the strong layer thickness and the layer/substrate strength contrast. On the basis of Fig. 6, 10 and 20 km thick strong layers require values of R = 100 and 10, respectively, to explain the width of the Beta Regio Rift. 5. Summary Unstable extension of a lithosphere consisting of a strong surface layer overlying a weaker viscous substrate results in a pattern of deformation that is consistent with the major morphological characteristics of rift zones. The surface topography consists of a central depression, beneath which the layer thins by necking, and flanking uplifts. The rift is nucleated by introducing a small amplitude thickness perturbation at the base of the layer. For an initial perturbation narrower than the dominant wavelength, deformation concentrates in a zone of width comparable to the dominant wavelength. For an initial thickness perturbation wider than the dominant wavelength, deformation develops periodically at the dominant wavelength in the region above the perturbation. The dominant wavelength is controlled by the layer thickness and by the growth rate spectrum of extensional instability, which is a function of the layer/substrate strength contrast, the density stratification, and the stress exponents describing flow in the layer and substrate. Extension of a strong surface layer that deforms by distributed faulting, idealized by plastic behavior (n 1 = oo), results in a rift zone width approximately four times the plastic layer thickness. If the base of the
strong layer corresponds to the maximum depth of brittle deformation, determined either by experimental flow laws or the observed depth of earthquakes, then the width predicted by this model is generally consistent with that of a typical continental rift zone. Dominant wavelength and therefore rift zone width for a plastic surface layer is essentially independent of the layer/substrate strength contrast. For a strong viscous surface layer (n~ = 3), dominant wavelength varies as the strength ratio to the one-third power. Extension of a viscous surface layer always produces a rift zone that is wider than that for a plastic surface layer. The great width of rift zones on Venus, relative to the spacing of other features of presumed tectonic origin, may be explained by ductile necking of a strong viscous layer. Acknowledgements This research was supported by NASA grant NSG-7605. We thank Ray Fletcher for helpful discussions and Philip England for a constructive review of the manuscript. References 1 F.A. Vening Meinesz, Les "graben" africains r6sultat de compression ou de tension dans la croBte terrestre?, K. Belg. Kol. Inst. Bull. 21, 539-552, 1950. 2 M.E. Artemjev and E.V. Artyushkov, Structure and isostasy of the Baikal Rift and the mechanism of rifting, J. Geophys. Res. 76, 1971. 3 M.H.P. Bott, Formation of sedimentary basins of graben type by extension of continental crust, Tectonophysics 36, 77-86, 1976. 4 A.H.C. SengSr and K. Burke, Relative timing of rifting and volcanism on earth and its tectonic implications, Geophys. Res. Lett. 5, 419-421, 1978. 5 E.W. Hart, Theory of the tensile test, Act. Metall. 15, 351-355, 1967. 6 R.C. Fletcher and B. Hallet, Unstable extension of the lithosphere: a mechanical model for Basin and Range structure, J. Geophys. Res. 88, 7457-7466, 1983. 7 R.C. Fletcher, Folding of a single viscous layer: exact infinitesimal amplitude solution, Tectonophysics 39, 593-606, 1977. 8 R.B. Smith, Unified theory of the onset of folding, boudinage, and mullion structure, Geol. Soc. Am. Bull. 86, 1601-1609, 1975. 9 S.H. Emerman and D.L. Turcotte, A back-of-the-envelope approach to boudinage mechanics, Tectonophysics 110, 333-338, 1984. 10 M.A. Biot, Folding instability of a layered viscoelastic
383
11
12
13 14
15
16
17
18
19
medium under compression, Proc. R. Soc. London, Ser. A 242, 444-454, 1957. R.B. Smith, Formation of folds, boudinage, and mullions in non-Newtonian materials, Geol. Soc. Am. Bull. 88, 312-320, 1977. M.A. Biot, The influence of gravity on the folding of a layered viscoelastic medium under compression, J. Franklin Inst. 267, 211-228, 1959. D. McKenzie, Some remarks on the development of sedimentary basins, Earth Planet. Sci. Lett. 40, 25-32, 1978. H.D. Lynch, Numerical Models of the Formation of Continental Rifts by Processes of Lithospheric Necking, pp. 164-216, Ph.D. Thesis, New Mexico State University, Las Cruces, N.M., 1983. W.R. Buck, Small-scale convection induced by passive rifting: the cause for uplift of rift shoulders, submitted to Earth Planet. Sci. Lett., 1985. C.E. Keen, The dynamics of rifting: Deformation of the lithosphere by active and passive driving mechanisms, Geophys. J. R. Astron. Soc. 80, 95-120, 1985. M.S. Steckler, Uplift and extension at the Gulf of Suez: Indications of induced mantle convection, Nature 317, 135-139, 1985. W.F. Brace and D.L. Kohlstedt, Limits on lithospheric stress imposed from laboratory experiments, J. Geophys. Res. 85, 6248-6252, 1980. I.B. Ramberg and P. Morgan, Physical characteristics and
20
21
22
23
24
25
26
evolutionary trends of continental rifts, Proc. 27th Int. Geol. Congr. 7, 165-216, 1984. M.T. Zuber, E.M. Parmentier and R.C. Fletcher, Extension of continental lithosphere: A model for two scales of Basin and Range deformation, J. Geophys. Res., in press, 1985. G.P. Eaton, The Basin and Range Province: origin and tectonic significance, Annu. Rev. Earth Planet. Sci. 10, 409-440, 1982. A.H. Lachenbruch and J.H. Sass, Models of an extending lithosphere and heat flow in the Basin and Range Province, in: R.B. Smith and G.P. Eaton, eds., Geol. Soc. Am. Mem. 152, 209-250, 1978. R.B. Smith, Seismicity, crustal structure, and intraplate tectonics of the interior of the western cordillera, in: R.B. Smith and G.P. Eaton, eds., Geol. Soc. Am. Mem. 152, 111-144, 1978. G.G. Schaber, Venus: limited extension and volcanism along zones of lithospheric weakness, Geophys. Res. Lett. 9, 499-502, 1982. D.B. Campbell, J.W. Head, J.K. Harmon and A.A. Hine, Venus: volcanism and rift formation in Beta Regio, Science 226, 167-170, 1984. S.C. Solomon and J.W. Head, Venus banded terrain: tectonic models for band formation and their relationship to lithospheric thermal structure, J. Geophys. Res. 89, 6885-6897, 1984.