Tectonophysics 484 (2010) 181–192

Contents lists available at ScienceDirect

Tectonophysics j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / t e c t o

Structural evolution of a three-dimensional, finite-width crustal wedge Jean Braun ⁎, Philippe Yamato 1 Géosciences Rennes, Université de Rennes 1, CNRS, Rennes Cedex CS 35042, France

a r t i c l e

i n f o

Article history: Received 20 February 2009 Received in revised form 23 June 2009 Accepted 31 August 2009 Available online 10 September 2009 Keywords: Wedge mechanics Three-dimensional Basal discontinuity Structural geology Surface processes Numerical model

a b s t r a c t We present the results of three-dimensional numerical experiments designed to study the response of a layer of crustal material subjected to convergence through an imposed basal velocity discontinuity and to surface processes (erosion/sedimentation). We focus on the three-dimensional response of the system arising from the finite width of the imposed velocity discontinuity. In particular, we describe the complex structures that develop around the wedge and their interactions with the loading/unloading produced by the surface processes. We show that the pro- and retro-shear zones that develop in a doubly-vergent twodimensional wedge curve around the end of the velocity line-discontinuity to merge into the strike-slip structure that naturally develop, i.e. as a consequence of the imposed boundary conditions, along the edge of the wedge. Along the retro-shear zone the stress orientation rotates along a vertical axis, which implies that the retro-shear zone is a pure thrust along all of its curved length, whereas, along the pro-shear zone stresses rotate along a horizontal axis, which, in turn, implies that the pro-shear zone progressively evolves towards an oblique thrust in its curved section. Furthermore, the outward motion (i.e. perpendicular to the direction of imposed shortening) along the curved section of the retro-shear zone is accommodated by oblique extension along a secondary, kinked structure antithetic to the retro-shear zone. We also show the complex evolution of the wedge when ductile flow and ductile strain softening is activated by increasing the imposed basal temperature. In this situation, the wedge is broader as structures develop at finite distances on either side of the line-discontinuity and its dynamics resembles more that of a ‘vise-like’ orogen. At the surface, a flat plateau forms that accumulates sediment from the surrounding actively deforming mountain ranges until a channel breaks through one of the sides and flushes the inward-draining basin of its sedimentary content. The complex numerical models provide a detailed insight into the geometry and complex threedimensional structures within and at the margins of an orogen, and how they can be affected by erosional processes at the surface. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Small orogens, such as the Southern Alps in New Zealand or Taiwan, are often regarded as doubly-vergent critical wedges that form above an imposed velocity discontinuity (or S-point) representing subduction of the underlying mantle lithosphere (Malavielle, 1984; Koons, 1990; Willett et al., 1993). As demonstrated in many analytical, analogue and numerical models, the behavior of a critical wedge is dictated by its internal rheological properties (i.e., internal and basal angles of friction and cohesion) as well as the nature and efficiency of the processes (erosion and deposition) acting on its surface (Dahlen, 1984; Willett et al., 1993). A wide variety of behaviors has been shown by varying the rheological and erosional parameters of wedge systems

(Beaumont et al., 1992; Willett, 1999a,b; Whipple and Meade, 2004; Konstantinovskaia and Malavieille, 2005). However, and except for some very few exceptions (Braun, 1993; Koons, 1994; Braun and Beaumont, 1995; Gerbault et al., 2002), these efforts have been limited to two-dimensional studies, mostly due to the high computational cost of three-dimensional numerical models. None have explored the dynamics of a finite-width orogenic wedge. Here we present results from a new fully-coupled thermo-mechanical and surface process three-dimensional model describing the behavior of a frictional and viscous critical wedge driven by a finite-width (or S-line) velocity discontinuity along its base.

2. Numerical model description and experimental setup

⁎ Corresponding author. Fax: +33 2 2323 6780. E-mail address: [email protected] (J. Braun). 1 Now at the Geological Institute and Department of Earth Sciences, ETH-Zürich, Switzerland. 0040-1951/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2009.08.032

A full description of the numerical methods used in the model, called DOUAR, can be found in Braun et al. (2008). Here, we will limit ourselves to the description of its main features. The model is centered on a finite element solution of the quasi-static force balance equations in three dimensions using an octree division of space that is dynamically adapted

182

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

to the evolving solution. An octree division of the unit cube is an efficient, regular, yet non uniform way to define a finite element grid using 8 noded, tri-linear elements. It suffers however from nodal incompatibilities along the interfaces between adjacent elements of different sizes. These incompatibilities are solved by introducing linear constraints between the nodal degrees of freedom that belong to both elements and those which don't. Rocks are assumed to behave like an incompressible visco-plastic fluid that obeys the Stokes flow law: ∇⋅μð∇v + ∇vT Þ  ∇p = ρg ∇⋅v = 0

ð1Þ

where v is the velocity field, μ is the viscosity, p the pressure, g the acceleration due to gravity and ρ the density. During deformation, rocks are also assumed to obey Mohr–Coulomb failure criterion: τ = c0  σn tan ϕ

ð2Þ

where τ is the magnitude of the shearing stress and σn is the normal stress, c0 is cohesion and ϕ is the angle of friction, and to flow according to a power-law thermally-activated creep law: 1 = n1 Q = nRT

μ = μ0 ε˙

