ABSTRACT
High-quality, normal-incidence seismic reflection data confirm that tectonic deformation in the central Indian Ocean occurs at two spatial scales: whole lithosphere folding with wavelengths varying between 100 and 300 km, and compressional reactivation of crustal faults with a characteristic spacing of c. 5 km. Faults penetrate through the crust and probably into the upper mantle. Both types of deformation are driven by regional large intraplate stresses originating from the Indo-Eurasian collision. Numerical modelling of the spatial and temporal relationships between these two modes of deformations shows that, in agreement with geophysical observations, crustal faults are reactivated first with stick-slip behaviour. Subsequent lithospheric folding does not start until horizontal loading has significantly reduced the mechanical strength of the lithosphere, as predicted by elasto-plastic buckling theory. Modelling suggests that lithospheric folding does not develop in the absence of fault reactivation. Crustal fault reactivation, therefore, appears to be a key facilitating mechanism for oceanic lithospheric buckling in the central Indian Ocean.
INTRODUCTION
Recent single-channel and multi-channel seismic reflection studies (e.g. Bull and Scrutton, 1990, 1992; Pilipenko and Sivukha, 1991; Chamot-Rooke et al., 1993) have demonstrated the occurrence of two spatial scales of active tectonic deformation in the central Indian Ocean (Fig. 1), confirming earlier findings from single-channel seismic reflection data and compilations of earthquake, heatflow, geoid and gravity data (Weissel et al., 1980; Geller et al., 1983). The tectonic deformation is characterised by long wavelength undulations of oceanic basement and overlying sediments, and superimposed small scale faulting and folding (Fig. 2). The active intraplate deformation in the Indian Ocean is thought to be driven by large regional compressional stresses (e.g. Curray and Moore, 1971; Patriat and Achache, 1984). The large stresses may originate from the intense plate boundary processes (Indo-Eurasion collision; Indonesian subduction) surrounding this intraplate section of the Indo-Australian plate (Cloetingh and Wortel, 1985, 1986), or be related with a possible (nascent) diffuse plate boundary separating the Australian plate from the Indo-Arabian plate (Petroy and Wiens, 1989; DeMets et al., 1990; Gordon et al., 1990; Royer and Chang, 1991).
Figure 1. Illuminated relief image of ERS/GEOSAT/SEASAT derived free-air gravity anomalies (Sandwell and Smith, 1992) within the central Indian Ocean, and location map of main features. Boxes denote areas with consistently oriented linear trends in the gravity anomalies coinciding with axes of undulations of oceanic basement. Short line in left box denotes position of long seismic reflection line shown in Fig. 2.
Several tectonic models have been proposed to explain the undulations
of the oceanic basement, from the interpretations of geoid and
gravity anomalies, seismic reflection and refraction studies,
and numerical and analogue modelling. Most support is given to
flexural folding/buckling of the mechanically strong upper part
of the lithosphere (Weissel et al., 1980; McAdoo and Sandwell,
1985; Karner and Weissel, 1990; Bull et al., 1992; Martinod
and Davy, 1992; Molnar et al., 1993). Other proposed models
are inverse boudinage of the oceanic crust caused by hydrodynamic
viscous flow (Zuber, 1987; Leger and Louden, 1990), whole-crustal
block faulting (Neprochnov et al., 1988), and a flexurally
folded oceanic crust decoupled from the subcrustal lithosphere
(Verzhbitsky and Lobkovsky, 1993). Theoretical studies (e.g. Timoshenko
and Woinowsky-Krieger, 1959; Martinod and Davy, 1992) show that
some kind of perturbation is required to destabilize the homogeneous
shortening of a plate under compression in order to initiate other
modes of deformation. In the previously mentioned modelling studies
of intraplate deformation, various perturbation mechanisms have
been assumed, such as line loading by the aseismic Ninety-East
Ridge (McAdoo and Sandwell, 1985), the growth of mantle instabilities
(Zuber, 1987), point loading by seamounts (Karner and Weissel,
1990), or surface sediment loading (Molnar et al., 1993).
However, regardless of the tectonic model and perturbation mechanism,
all studies concentrate on reproducing only the long wavelength
undulations of oceanic basement, neglecting the crustal faulting.
This paper presents results of a numerical analysis using a nonlinear,
elasto-plastic finite element method to analyse the temporal and
spatial relationship between the two modes of tectonic deformation
(lithospheric buckling and crustal fault reactivation) observed
in the central Indian Ocean. The main objective is to investigate
the role of crustal fault reactivation as a potential perturbation
mechanism.
INTRAPLATE DEFORMATION IN THE CENTRAL INDIAN OCEAN
In the central Indian Ocean Basin, undulations of oceanic basement
and, if present, sedimentary cover (Figs 2 and 3) trend E-W to
NE-SW, with an amplitude of 1-2 km and a wavelength that increases
from less than 100 km in the south to more than 300 km towards
the north, averaging around 200 km (Weissel et al., 1980;
Zuber, 1987; Neprochnov et al., 1988; Bull and Scrutton,
1992). The increase in wavelength is possibly related to the northward
increase in age from 40 to 70 Ma of the oceanic lithosphere, and
the associated increase in flexural rigidity (McAdoo and Sandwell,
1985). Within the Wharton Basin the undulations have similar wavelengths
and amplitudes, but strike more consistently NE-SW (Stein et
al., 1989; Pilipenko and Sivukha, 1991). In both basins, the
axes of undulations are roughly perpendicular to directions of
maximal horizontal compressional stress, as inferred from earthquake
focal mechanism studies (Bergman and Solomon, 1985; Levchenko,
1989; Petroy and Wiens, 1989) and numerical computations of the
Indo-Australian plate stress field (Cloetingh and Wortel, 1985,
1986). The basement undulations are not symmetrically shaped,
but have broad lows and pronounced crests. In the southern parts
of the area of intraplate deformation there is negligible sediment
thickness. To the north, there is a progressive infill of the
undulation lows by the sediments of the Bengal and Nicobar fans
(Fig. 1).
Geoid and free-air gravity anomalies contour maps also exhibit
undulating character (Figs 1 and 3), correlating with oceanic
basement undulations. Folding axes appear to be discontinuous
and even seem to undergo rotation across old fracture zones (Bull,
1990; Bull and Scrutton, 1992). The offset of undulations across
the fracture zones is consistent with N-S compression. Some rotation
is to be expected as the 3-dimensional intraplate stress field
is unlikely to be uniformly oriented over this wide area. This
is compatible with stress orientations derived from earthquake
focal mechanism studies (Bergman and Solomon, 1984,1985; Petroy
and Wiens, 1989). Numerical computations of the Indo-Australian
plate stress field (Cloetingh and Wortel, 1985, 1986) show how
the unique geometry of the plate boundaries and the distribution
of dynamic plate boundary processes (collision, subduction, spreading)
led to rotated directions of principal stresses with variable
magnitudes in the interior of the Indo-Australian plate.
Superimposed on the long wavelength undulations of oceanic basement
is intense short wavelength folding and faulting in the crust
and overlying sediments. Recent multi-channel seismic profiles
(Bull and Scrutton, 1990, 1992; Chamot-Rooke et al., 1993)
provide new evidence that the faulted blocks are bounded by high-angle
reverse faults penetrating through the crust and possibly into
the upper mantle (Fig. 2). Fault activity is recognised on the
seismic images by displacements of oceanic basement and sedimentary
horizons. Some unreactivated pre-existing basement faults have
also been identified, as well as some intra-crustal low-angle
faults and sub-horizontal reflectors. The new seismic data also
reveals the widespread existence of open anticlinal folds of 5-10
km wavelength within the sedimentary cover, most of them occurring
with limited faulting in the deeper parts or with faulting extending
from seafloor to basement.
Figure 2. Migrated multichannel seismic reflection profile (from Charles Darwin cruise 28-CD28; Bull and Scrutton, 1990, 1992). Location of the line is displayed in Fig. 1. Reverse faults cutting through the crust are indicated by arrows and possible Moho reflections are highlighted by dots. Vertical scale is in seconds, two-way travel time.
The fault pattern is very complex with variations in fault type, fault dip, dip direction, penetration depth, fault spacing and length, and fault activity. This is, for example, illustrated by the north-south running profile of acoustic basement (Fig. 3), where the small scale effects of varying fault throws are seen to be superimposed on the long wavelength undulations of the oceanic basement. The vast majority of high-angle faults display reverse displacements (Bull and Scrutton, 1990; Pilipenko and Sivukha, 1991). Basement faults can be divided into those that dip 40 to the south, and those that dip 35-40 to the north. If penetrating into the sedimentary cover, faults dramatically steepen (40-90), presumably due to the change in rheology (Melosh, 1990). Some of the N-dipping faults can be resolved to Moho depth (Bull and Scrutton, 1990). The S-dipping faults are less well resolved, as they dip more steeply. Both sets strike 90-100 E, roughly perpendicular to the strike of fracture zones in the area. The fault length averages 10 km, with lengths up to 40 km. Average fault spacing is 5-6 km. The fault amplitude, in terms of the vertical throw of the basement/sedimentary cover interface, averages several hundreds of meters, ranging from very small (detection limit) to over more than 1000 meters (Chamot-Rooke et al., 1993). Some of the basement faults on the seismic profiles seem to penetrate the entire crust and possibly into the subcrustal lithosphere. Bull and Scrutton (1992) suggest that these faults may nucleate in the upper mantle at the brittle-ductile transition and propagate upwards, reactivating the ridge-parallel fault fabric in the oceanic crust.
Figure 3. North-south running profiles (for location see Fig. 1) of digitised acoustic basement (solid line) and free-air gravity anomaly (dashed line), demonstrating both spatial modes of tectonic deformation: long wavelength folding and small scale faulting.
These geometrical characteristics, together with other factors,
have led to the interpretation that the faults result from the
reactivation of the oceanic structural fabric that originated
at the midocean spreading centres (Bull, 1990). Steep dipping
syn- and antithetic faults are formed in the transition of the
rift valley into the rift mountains, creating fault blocks with
widths ranging from 1 to 5 km (Macdonald, 1982; Kaz'min and Borisova,
1992). A similar complex fault fabric (in pristine state), although,
in general, with lower dips, has also been observed on seismic
images of oceanic crust in the Atlantic (White et al.,
1990).
Excellent resolution of faults is attributed to the presence
of hydrothermal alteration fronts along the fault planes (Shipboard
Scientific Party, 1990; Bull and Scrutton, 1992). The hydrothermal
activity within the crustal fault blocks is thought to be responsible
for the locally anomalous high surface heat flow measured within
the intraplate deformation area (Geller et al., 1983; Stein
and Weissel, 1990). Regionally, heat flow is normal (as predicted
by cooling plate models), indicating that there is no deep thermal
origin for the driving mechanism of tectonic deformation. Furthermore,
the presence of fluids in the fault planes may have a weakening
effect on the frictional resistance against sliding (Brace and
Kohlstedt, 1980) and facilitate fault reactivation.
Seismic stratigraphy and DSDP and ODP drilling have demonstrated
the presence of two major unconformities within the sediments
of the Bengal and Nicobar distal fans over much of the central
and northeastern part of the Indian Ocean (Cochran, 1990; Curray
and Munasinghe, 1989). An early Eocene (55 Ma) unconformity is
correlatable with the first stage of collision of the Indo-Australian
plate with the Eurasian plate. A latest Miocene (7 Ma) unconformity
has been correlated with the collision of the Indian continent
with the Eurasian plate (Peirce, 1978; Patriat and Achache, 1984;
Klootwijk et al., 1985; Cochran, 1990; Molnar et al.,
1993). The Miocene unconformity separates pre- and syn-deformational
sediments, dating the onset of the intraplate deformation at 7
Ma. There is little evidence for deformation before this time.
Sedimentary onlap patterns suggest that at 7 Ma there was a major
shortening event followed by several smaller pulses of intense
fault block rotation and fault motion (Bull and Scrutton, 1992).
The relative timing of crustal faulting versus long wavelength
folding is unclear, and requires more detailed stratigraphic analysis.
Analysis of offset sedimentary horizons shows that on several
faults displacement decreases steadily upwards: these faults grow
by upwards propagation. Other faults have a roughly constant offset
through time, and either formed rapidly at one of the deformation
pulses, or are faults on which propagation is complete.
Statistical analysis shows that there is no relationship between
magnitude of fault throw and position relative to the long wavelength
basement undulations (Bull, 1990). Chamot-Rooke et al.
(1993) report a gradual increase in fault throw from south to
north over the area of intraplate deformation in the Central Indian
Ocean Basin, likely related with the northward increase in sediment
loading by the Bengal fan. These authors find for the same area
a very long wavelength (600-700 km) trend in fault downthrow direction.
The southernmost area is affected primarily by north-dipping faults,
the central area by south-dipping faults, whereas the more extensively
deformed northernmost area has a mixture of both polarities, although
with a majority of north-dipping faults. There is, however, no
general relationship between fault downthrow direction and basement
undulations (wavelength of 200 km), suggesting a flexural origin
for the undulations. In some cases faults tend to downthrow downslope,
accentuating the crests of the undulations.
Shortening of this part of the Indo-Australian plate, subject
to large compressional stresses, is accommodated by both long
wavelength folding and short wavelength faulting and folding.
Digitization of oceanic basement from seismic profiles shows that
the long wavelength contribution is equivalent to less than 0.1%
shortening. Bull and Scrutton (1992) estimate from seismic sections
the horizontal shortening due to the reverse faulting to be 1.2
(0.4)%, corresponding to 18 (6)km of shortening over 1500 km north-south
extent of deformation. Chamot-Rooke et al. (1993) also
assess a possible contribution from very small-scale faulting.
Their estimates of the shortening range from 2.5% to maximal 4.3%
when small-scale faulting is taken into account. When distributed
over the 900 km of deformed basement observed along their seismic
profile, these estimates correspond with absolute shortening values
of 22 km and 62 km, respectively. At the same longitude, plate
tectonic reconstructions (Royer and Chang, 1991) predict a finite
shortening of 23 (45)km. Assuming that the shortening has taken
place in a steady way from the onset of deformation (7 Ma) until
present, corresponding rates of shortening are 2.5 (0.9) mm.a-1
(Bull and Scrutton, 1992), 3 mm.a-1 (Royer and Chang, 1991), and
6 (3) mm.a-1 (Chamot-Rooke et al., 1993). This last rate
is consistent with the 6 to 7 mm.a-1 rate inferred from plate
kinematic models (Gordon et al., 1990). The pulselike character
of shortening through time, as suggested by the deformation patterns
in the sedimentary cover on the seismic profiles, may have led
to higher shortening rates at certain times.
Active intraplate deformation is still taking place in the central
Indian Ocean, as evidenced by the present-day occurrences of large
intraplate earthquakes (Bergman and Solomon, 1980, 1985; Levchenko,
1989; Petroy and Wiens, 1989) and reported substantial micro-seismicity
in this area (Levchenko and Ostrovsky, 1993).
A FINITE ELEMENT MODEL
In order to investigate the temporal and spatial relationship between the two modes of tectonic deformation, numerical analyses are carried out using the finite element technique. The numerical analysis assumes a 2-dimensional plane-strain approach, which is validated by the east-west continuity of the axes of the undulations of oceanic basement (Weissel et al., 1980; Bull and Scrutton, 1992). A finite element model (Fig. 4a) is constructed comprising a 2000 km long, N-S section through the central Indian Ocean part of the Indo-Australian plate. To correctly compute displacements, finite elements in the faulted upper part of the plate must be small enough to allow deformation within the faulted blocks, and are therefore sized 1x1 km triangular in the uppermost 10 km, increasing in size to 1x5 km rectangular in the lower parts (Fig. 4b). For the full mesh with dimension 2000x40 km this becomes, including spring elements, 60000 finite elements, connected by more than 42000 nodal points, and with a total of more than 80000 degrees of freedom. Correct behaviour of the elements is ensured by using higher-order elements, in this case linear triangulars and quadrilaterals which have been expanded quadratically using "incompatible modes of deformation" (e.g. Beekman, 1994).
Figure 4. (a) Geometry and boundary conditions of the mechanical model. (b) Enlarged part of the finite element mesh showing geometry of the pre-existing crustal faults. (c) Geotherm (dashed) and yield-strength envelopes (continuous) for a 50 Ma oceanic lithosphere. Assuming hydrostatic pore pressure reduces the strength of the lithosphere.
The complex fault geometry in the oceanic crust is approximated
by including a simplified set of faults in the upper part of the
plate model, which dip 45 and have a penetration depth and fault
spacing of 5 km (average values). For the computations the finite
element code TECTON is used (Melosh and Raefsky, 1981), which
contains the slippery node technique to model fault behaviour
(Melosh and Williams, 1989), and which is modified to also account
for friction along fault planes and for elasto-plastic deformation
(Beekman, 1994).
The model boundaries are allowed to move freely in both the horizontal
and vertical direction, except for the vertical outer edges which
are restricted to move horizontally. The vertical buoyancy of
the asthenosphere, which can be modelled as a fluid on geological
timescales and which counteracts vertical motions of the overlying
lithosphere, is included by attaching an array of linear elastic
Hookean springs to the base of the plate. The springs will exert
forces on the base of the plate, proportional to, but directed
opposite to, any vertical deflection of the base of the plate.
Horizontal boundary conditions at the N and S edge imitate the
effects of the India-Himalayan collision and resistance associated
with subduction along the Banda and Sunda arcs, and mid-ocean
ridge push, respectively (Cloetingh and Wortel, 1985, 1986). The
plate is shortened with a velocity of 6.3 mm.a-1, according to
the currently most accurate available shortening estimate for
the area of deformation (Chamot-Rooke et al., 1993). For
a 2000 km long plate, this velocity is equivalent to a strain-rate
of c. 10-16 s-1, which is a value
characteristic for the slow deformation
of oceanic lithosphere (Carter and Tsenn, 1987). Furthermore,
the model is loaded with gravitational body forces.
The Indian oceanic lithosphere comprises a 6 km thick basaltic
crust, overlying an olivine upper mantle. The age of the lithosphere
in the areas of intraplate deformation in the Indian Ocean varies
from 40 Ma in the south to 70 Ma in the north. An average age
of 50 Ma is assumed for the entire plate model. The temperatures
applied to the model are computed using this model geometry for
which an initial steady-state temperature solution can be derived
from the heat-conduction equation. Subsequently, the structure
is allowed to cool through conduction of heat. Transient temperature
fields are computed by solving the heat-conduction equation using
a finite difference scheme (3rd order Runge-Kutta). Thermal properties
are listed in Table 1. The assumed thermal boundary conditions
are a constant surface temperature of 0 oC, and an initial surface
heatflux of 80 mW.m-2, which is a characteristic value
for young oceanic lithosphere (Sclater et al., 1980). The
geotherm after 50 Ma of conductive cooling within the multi-layered
oceanic lithosphere is computed for the adopted thermal initial
and boundary conditions and plotted in Fig. 4c (dashed line).
The rheological behaviour of the rocks which constitute the
lithosphere is derived from laboratory experiments, and subsequently
extrapolated to geologically relevant timescales (assuming that
the empirically derived rheological laws remain valid). At low
confining pressures and low temperatures, brittle failure is predominant.
The corresponding frictional yield criterion is given by Byerlee's
law (Byerlee, 1978), which can be rewritten in terms of effective
principal stress difference, lithostatic overburden pressure,
and pore fluid pressure (Ranalli, 1987):
![]() | (1) |
and where sigma1 and sigma3 are the maximum and minimum principal stress, respectively, rho is density, g is gravitational acceleration, z is depth, and where
![]() |
Steady-state creep of a wide variety of rocks is empirically described by a ductile flow law which relates the critical principal stress difference necessary to maintain a given steady-state strain-rate of deformation to a power of the strain-rate, and which, therefore, is called power-law creep (Kirby, 1983):
![]() | (2) |
where AP, N, and EP are
empirically determined
material "constants", assumed not to vary with stress
and (p,T)-conditions (Ranalli, 1987), and where is the
strain-rate. Micro-physically this type of creep mainly involves
dislocation climb. The ductile flow law (eqn. 2) shows that the
variation of the creep stress with depth is strongly controlled
by rock-type and temperature. A steady-state bulk strain-rate
of 10-16 s-1 is assumed, induced from the
applied velocity boundary condition.
Above a critical stress level (200 MPa), power law creep breaks
down into high-stress, low-temperature plasticity (dislocation
glide), which is phenomenologically described by a Dorn law (Goetze
and Evans, 1979; Tsenn and Carter, 1987):
![]() | (3) |
with AD, ED, and D material
properties, assumed, as before, to be constant.
A depth-varying rheology has been incorporated in the finite
element model with elasto-brittle deformation in the upper parts
of the oceanic lithosphere and ductile deformation in the lower
parts. The associated yield envelope has been constructed using
the rheological equation (1) for brittle failure, and equations
(2) and (3) for creep flow (at a constant strain-rate of 10-16
s-1) of rocks, respectively, and assuming a basaltic crust and
a olivine subcrustal lithosphere. Adopted material properties
and physical parameters are listed in Table 1. The empirical rheological
equations and laboratory derived material properties are used
under the assumption they are also valid when extrapolated to
geological strain-rates and macroscopic scales. The brittle-ductile
transition lies at a depth of 24 km. The depth at which the yield
envelope has reduced to an arbitrary small stress level in the
order of 10 MPa (Ranalli, 1987) indicates the base of the mechanically
strong upper part of the lithosphere, and is used in the modelling
to control the depth of the lower boundary of the numerical plate
model, here selected at 40 km depth.
DEFINITION | SYMBOL | UNITS | VALUE | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
gravitational acceleration | g | m.s-2 | 9.81 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
universal gas constant | R | J.mole-1.K-1 | 8.314 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
seawater density | rhow | kg.m-3 | 1030 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
sediment density | rhos | kg.m-3 | 2200 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Moho depth | ZM | km | 6 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
mantle melt temperature | Tm | oC | 1300 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
mantle basal heat flux | qb | mW.m-2 | 23 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
static friction coefficient | fs | -- | 0.6 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
bulk strain-rate of deformation | epsilon | s-1 | 10-16 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
CRUST | MANTLE | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
petrology | --- | --- | basalt | olivine (dry)2 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
density | rho | kg.m-3 | 2950 | 3300 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Young's modulus | E | GPa | 90 | 70 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Poisson's ratio | nu | --- | 0.25 | 0.25 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
thermal conductivity | k | W.m-1.K-1 | 3.1 | 3.5 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
specific heat | cP | J.kg-1.K-1 | 1050 | 1050 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
power law exponent | N | --- | 3.05 | 3.0 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
power law activation energy | EP | kJ.mole-1 | 276 | 510 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
power law strain rate | epsilonP | Pa-N.s-1 | 6.31x10-20 | 7.0x10-14 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
dorn law activation energy | ED | kJ.mole-1 | --- | 535 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
dorn law strain rate | epsilonD | s-1 | --- | 5.7x1011 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
dorn law stress | sigmaD | GPa | --- | 8.5 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
FAILURE/CREEP EQUATIONS | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
brittle failure | sigmabrittle=alpha.rho.g.z(1-lambda) | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
power law creep | sigmacreep=(epsilon/AP)1/Nexp[EP/NRT] | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
dorn law creep | sigmacreep=sigmaD{1-[1-(RT/ED)ln(epsilon/AD)]1/2} |
Based on the presence of hydrothermal fluid convection in the
upper part of the oceanic lithosphere (Shipboard Scientific Party,
1992) a hydrostatic pore fluid pressure is assumed, which significantly
weakens the brittle strength of the lithosphere (equation 1: parameter).
Hydrothermal fluid alteration fronts seem to coincide with
pre-existing fault planes in the crust, and are thought to lower
the frictional resistance against renewed sliding along the fault
planes (Brace and Kohlstedt, 1980).
Brittle deformation in the depth range over which the pre-existing
faults extend is assumed to be carried out explicitly by these
predefined faults. Therefore, no (brittle) yield stress is assigned
to the finite elements in this upper part of the plate. Over the
remaining depth range down to the brittle-ductile transition there
is no information on whether brittle deformation occurs as fracturing
of rocks or as fault reactivation. What is known, as evidenced
by the existence of deep earthquakes, is that it does occur. Below
the resolution depth of the seismic profiles, which is approximately
Moho depth, there is no geometrical information of faults (if
existing). Therefore, no faults are incorporated at depths exceeding
Moho depth. The elements in this depth range are nevertheless
characterised by a brittle yield stress assigned on the basis
of the empirical brittle deformation equation (1), although any
associated brittle failure is modelled numerically as plastic
flow (Beekman, 1994).
RESULTS
At time zero the plate undergoes an instantaneous displacement
downwards in response to the gravitational loading of the model
(see, for example, Fig. 6a). This is followed by a slow linear
uplift of the plate surface in response to continuing horizontal
shortening. After several million years of shortening (0.3% of
shortening occurs in a period of 1 Ma), the models start to develop
a folding phase, as illustrated by the evolution of the surface
deflection of the plate (Fig. 5a). Fourier spectral analysis of
the calculated surface folding demonstrates a dominant wavelength
between 200 and 350 km (Fig. 5b).
The folding evolves nonlinearly through time: after several Ma,
surface deflections (Fig. 6a) begin to deviate from the homogeneous
uplift trend and accelerate upward or downward depending on their
position inside a fold (crest or trough). The onset of the deviations
corresponds to the onset of lithospheric buckling.
The pre-defined faults in the models exhibit typical stick-slip
behaviour through time (Fig. 6b). Jumps in fault slip (at the
surface of the plate) are in the order of several meters, accumulating
to throws of more than 100 m (Fig. 6b). Comparison of Figs 6a
and 6b shows that fault activity starts before the onset of the
large scale folding. After long wavelength folding has started
to develop, fault activity becomes more complex with some faults
exhibiting accelerated reversed faulting, some faults stay locked,
and some faults even show normal faulting.
The spatial relationship between long wavelength folding and
intensity of faulting is demonstrated by plotting the surface
deflection of the plate in combination with the course of the
amount of surface fault slip over the plate, both after 8 Ma of
shortening (Fig. 7a). In general, maxima in fault throw coincide
with fold lows, whereas minima in fault throw coincide with fold
highs. This is explained by analysing the fault slip at the plate
surface for all faults within the latitudinal distance running
from 900 km to 1100 km, which roughly comprises one long-wavelength
fold. Figure 7b depicts the surface fault slip for all faults
located between 900 km and 1000 km, roughly coinciding with a
fold low, and Figure 7c shows the surface fault slip for faults
between 1000 km and 1100 km, coinciding with a fold high. All
"low positioned" faults exhibit accelerating reverse
faulting, thus producing the maxima in fault slip, whereas the
"high positioned" faults show fault locking or even
an inversion to normal faulting, thus leading to the minima in
surface fault slip.
Figure 5. (a) Vertical deflection of the plate model surface through time (1 Ma equals 0.3% shortening). Wavelength of developing folding is in agreement with observations. (b) Normalised Fourier power spectra of the plate surface at several times of loading.
Fig. 8a depicts vertical effective stress profiles inside the plate (at the position denoted by the arrow in Fig. 5a) at several stages of compression. The adopted yield stress envelope limits the effective stresses in the upper brittle and lower ductile part of the lithosphere. In the faulted upper part of the plate, stress profiles develop that have a steeper slope than the theoretical brittle yield envelope. Furthermore, high effective stresses appear at a depth of 5 km (Fig. 8a), which is the depth at which the predefined faults in the model end. The steeper slope, equivalent to a higher "fault reactivation" stress, arises from pre-setting the dip angle of the predefined fault planes. This results, in general, in an angle between fault planes and maximal principal stress that is not the optimal angle for faulting (according to a Mohr diagram) (Jaeger and Cook, 1976). More effective stress is required for reactivation than in the case of an optimal angle on which the yield stress value is based.
Figure 6. (a) Vertical deflection through time of distinct points at the surface of a horizontally compressed oceanic lithospheric plate with pre-existing faults in its crust. Curves are for surface points separated by a horizontal distance of 50 km. The nonlinear time response is typical for a buckling mode of plate deformation. Inset diagram illustrates behaviour of a plate with locked faults: no buckling. (b) Amount of net slip through time at the surface appearance of a fault. Curves are for every tenth fault (horizontal separation of 50 km). Reverse slip is positive, normal slip is negative. Note the typical stick-slip response. Comparison of Figs 6a and 6b indicates that crustal fault reactivation occurs before lithospheric buckling starts.
Figure 7. (a) Deflection of the surface of the plate model (continuous line), and the distribution of surface fault slip over the plate (dashed line), both after 8 Ma (=2.4 %) of shortening. (b) Fault slip at the plate surface versus time for all faults located in the long wavelength fold low between 900 km and 1000 km, and (c) for all faults located in the long wavelength fold high between 1000 km and 1100 km.
Fig. 8a depicts vertical effective stress profiles inside the
plate (at the position denoted by the arrow in Fig. 5a) at several
stages of compression. The adopted yield stress envelope limits
the effective stresses in the upper brittle and lower ductile
part of the lithosphere. In the faulted upper part of the plate,
stress profiles develop that have a steeper slope than the theoretical
brittle yield envelope. Furthermore, high effective stresses appear
at a depth of 5 km (Fig. 8a), which is the depth at which the
predefined faults in the model end. The steeper slope, equivalent
to a higher "fault reactivation" stress, arises from
pre-setting the dip angle of the predefined fault planes. This
results, in general, in an angle between fault planes and maximal
principal stress that is not the optimal angle for faulting (according
to a Mohr diagram) (Jaeger and Cook, 1976). More effective stress
is required for reactivation than in the case of an optimal angle
on which the yield stress value is based.
Immediately below the layer which contains the faults effective
stresses reach the yield stress limit over an increasingly larger
depth range. At depths where this occurs, rocks will deform in
a brittle way. Stress also increasingly reaches the yield limit
in the lowermost part of the plate model, where rock deformation
will occur by ductile flow.
Also illustrated by Fig. 8a is the decrease in thickness of the
strong elastic core of the oceanic lithosphere and the associated
saturation of the yield envelope, until after c. 8 Ma the
elastic core has almost vanished and is not able to mechanically
sustain the stresses anymore. The lithosphere starts to fail completely,
which will result in extensive deformation, as is numerically
indicated by the acceleration of the surface deformation of the
plate model and by the increase in magnitude of fault slip (Figs
6a and 6b).
The stick-slip behaviour of faults is also reflected in the temporal
fluctuations of effective stress within crustal finite elements
(Fig. 8c). Curves are for the first 5 upper-triangular elements
in the column below the location denoted by the arrow in Fig.
5a. The "lowest stress" curve belongs to shallowest
element (depth range 0-1 km), whereas the "highest stress"
curve belongs to the deepest considered element (depth range 4-5
km). Repeated periods of stress buildup are followed by unlocking
of the fault and displacement of the fault blocks resulting in
stress relaxation. Stress drops associated with the stick-slip
fault activity are in the order of a few tens of megapascals.
Comparison of stress drops from the same event yields information
on propagation direction of fault reactivation. Although not always
clearly distinguishable it is possible to recognise downwards
propagating stress drop fronts (for example, between 3 and 3.5
Ma) as well as upwards propagating stress drop fronts (e.g. just
after 6 Ma), while the event after 8 Ma seems to reactivate the
entire fault almost instantaneously. The upwards propagating event,
occurring slighly after 6 Ma, starts in the middle of the crust
(not originating from a deeper level).
An additional numerical experiment has been undertaken to investigate
whether it is indeed the reactivation of the faults, and associated
local changes in stress and strain, that ultimately cause the
large scale buckling of the plate. To ensure that the faults do
not reactivate they have been assigned a very high friction coefficient.
Furthermore, the elements in the "faulted" upper part
of the plate have now a yield stress, according to the regular
yield envelope (Fig. 4c), to avoid this layer behaving as an elastic
layer. Thus, failure is still included in this layer with a brittle
yield stress derived from the brittle failure equation (1), but
the associated deformation is now simulated numerically by plastic
flow. The results show that in this case even after 10 Ma of convergence
(equivalent to a horizontal shortening of 3%), the initially flat
plate is still flat and has not developed a folding mode of deformation.
Deflections and power spectra of the plate surface at several
times are included in Figs 5a and 5b, but no harmonic deviations
are evident. This is also demonstrated by Figs 6a and 6b, where
the temporal course of surface point deflection and surface fault
slip are shown in the inset diagrams. The slope of the vertical
deflection curves in Fig. 6a indicates the homogeneous thickening
of the plate in response to the horizontal shortening. Folding/buckling
is absent, even when the elastic core of the oceanic lithosphere
has almost vanished. This demonstrates that it is the crustal
fault reactivation and associated changes in stress and strain,
that in these finite element models acts as the perturbing mechanism
that facilitates the large scale mode of tectonic deformation.
Figure 8. (a) Vertical stress profiles at several stages of compression for the position marked by arrow in Fig. 5a. Solid line is the limiting yield envelope. (b) Same as (a), but now for a plate with locked faults. (c) Effective stress through time in the uppermost five "crustal" upper-triangular finite elements in the column at the location denoted by the arrow in Fig. 5a. The "lowest stress" curve belongs to shallowest element (depth range 0-1 km), whereas the "highest stress" curve belongs to the deepest considered element (depth range 4-5 km). Fault response is typical stick-slip with cycles of stress buildup followed by stress relaxation as the fault becomes active. Lining up stress drops belonging to the same event, shows that activation fronts propagate downwards as well as upwards along a fault plane. Both the stick-slip behaviour and up- and downward propagating reactivation of crustal faults are confirmed by seismic stratigraphic observations.
DISCUSSION
Although quite advanced in the treatment of frictional faulting
and the incorporation of a nonlinear depth-varying rheology, the
numerical models that have been used are relative simple when
considering the complex features and processes within the oceanic
lithosphere of the central Indian Ocean. Nevertheless, the models
successfully predict stress and strain patterns associated with
the long wavelength folding and small scale faulting occurring
in this intraplate part of the Indo-Australian plate, which is
subjected to considerable compressional stress loading and associated
tectonic shortening.
The dominant wavelength of large scale folding that is developed
in the models lies well within the range of observed wavelengths
in the central Indian Ocean (Zuber, 1987). It is slightly larger
than the average observed wavelength, but this can be attributed
to the simplified character of the model, not accounting for weakening
effects introduced by intra-lithospheric inhomogeneities. The
yield envelope, which exerts a major control on the flexural rigidity
of the lithosphere, is also a theoretical idealisation, probably
resulting in a mechanically too strong lithosphere. The sharp
brittle-ductile transition, for example, will probably be transitional,
as argued by Carter and Tsenn (1987). Furthermore, recent rock
mechanic experiments (Ord and Hobbs, 1989) indicate that the linear
brittle branch of the yield envelope may even break down, and
maximum sustainable stresses may be as low as 300 MPa, substantially
reducing the maximum strength of the lithosphere.
The nonlinear behaviour through time is consistent with the "buckling"
response of a horizontally compressed plate, characterised by
rapidly accelerating vertical displacements prior to full lithospheric
failure (Cloetingh et al., 1989; Stephenson et al.,
1990; Stephenson and Cloetingh, 1991). Analysis of vertical effective
stress profiles shows that buckling starts when the remaining
elastic core of the oceanic lithosphere has almost vanished, as
predicted by elasto-plastic plate deformation theory (McAdoo and
Sandwell, 1985).
Fault activity starts before the onset of the large scale folding.
The increasingly complex activity of the pre-defined faults after
long wavelength deformation has started to develop, seems to be
caused by an interaction of developing bending stresses with the
tectonic compressional stress. Support for this is given by the
relationship of fault type and folding magnitude with fault position
within a long wavelength fold. Bending would produce additional
compressive stress in the downflexed areas and tensile stress
in the upflexed areas. These compressive or tensile bending will,
respectively, assist or counteract the tectonic stresses. This
leads to a decrease in effective stress in upflexed regions, accompanied
by permanent locking of faults, or even leads to net tensile effective
stresses and normal fault reactivation, as indeed demonstrated
by the model results. Statistical analysis of bathymetric and
seismic data is required to see if the effects of bending are
present.
The pre-defined faults in the models exhibit typical stick-slip
behaviour through time. Jumps in fault slip at the plate surface
are in the order of several meters, accumulating to throws at
the plate's surface of more than 100 m. These values fit with
measured fault slips on seismic profiles (Bull and Scrutton, 1992;
Chamot-Rooke et al., 1993). The stick-slip fault activity
is in agreement with observations of substantial microseismicity
in the crust of the central Indian Ocean (Levchenko and Ostrovsky,
1993). Deep earthquakes, which have been observed in this intraplate
area (Bergman and Solomon, 1980, 1985), are also correctly predicted
by the models with stress within the plate reaching the brittle
failure limit over considerable depth ranges down to the brittle-ductile
transition (Govers et al., 1992). Reactivation fronts along
a single fault can propagate downwards as well as upwards, or
reactivate the entire fault instantaneously. Upwards propagating
faults and possibly rapidly formed faults have been identified
on the seismic stratigraphic record (Bull and Scrutton, 1992).
Further analysis is required to check if there are relationships
of type of fault reactivation with the position of the fault concerned
relative to the long wavelength folds.
The numerical models develop high effective stresses at the tips
of the faults. These stress concentrations can be interpreted
physically as being responsible for propagation of faults. Such
stress concentrations at fault tips are predicted theoretically
by the Griffith Crack Theory (Scholz, 1990), and are consistent
with rock-mechanics experiments showing that stress concentrations
do occur at fault tips (Jaeger and Cook, 1976).
Modelling also shows that no long wavelength mode of tectonic
deformation develops when there is no fault activity. This demonstrates
that in this study the crustal fault reactivation and associated
changes in stress and strain act as the essential perturbing mechanism
that facilitates the large scale mode of tectonic deformation.
The complex fault fabric in the Indian Ocean crust has been approximated
by a simplified set of faults with a geometry based on average
values. Additional numerical analyses have to be carried out to
investigate if variations in one or more of the fault parameters
(such as dip, dip direction, penetration depth) have a major influence
on the long wavelength deformation of the Indian Ocean lithosphere.
The successful predictions of the numerical model with the "average"
fault set do suggest, however, that it is merely the reactivation
of faults that is of essential importance rather than the geometrical
features of the faults.
Another model simplification which may affect the model results
is the assumption of a constant lithospheric thermal age (and
thus lithospheric thickness) for the entire model, whereas the
the thermal age of the area of interest ranges from approximately
40 Ma in the south to 70 Ma in the north. The increase in age
is probably to gradual to initiate the intraplate deformation.
It may, however, affect the wavelength of the large scale folding,
as with age also the flexural rigidity of the lithosphere increases,
which will lead to a larger wavelength of deformation (McAdoo
and Sandwell, 1985). Zuber (1987) indeed observed a gradual northward
increase in the wavelength of lithospheric folding in the central
Indian Ocean.
CONCLUSIONS
Substantial compressional stresses in the central Indian Ocean
reactivate pre-existing crustal faults, and, in a later phase,
initiate folding and buckling of the entire oceanic lithosphere.
Stress and strain patterns, spatial as well as temporal, are successfully
predicted by the numerical models, and are in agreement with geophysical
observations.
The modelling also demonstrates that in the absence of fault
reactivation no lithospheric folding develops. Crustal fault reactivation,
therefore, appears to be essential to facilitate long-wavelength
oceanic lithospheric buckling in the central Indian Ocean.
ACKNOWLEDGEMENTS
We thank G. Ranalli and R.A. Stephenson for their constructive reviews of an early draft of this paper. Comments by F. Nieuwland and an anonymous reviewer are also gratefully acknowledged. This is NSG publication 950408.
REFERENCES
Institute of Earth Sciences, Vrije Universiteit
De Boelelaan 1085
1081 HV Amsterdam
The Netherlands
Phone: +31-20-4449801
Fax: +31-20-4449943
Last Modified: August 22, 1996.