In an elasto-plastic material, that undergoes an almost infinitesimal dilation during the early stages of (elastic) deformation, and as commonly observed in rocks, faults (brittle shear bands) form in a plane containing σ2, the intermediary principal stress, and at 45° ± (ϕ/4 + ψ/4) to the direction of σ1, where ϕ is the angle of friction and ψ is the dilatancy angle, commonly assumed to be identical to the friction angle. Due to the incompressibility condition, in our model, shear bands that form in either of the two rheologies (brittle and viscous), form in a plane that is at 45° angle to the direction of σ1. This difference between the model behavior and that of the natural system must be kept in mind during the interpretation of the results shown here. Finite strain, ε is computed and stored on a three-dimensional cloud of points of self-adapting density and is used to include the effects of strain weakening which is parameterized here by progressive reduction in the internal angle of friction, ϕ, or in the intrinsic viscosity, μ0, with deformation:

ð3Þ

e

where μ is the viscosity (in Pa s), μ0 the intrinsic viscosity, ε̇ the second invariant of the strain rate tensor, n the power-law exponent, Q the activation energy, R the gas constant and T the absolute temperature (values and units are presented in Table 1). We neglect the infinitesimal elastic contribution to the deformation.

Table 1 Numerical model parameters. Parameter

Symbol

Value (units)

Model width Initial crustal layer height Imposed velocity Basal temperature (Exp.1 and 2) Basal temperature (Exp.3) Thermal diffusivity Heat production Cohesion Initial friction angle Reduced friction Frictional weakening strain range Activation energy Power-law exponent Initial pre-exponential factor Reduced pre-exponential factor Viscous weakening critical strain Viscous weakening strain range Rock density Tectonic time step length Surface processes time step length Hillslope diffusivity (Exp.1) Hillslope diffusivity (Exp.2) Hillslope diffusivity (Exp.3) Fluvial erosion/transport coefficient (Exp.1) Fluvial erosion/transport coefficient (Exp.2) Fluvial erosion/transport coefficient (Exp.3) Bedrock erodibility length scale Channel width

L h v0 T0 T0 k/ρc H c0 ϕi ϕf Δεϕ Q n μ0,i μ0,f εμi Δεμ ρ Δtt Δts Kd Kd Kd Kf

100 km 10 km 10 km Myr− 1 250 °C 400 °C 25 km2 Myr− 1 0 20 MPa 15° 5° 100% 223 × 103 J 4 2.92 × 106 Pa s1/4 2.92 × 105 Pa s1/4 20% 5% 2800 kg m− 3 104 yr 10 yr 32 × 10− 2 m2 yr− 1 0 m2 yr− 1 16 × 10− 2 m2 yr− 1 4 × 10− 2 m yr− 1

Kf

8 × 10− 2 m yr− 1

Kf

2 × 10− 2 m yr− 1

lBR wc

100 km 100 m

Reference

ϕ = ϕi +

Landscape evolution parameters are chosen within ranges used in other numerical studies: van der Beek and Braun (1998) and Kooi and Beaumont (1996). The exact value of these parameters is poorly constrained. It is the variation between the two first experiments that is of importance here.

Δεϕ

ðϕf  ϕi Þ

ð4Þ

where ϕi and ϕf are the initial and final values for the angle of friction (i.e. before and after strain softening), respectively, and εϕi and Δεϕ are the deformation level at which brittle softening is initiated and the deformation increment over which brittle softening takes place, respectively, and: μ

μ0 = μ0;i +

0:5 + arctan

ε  εi Δεμ

!

=

! π ðμ0;f  μ0;i Þ

ð5Þ

where μ0,i and μ0,f are the initial and final values for the intrinsic viscosity (i.e. before and after strain softening), respectively, and εμi and Δεμi are the deformation level at which viscous softening is initiated and the deformation increment over which viscous softening takes place, respectively. The temperature, T, is computed by including the effects of conduction and advection, according to: ρc

Tullis (2002) idem idem

ε  εϕi

  ∂T + v⋅∇T = ∇⋅k∇T + ρH ∂t

ð6Þ

where k is the thermal conductivity, ρ is density, c is heat capacity and H is heat production per unit mass, and is used to compute the viscosity. We track the free upper surface as well as internal interfaces defining regions of varying mechanical properties by advecting a selfadapting set of points and normals attached to each interface. At every time step, a Delaunay triangulation is computed among the points on each interface using a local metrics (or measure of distance) that takes into account the surface curvature. This triangulation is then used to compute a level set function on the nodes of the finite element mesh (the octree). The level set function is defined as the signed distance to the interface such that its zero iso-surface corresponds to the interface. The level set function is later used to localize any geometrical point, X, with respect to the interfaces by interpolation of the level set function from the node of the octree to X and inspecting its sign. For example, the use of level set functions provides an efficient way to determine the position of any integration point with respect to each of the potentially numerous material interfaces during the computation of the finite element matrix even in situation of complex intersections among the interfaces in a single element. Like all other surfaces, the upper surface is advected by tectonic transport, i.e. using the velocity field derived from the solution of the Stokes equation (Eq. (2)).The advection of the surface nodes may lead to bow-tying of the Delaunay triangles used to define the level set

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

183

function. To prevent this undesirable result, we apply a small correction to the position of the nodes by moving each node towards the mean position (or center of mass) of its neighboring nodes (i.e. those to which it is connected by the Delaunay triangulation). We have tested that this slight perturbation of the nodes' position does not affect the solution because it perturbs the geometry of the surface by an amount that is almost negligible with respect to the perturbation brought upon by the surface processes. At every time step, the geometry of the free upper surface is also subjected to a surface processes model (CASCADE, Braun and Sambridge (1997)) that includes the effects of fluvial and hillslope processes. In CASCADE, we assume that river channels have a carrying capacity (volume that can be carried per unit time), qc, that is proportional to local slope, S, and drainage area, A, a proxy for local discharge (Kooi and Beaumont, 1994): qc = Kf SA:

ð7Þ

Kf is a constant of proportionality (in m yr− 1) that is assumed to vary with climate, mostly with precipitation rate (higher Kf values corresponding to a higher precipitation rate). Where sediment load, qs, obtained by integrating the volume of rock eroded from the upstream area, is smaller than the river capacity, erosion takes place, at a rate set by: ∂h q  qc = s wc lBR ∂t

ð8Þ

where lBR is a length scale characterizing the erodibility of bedrock (Kooi and Beaumont, 1994) (higher lBR values corresponding to more resistant rocks) and wc is the set river channel width. Where sediment load exceeds capacity deposition takes place at a rate set by: ∂h q  qc = s Ω ∂t

ð9Þ

where Ω is the surface area attached to each integration point and defined by the spatial discretization. Hillslope processes (landslides, soil creep, rainsplash, etc.) are represented by a simple diffusion equation that assumes that mass transport, q, is linearly proportional to slope, S (Kooi and Beaumont, 1994; Braun et al., 2001): q =  Kd S:

ð10Þ

Kd is a proportionality constant or diffusivity (in m2 yr− 1) that is a measure of the efficiency of the hillslope transport processes and is mostly proportional to precipitation rate (higher values of Kd corresponding to wetter climates). Here we present a series of numerical experiments in which shortening of an initially uniform crustal layer (thickness L and height h) is driven by an imposed basal velocity discontinuity (from v0 to 0, Fig. 1). This discontinuity is meant to represent shortening of the underlying substratum in a subduction-like process (as used by Malavielle (1984) or Willett et al. (1993) in two dimensions) or across any mechanical or lithological discontinuity (as used in Braun and Sambridge (1994) for example). Here we are interested in generalizing the work done on doubly-vergent critical wedges to the third dimension by assuming that the discontinuity has a finite width (in the third dimension). As shown in Fig. 1, this discontinuity is imposed at y =L/2 but along one half of the model only (from x= 0 to x =L/2) giving to the experiment a threedimensional character. The temperature is fixed at the base of the model (T =T0) and the initial temperature structure is obtained by assuming conductive transport only. In the subsequent paragraphs, we present results of our computations for a small range of model parameters and discuss primarily the mechanical implications of considering the third dimension.

Fig. 1. Experimental setup and boundary conditions. On all sides of the model, all velocity components are set to zero, except along the grey shaded side boundary where the normal component (parallel to the y-axis) is set to v0. Along the base of the model, the velocity is also set to zero, except along the grey shaded area where the velocity component parallel to the y-axis is set to v0. The thick black line represents the location of the S-line velocity discontinuity, i.e. where the velocity changes from v0 to 0.

The model general parameters and geometrical setup are shown in Fig. 1 and Table 1. Note that the thickness of the model, h, is relatively small, i.e. 10 km, because the base of the model being fixed, we cannot enforce the isostatic compensation. It is thus more appropriate to assume that the crustal layer being deformed is relatively thin such that the horizontal stress gradient created by its deformation and thickening can be sustained by the much thicker underlying lithosphere. In a first step, we will show results from two separate model runs or experiments, differing by the assumed efficiency of surface processes. The value of the parameters characterizing the surface processes is arbitrary and has been selected such that the system is capable of reaching an approximate flux steady-state in a time that corresponds to shortening by an amount of the order of once to twice the model thickness. In the first of the two experiments, both fluvial and hillslope processes are activated. In the second experiment, fluvial efficiency is twice as large as in the first model, while hillslope processes (and thus the system's ability to erode drainage divide) have been turned off. In the last part of this paper, we will show results for a third, hotter model, i.e. one in which the basal temperature has been increased from 250 to 400 °C and the erosion parameters have been reduced by a factor 2 in comparison with the first experiment. Strain softening is activated in all the experiments shown here. We have performed a few model experiments in which the strain softening was turned off; they show that the strain softening controls the width of the shear zones, but not their geometry (orientation). 3. Wedge development The evolution of the system is shown in Fig. 2 for the two different model setups described above, the low erosional efficiency experiment on the left, and the high efficiency experiment on the right. For each experiment five snap shots are shown at the same times in the model evolution, i.e. for the same amount of convergence at t = 1, 2, 3, 4 and 5 Myr. We show the evolution of the surface topography, much larger in the first of the two experiments as expected if the surface processes efficiency is diminished, as well as the distribution of strain rate along a vertical cross-section parallel to the y-axis and located next to the axis of symmetry of the model (the x =0 axis). In both cases, the model's response to basally-driven shortening is the formation of two oppositely-dipping shear zones accommodating reverse sense of shear. The resulting surface uplift is, at first, symmetrical, that is of equal width on

184

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

either side of the newly formed divide. As shortening proceeds, the system becomes asymmetrical with deformation accumulating along one of the two shear zones (the so-called retro-shear zone, (Willett et al., 1993)) while the other (the pro-shear zone) becomes a transient feature that is continuously advected through the V-shaped wedge that has developed between the two shear zones. This asymmetry is the result of the basal boundary condition, which imposes that the point where the velocity changes from v0 to zero (the so-called discontinuity or “Spoint”) is attached to one of the two converging plate, here the one on the retro-side. The surface topography simultaneously evolves into an asymmetrical shape with the divide being advected towards the retro-shear zone. The mean elevation is approximately 3 km and 1 km in the first and second experiments, respectively. This evolution is identical to that evidenced in many two-dimensional numerical (Willett et al., 1993; Batt and Braun, 1997) and analogic (Malavielle, 1984) experiments and this behavior is now well understood as that of a doubly-vergent, twodimensional critical wedge, where rocks can be considered to be at or close to failure, according to Coulomb's criterion. This evolution leads to a steady-state situation where the tectonic flux is compensated by the surface erosion flux as postulated and observed in many orogenic settings (Jamieson and Beaumont, 1988; Bernet et al., 2001; Willett et al., 2001; Willett and Brandon, 2002; Whipple and Meade, 2004). The early evolution of the model is faster, and thus steady-state is reached more rapidly, in the second model experiment where surface processes are more efficient, in agreement with Whipple and Meade (2006) that it is the efficiency of surface processes that determines the rate at which an orogenic wedge goes back to steady-state conditions following a perturbation, whether it is climatic or tectonic in origin. Note, however, that although a flux steady-state is reached in both experiments, a geomorphic steady-state is never reached with the geometry of the surface continuing to evolve with time. Indeed, while, at flux steadystate, the mean elevation, the width of the mountain range, the spacing between the river channels and the mean position of the divide are all relatively stable, the shape of the surface is very dynamic, mostly in response to advection of the landforms by the horizontal component of the velocity. This important and characteristic feature of natural systems (Willett et al., 2001) is commonly neglected in quantitative models of the coupling between tectonics and landscape evolution, and can only be properly modeled in three dimensions. The two experiments also differ in their structural evolution in that, in the case where surface processes are less efficient, the relatively large topographic gradient that has formed leads to a substantial perturbation of the stress state in the underlying crust. This perturbation leads to a local rotation of the direction of the principal stress axes along the x-axis which, in turn, causes a kink to develop in the retro-shear zone at mid-layer depth as well as a hinge structure to accommodate the bending of the wedge along that kink. To demonstrate this point, we show in Fig. 3 the orientation of the principal strain rate (and thus stress) axes as computed by the model in both experiments. In both experiments, the stress orientation is mostly determined by the horizontal compression resulting from the imposed basal shortening. As explained in the model description section, shear bands form at 45° to the direction of σ1 which is horizontal and aligned with the y-axis in most of the deforming wedge. σ3, the least compressive principal stress component is vertical and the intermediary stress, σ2 is horizontal, in the x-direction. In

185

Fig. 3, the retro- and pro-shear zones are labeled ‘1’ and ‘2’, respectively. In the first experiment (top panel), the topographic gradient that has developed in front of the retro-shear zone has led to a counter-clockwise rotation of the σ1 and σ3 axes around σ2 which has led to the formation of a kink in the retro-shear zone. A secondary structure has developed (labeled ‘1a’ in Fig. 3), originating at the kink in a direction parallel to the pro-shear zone and thus at 45° to the lower section of the retro-shear zone. In this shear band, the most compressive stress axis is vertical and the least compressive stress axis is horizontal. This structure accommodates extensional deformation. However, it is not sensu strictu a normal fault but a hinge structure because, although it corresponds to a discontinuity in the velocity field, velocity does not change in amplitude but only in direction. The formation of this structure is clearly related to the presence of the surface topographic gradient as it is not observed in the second, high erosional efficiency experiment. 4. Three-dimensional mechanical response The three-dimensional nature of the experiment we propose here comes from the existence of a free edge (at x = 0.5) between the deforming wedge and the adjacent ‘intact’ crust. To appreciate the complex nature and dynamics of this boundary, we show in Fig. 4 a horizontal cross-section, half-way through the thickness of the model, in which we have superimposed the principal stress orientations on a color map of the computed strain rate. The retro- and pro-shear zones (labeled ‘1’ and ‘2’ in Fig. 4) appear as two zones of intense deformation separated by the undeformed, yet uplifting plug. A third structure (labeled ‘4’ in Fig. 4) accommodates the differential motion between the advancing crustal layer subjected to the imposed basal velocity and the neighboring, static crust. Along all of its length this structure is a pure strike-slip shear zone characterized by horizontal σ1 and σ3 axes both at 45° to the strike of the shear zone. The most compressive axis lies in the plane x = y. The retro-shear zone (labeled ‘1’ in Fig. 4) does not connect in a simple manner with the strike-slip shear zone (labeled ‘4’ in Fig. 4): it forms a curved arc along the outside corner of the model, above the termination point of the imposed basal discontinuity. Along this external section of the retroshear zone (labeled ‘1b’ in Fig. 4), the maximum compressive stress orientation rotates by almost 90° from parallel to the y-axis to subparallel to the x-axis, and remains perpendicular to the strike of the curved shear zone, indicating that all along its length, this shear zone ‘1’ + ‘1b’ is a pure thrust, i.e. it accommodates a pure reverse motion. The strike-slip deformation is entirely accommodated by the through-going strike-slip shear zone that runs parallel to the y-axis (labeled ‘4’ in Fig. 4). This partitioning of the deformation between a pure thrust and a pure strike-slip shear zone is most visible in the first of the two experiments (Fig. 4a). In contrast, the pro-shear zone (labeled ‘2’ in Fig. 4) connects in a simpler manner to the strike-slip fault (labeled ‘4’ in Fig. 4): away from the x = 0 axis of symmetry and in the vicinity of the strike-slip shear zone, it also becomes curved, and more so in the first of the two experiments, but the direction of the most compressive stress remains parallel to the y-axis, implying that it is an oblique thrust, because the direction of the most compressive stress is not orthogonal to the strike of the shear zone. The deformation is not partitioned; there is only one structure (labeled ‘2b’ in Fig. 4) accommodating the combined motions.

Fig. 2. Evolution of two model experiments differing by the efficiency of the surface processes (panels a to e correspond to low erosional efficiency and f to j to high efficiency). Surface color is proportional to the predicted height (or elevation); color range is dynamically adjusted to the range of elevation in each panel; colors on the vertical cross-section at x = 0 represent the second invariant of the computed strain rate tensor; for ease of comparison, the color range is the same for all panels. The black arrows represent the computed velocity field. The two planes of high strain rate intensity correspond to major structures (the retro- and pro-shear zones) accommodating the shortening imposed by the basal velocity discontinuity. The two model experiments differ by the height of the predicted topography but also by the geometry of the retro-shear zone. Note that the axes given in panel a indicate the direction not the absolute position of the system of reference used to describe the solution. The high strain rate values shown along the top surface (i.e. at the interface between rock and atmosphere) in this figure and in Figs. 3, 5 and 7 are artifacts of the plotting package we have used.

186

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

Fig. 3. Vertical cross-section in the [y − z] plane at x = 0 showing (in subdued colors) the computed strain rate and the orientation of the principal strain rate (and thus stress) axes for the two model experiments (a low and b high erosional efficiency). Labeled structures are: 1 = retro-shear zone; 2 = pro-shear zone; 1a = upper section of the retro-shear zone forming a kink with the lower section; and 3 = hinge. The double white arrow and double black arrow indicate the orientation of the most compressive principal stress (σ1) and least compressive principal stress (σ3), respectively.

Where the two shear zones (‘2b’ and ‘4’) connect, the principal stress triad undergoes a complex rotation along both the vertical and y-axis such that both σ1 and σ3 lie in the horizontal plane. Looking at the deformation in a plane perpendicular to the y-axis, half-way through the length of the model, i.e. along the basal velocity discontinuity (Fig. 5), one notices that thrusting along the curved thrust is accompanied by deformation along a secondary, oppositelydipping shear zone (labeled ‘5’ in Fig. 5) that is best seen in the first of the two experiments where surface processes are less efficient and substantial surface topography has developed (Fig. 5a). This structure is made of two sections: one, in the lower half of the model, that is nearly vertical and originates from the end of the basal velocity discontinuity, the other, in the upper half of the model, that is dipping at approximately 45° in a direction opposite to that of the main thrust plane (labeled ‘1b’ in Fig. 5). Deformation on this structure is dominated by oblique normal movement in the upper part and strike-slip deformation in the lower part. To better understand how

this structure connects to the others, we have plotted the strain rate along two cross-sections (Fig. 6), a vertical one passing through the velocity discontinuity and a horizontal one. Contrary to all other figures, we only show here results for the first of the two experiments. In the first of the two panels (Fig. 6a) the horizontal cross-section is placed half-way through the model, in the second panel (Fig. 6b), the horizontal section is closer to the surface. In the lower half of the model (Fig. 6a), the secondary structure (labeled ‘5’) connects to the vertical, purely strike-slip shear zone (labeled ‘4’ in Fig. 6a), parallel to the y-axis; in the upper half (Fig. 6b), the secondary structure (labeled ‘5’) connects to the hinge structure (labeled ‘3’ in Fig. 6b) that is antithetic to the main thrust (or retro-shear zone). Variations in the velocity field (not shown here for clarity) across this structure demonstrate that, in its lower section, it is a pure strike-slip shear zone while, in its upper section, it is a hinge that accommodates the change in velocity direction (not amplitude) from parallel to the y-axis to parallel to the x-axis. In other words it accommodates the

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

187

Fig. 4. Horizontal cross-section through the model half-way through the model thickness for the two model experiments (a low and b high erosional efficiency). Principal stress axes are superimposed on a color map of the second invariant of the computed strain rate. Labeled structures are: 1 = retro-shear zone; 2 = pro-shear zone; 1b = curved section of the retro-shear zone; 2b = curved section of the pro-shear zone; and 4 = strike-slip shear zone. The curved section of the retro-shear zone is a pure thrust while the curved section of the pro-shear zone is an oblique thrust.

change in velocity direction resulting from the ‘lateral escape’ motion accommodated by the curved section of the retro-shear zone. 5. Effect of basal viscous creep Increasing the temperature at the base of the crustal layer from 250 to 400 °C leads to a very different mechanical and topographical response, with the formation of a thin decollement at the base of the crust, widening of the deforming region and uplift of a broad topographic plateau (Fig. 7). One may regard this imposed basal temperature as ‘extreme’ and difficult to justify at the base of a 10 km thick crustal layer; however, our purpose here is to demonstrate the importance of the mechanical behavior of the wedge material, from purely brittle to thermally-activated creep, on its three-dimensional response to basally-driven shortening. The model dimensions are thus kept identical for easy comparison with the previous two model runs. Surface erosion efficiency is twice as low as in experiment 1 (see Table 1 for the exact value of the erosional coefficient).

The 2D evolution is different to the one described above, but similar to that evidenced in other studies (Ellis et al., 1998; Beaumont et al., 2000) where a temperature increase has led to thermallyactivated dislocation creep becoming the dominant deformation mechanism. In our simulation, the material is also subject to strain softening in the ductile regime. The system first evolves in a similar way to that of the pure brittle case, with the formation of oppositelydipping crustal-scale shear zones (labeled ‘1’ and ‘2’ in Fig. 7a) that accommodate the basally-driven shortening and uplift of a triangular plug. However, with accumulating deformation and the growth of steep surface topography, a horizontal decollement (labeled ‘3’ in Fig. 7) develops along the base of the model, first in the pro-direction (panels a, b and c in Fig. 7), then in the retro-direction (panels c and in Fig. 7). The development and evolution of this horizontal basal decollement is facilitated by the assumed strain softening component of the viscous rheology. When this decollement reaches the far-side (y = 0) boundary of the model, it cannot propagate any further, a thrust forms (labeled ‘5’ in Fig. 7) and uplift takes place along the

188

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

Fig. 5. Vertical cross-section through the model in a plane passing through the imposed basal discontinuity (i.e. at y = 0.5) for the two model experiments (a low and b high erosional efficiency). Principal stress axes are superimposed on a color map of the second invariant of the computed strain rate. Labeled structures are: 1b = curved section of the retro-shear zone; and 5 = oblique hinge. The oblique hinge is made up of two sections; the upper section dips at approximately 45° and the lower section is nearly vertical.

y = 0 far boundary. The location of this structure is thus imposed by the finite extent of the numerical experiment, and the interaction of the internal, mechanical response of the system with an artificially imposed boundary condition. In a natural system, such decollements propagate along distances that are dictated by the viscosity of the decollement, itself determined by the temperature and rheology at the base of the crustal layer, and the depth of the decollement (Braun and Beaumont, 1989; Vanbrabant et al., 2002). In many instances, horizontal heterogeneity in the crust will also dictate the distance over which such decollements propagate (Beaumont et al., 2000; Braun and Pauselli, 2004), with strong or cold inclusions inhibiting propagation and facilitating the nucleation of thrust faults. Along the retro-side of the wedge, a secondary retro-shear zone develops (labeled ‘4’ in Fig. 7) outboard of the initial wedge, that is connected to the imposed discontinuity (the ‘S’ point in Fig. 7) by the propagation of the basal decollement in the retro-direction. This propagation of the deformation in the retro-direction is likely to be

the result of the inability of the system to propagate further in the other direction due to the finite length of the model. At this stage of its evolution, the wedge mechanics becomes disconnected from the linediscontinuity; its dynamics resembles more that of a side-driven or ‘vise-type’ compressional orogen (Ellis et al., 1998). Movement along the most external structures (4 and 5 in Fig. 7) leads to the formation of a broad, plateau-like surface that is composed of actively uplifting rugged mountain ranges along its boundaries and a very flat interior, that forms by sediment transport and deposition. Sediment accumulation lasts until a channel draining along the side of the plateau-like topographic feature is capable of breaking through the local drainage divide (between the plateau and the low-lying, underformed adjacent areas) and rapidly drains the thick sediment pile. This rapid unloading is concomitant with the formation of yet another structure (labeled ‘6’ in Fig. 7), illustrating the complexity of interactions that one develops in a fully-coupled three-dimensional model of crustal deformation and surface processes.

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

189

Fig. 6. Double cross-section through the model for the first of the two model experiments, i.e. the low erosional efficiency, high topography one. Vertical section is along a plane passing through the basal discontinuity (i.e. at y = 0.5); the horizontal plane passes half-way through the model thickness in panel a and is closer to the surface (at 0.116 × h, the model initial thickness, from the base of the model) in panel b. Labeled structures are: 1 = retro-shear zone; 2 = pro-shear zone; 1b = curved section of the retro-shear zone; 3 = hinge; 4 = strike-slip shear zone; and 5 = oblique hinge. Along its lower section the oblique hinge (labeled 5) connects to the strike-slip shear zone (labeled 4), while in its upper section, it connects to the hinge (labeled 3).

We complete this already dynamic evolution of the system by examining the three-dimensional nature of the solution (Fig. 8) and considering the evolution of the structures in a horizontal plane located near the surface (i.e. at z = 0.9) with a special emphasis given to the external boundary of the deforming wedge (i.e. at x = 0.5). In the initial stages of deformation (panel a in Fig. 8), the threedimensional strain rate and its partitioning between pure thrust, pure strike-slip and oblique structures are similar to those of the purely brittle/frictional case shown above (Fig. 4). However, as deformation progresses, the pro-shear zone is abandoned and another structure (or secondary pro-shear zone) develops along the bottom part of the model (i.e. near the y = 0 boundary in this top view of the model). This transfer is progressive (panels b to e in Fig. 8) and evolves from right to left, i.e. from the external boundary of the deforming wedge towards its interior. The reason for this progressive transfer is the preferential viscous strain weakening of the wedge material along its side boundary, in response to the accumulated strike-slip deformation.

Note that this structure is, at all times, strongly controlled by, and thus pinned to, the edge of the model which is an artifact of the finite size of the experiment. With further deformation, another or secondary retro-thrust develops outboard of the retro-shear zone which progressively accommodates thrusting. At the end of the model run (panel h), the thrusting is equally distributed among the two structures. Note that the outmost curved section of this secondary structure is, at all times in its evolution, a pure thrust, i.e. the most compressive stress component is everywhere perpendicular to the strike of the shear zone. On the contrary, the connection between the first retro-shear and the strike-slip shear zones on the right-hand side of the wedge (at x = 0.5) is at first a curved pure thrust but evolves into an oblique structure. The strike-slip shear zone evolves through time too as it curves half-way through its length to connect the edge of the oblique retro-thrust with the imposed discontinuity in horizontal velocity along the bottom boundary (at x = 0.5 and y = 0).

190

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

Fig. 7. Evolution of a model experiment characterized by a higher basal temperature and in which thermally-activated creep is thus activated. Surface color is proportional to the predicted height (or elevation); colors on the vertical cross-section at x = 0 represent the second invariant of the computed strain rate tensor. White arrows and labels indicate structures discussed in the text. The ‘S’ label indicates the position of the imposed, fixed velocity discontinuity.

Finally, it is worth noting the formation of the final structure, antithetic to the secondary pro-thrust that developed earlier at the bottom (y = 0) of the model (panels f to h). This structure evolves also from right to left, but it is not a ‘fixed’ structure like the bottom

thrust; it is also progressively advected by the horizontal tectonic movement accommodated by thrusting along the two retro-shear zones near the top of the model. This explains its curved shape (panel h).

Fig. 8. Evolution of the model experiment shown in Fig. 7 but along a horizontal plane at z = 0.9. Colors are proportional to the second invariant of the stain rate tensor and small black arrows show the direction and amplitude of the principal axis of the strain rate (and thus stress) tensor. Large white and black arrows indicate σ1 and σ3 directions, respectively.

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

191

192

J. Braun, P. Yamato / Tectonophysics 484 (2010) 181–192

6. Conclusions We have shown that three-dimensional effects arising from the finite width of an orogenic wedge are complex and certainly not intuitive. Most noticeable are two structures: (a) the curved section of the retro-shear zone around the edge of the imposed velocity discontinuity and (b), antithetic to it, the complex hinge structure that accommodates the change in velocity direction resulting from a kink in the retro-shear zone and the lateral escape of crustal material, which, in turn, are both the product of the finite amplitude topography that has developed in the model experiment characterized by low surface processes efficiency. The geometry and nature (sense of movement) of these structures could not be predicted from a simple stress balance argument as may be the case in a simple two-dimensional critical wedge. Our results therefore demonstrate the power of the numerical model we have used which provides unprecedented insight into the dynamics of orogenic systems. To be able to predict the formation and evolution of complex three-dimensional structures, to display their three-dimensional geometry, their interactions with each other and the free evolving surface and, finally, to compute and display the orientation of the stress or strain rate principal axes with accuracy is an exciting novel development. The theoretical insight that is gained complements what could be observed in the field or in an analogue experiment. Our numerical model, and the results of the computations made with it, should therefore be regarded as tools to help the structural geologist interpret observations in similar geological settings but in a broad range of environments and geometries/scales. Although some scaling issues regarding the balance between topographically generated gravitational stress gradients and the intrinsic strength of crustal rocks should be considered before applying our results to any given geological setting, one could consider using them to interpret the formation and evolution of complex three-dimensional structures at the scale of a small orogen down to that of a single sedimentary decollement. We also note that major differences among model results originate in the assumed efficiency of the surface processes. Where surface processes are not efficient and relatively high topography develops at the surface of the model, the main thrust structures are kinked and antithetic hinge structures form. On the contrary, thrust faults are planar and no hinge develops when surface processes are efficient and the topography remains modest. Lowering the intrinsic strength of the crust by increasing the basal temperature and activating dislocation creep, makes the system more responsive to the loading/unloading effect of surface topography. In this situation, we have shown that the system rapidly evolves towards a widening of the region of active deformation with the formation of offset structures on either side of the original thrusts (the retro- and pro-shear zones). This, in turn, causes uplift over a wider area and the formation of a relatively flat inward-draining plateau surrounded by actively eroding mountain belts. Our computations do clearly show the strong couplings that exists between surface processes and crustal deformation, at least at the scale considered here, i.e. 10–100 km. Acknowledgments The work described in this paper has been supported by a Postdoctoral Fellowship (P. Y.) and a Chaire d'Excellence (J. B.) awarded by the Agence Nationale de la Recherche (ANR). We thank T. Yamasaki and S. Cox for very useful reviews of this manuscript. References Batt, G.E., Braun, J., 1997. On the thermomechanical evolution of compressional orogens. Geophys. J. Int. 128, 364–382.

Beaumont, C., Fullsack, P., Hamilton, J., 1992. Erosional control of active compressional orogens. In: McClay, K.R. (Ed.), Thrust Tectonics. Chapman, New York, pp. 1–18. Beaumont, C., Munoz, J.A., Hamilton, J., Fullsack, P., 2000. Factors controlling the Alpine evolution of the central Pyrenees inferred from a comparison of observations and geodynamical models. J. Geophys. Res. 105, 8121–8145. Bernet, M., Zattin, M., Garver, J.I., Brandon, M.T., Vance, J.A., 2001. Steady-state exhumation of the European Alps. Geology 29, 35–38. Braun, J., 1993. Three-dimensional numerical modelling of compressional orogens: thrust geometry and oblique convergence. Geology 21, 153–156. Braun, J., Beaumont, C., 1989. Dynamical models of the role of crustal shear zones in asymmetric continental extension. Earth Planet. Sci. Lett. 93, 405–423. Braun, J., Beaumont, C., 1995. Three-dimensional numerical experiments of strain partitioning at oblique plate boundaries: implications for contrasting tectonic styles in the southern Coast Ranges, California and central South Island, New Zealand. J. Geophys. Res. 100, 18,059–18,074. Braun, J., Sambridge, M., 1997. Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretization. Basin Res. 9, 27–52. Braun, J., Pauselli, C., 2004. Tectonic evolution of the Lachlan Fold Belt (Southeastern Australia): constraints from coupled numerical models of crustal deformation and surface erosion driven by subduction of the underlying mantle. Phys. Earth Planet. Inter. 141, 281–301. Braun, J., Sambridge, M., 1994. Dynamical Lagrangian Remeshing (DLR): a new algorithm for solving large strain deformation problems and its application to faultpropagation folding. Earth Planet. Sci. Lett. 124, 211–220. Braun, J., Heimsath, A.J., Chappell, J., 2001. Sediment transport mechanisms on soilmantled hillslopes. Geology 29, 683–686. Braun, J., Thieulot, C., Fullsack, P., Dekool, M., Beaumont, C., Huismans, R., 2008. DOUAR: a new three-dimensional creeping flow numerical model for the solution of geological problems. Phys. Earth Planet. Inter. 171, 76–91. Dahlen, F.A., 1984. Noncohesive critical coulomb wedges: an exact solution. J. Geophys. Res. 89, 10,125–10,133. Ellis, S., Beaumont, C., Jamieson, R.A., Quinlan, G., 1998. Continental collision including a weak zone: the vise model and its application to the Newfoundland Appalachians. Can. J. Earth Sci. 35, 1323–1346. Gerbault, M., Davey, F., Henrys, S., 2002. Three dimensional lateral crustal thickening in continental oblique collision: an example from the Southern Alps, New Zealand. Geophys. J. Int. 150, 770–779. Jamieson, R.A., Beaumont, C., 1988. Orogeny and metamorphism—a model for deformation and pressure–temperature–time paths with application to the Central and Southern Appalachians. Tectonics 7, 417–445. Konstantinovskaia, E., Malavieille, J., 2005. Erosion and exhumation in accretionary orogens: experimental and geological approaches. Geochem. Geophys. Geosyst. 60 (2). doi:10.1029/2004GC000794. Kooi, H., Beaumont, C., 1994. Escarpment evolution on high-elevation rifted margins: insights derived from a surface processes model that combines diffusion, advection and reaction. J. Geophys. Res. 99, 12,191–12,209. Kooi, H., Beaumont, C., 1996. Large-scale geomorphology: classical concepts reconciled and integrated with contemporary ideas via a surface processes model. J. Geophys. Res. 101, 3361–3386. Koons, P.O., 1990. Two-sided orogen: collision and erosion from the sandbox to the Southern Alps, New Zealand. Geology 18, 679–682. Koons, P.O., 1994. Three-dimensional critical wedges: tectonics and topography in oblique collisional orogens. J. Geophys. Res. 99, 12,301–12,315. Malavielle, J., 1984. Modélisation expérimentale des chevauchements imbriqués: application aux chaînes de montagnes. Bull. Soc. Geol. Fr. 26, 129–138. Tullis, J., 2002. Deformation of granitic rocks: experimental studies and natural examples. In: Karato, S.-I., Wenk, H.-R. (Eds.), Plastic Deformation of Minerals and Rocks: Review in Mineralogy and Geochemistry, vol. 51, pp. 51–95. Vanbrabant, Y., Braun, J., Yongmans, D., 2002. Models of passive margin inversion: implications for the Rhenohercynian fold-and-thrust belt, Belgium and Germany. Earth Planet. Sci. Lett. 202, 15–29. van der Beek, P.A., Braun, J., 1998. Numerical modelling of landscape evolution on geological time scales: a parameter analysis and comparison with the southeastern highlands of Australia. Basin Res. 10, 49–68. Whipple, K.X., Meade, B.J., 2004. Controls on the strength of coupling among climate, erosion, and deformation in two-sided, frictional orogenic wedges at steady state. J. Geophys. Res. 109, F01011. doi:10.1029/2003JF000019. Whipple, K.X., Meade, B.J., 2006. Orogen response to changes in climatic and tectonic forcing. Earth Planet. Sci. Lett. 243, 218–228. Willett, S.D., 1999a. Orogeny and orography: the effects of erosion on the structure of mountain belts. J. Geophys. Res. 104, 28,957–28,981. Willett, S.D., 1999b. Rheological dependence of extension in wedge models of convergent orogens. Tectonophysics 305, 419–435. Willett, S.D., Brandon, M.T., 2002. On steady states in mountain belts. Geology 30, 175–178. Willett, S.D., Beaumont, C., Fullsack, P., 1993. Mechanical model for the tectonics of doubly-vergent compressional orogens. Geology 21, 371–374. Willett, S.D., Slingerland, R., Hovius, N., 2001. Uplift, shortening, and steady state topography in active mountain belts. Am. J. Sci. 301, 455–485.