Comparing thin-sheet models with 3-D multilayer models for continental collision

Summary

Various models have been proposed to explain tectonic deformations during continent collision. A frequently applied model is the thin viscous sheet model which is however not fully 3-D and assumes a priori diffuse thickening as the dominant deformation style.

We compare a fully 3-D multilayer numerical model with a corresponding thin viscous sheet numerical model for the scenario of continent indentation. In our comparison we focus on the three basic viscous deformation styles thickening, buckling (folding) and lateral crustal flow. Both numerical models are based on the finite element method (FEM) and employ either a linear or power-law viscous rheology. The 3-D model consists of four layers representing a simplified continental lithosphere: strong upper crust, weak lower crust, strong upper mantle and weak lower mantle. The effective viscosity depth-profile in the 3-D model is used to calculate the depth-averaged effective viscosity used in the thin-sheet model allowing a direct comparison of both models.

We quantify the differences in the strain rate and velocity fields, and investigate the evolution of crustal thickening, buckling and crustal flow resulting from the two models for two different phases of deformation: (1) indentation with a constant velocity and (2) gravitational collapse after a decrease of the indenting velocity by a factor of 5. The results indicate that thin-sheet models approximate well the overall large-scale lithospheric deformation, especially during indentation and for a linear viscous rheology. However, in the 3-D model, additional processes such as multilayer buckling and lower crustal flow emerge, which are ignored in the thin-sheet model but dominate the deformation style in the 3-D model within a range of a few hundreds of kilometres around the collision zone and indenter corner. Differences between the 3-D and thin-sheet model are considerably larger for a power-law viscous than for a linear viscous rheology. Buckling and lower crustal flow are significant in the 3-D model with power-law viscous rheology. For example, fibre strain rates due to buckling can be several hundred per cent different to the depth-averaged strain rate and the lateral mass flow of lower crustal material can be up to six times more than the flow of upper crustal material. Our results also show that the horizontal velocity fields of the upper crust and upper mantle remain nearly identical in the 3-D model during indentation despite their mechanical decoupling due to an intermediate weak lower crust. This result questions the validity of using similarities between velocities from the surface global positioning system (GPS) and mantle shear wave splitting data as evidence for a mechanically coupled lithosphere.

3-D multilayer models provide a more complete picture of continental collision than thin-sheet models as they enable studying the timing, locality and relative importance of different processes simultaneously which is especially important for the hundreds of kilometre scale around the collision zone and indenter corners. 3-D models are, however, still computationally challenging and we, therefore, also present results of a computational performance test of several solution algorithms.

Numerical solutions, Numerical approximations and analysis, Continental tectonics: compressional, Dynamics of lithosphere and mantle

1 Introduction

Continental tectonic deformations during, for example, the India–Asia collision are frequently explained with thin viscous sheet models (e.g. Bird & Piper 1980; England & McKenzie 1982; England & Houseman 1986; Houseman & England 1986; Jimenez-Munt & Platt 2006; Vergnolle et al. 2007). Such thin-sheet models were first applied to the India–Asia collision by England & McKenzie (1982; 1983) and to crustal deformations in southern California by Bird & Piper (1980). Thin-sheet models have been applied to explain certain large-scale features of the India–Asia collision. For example, Houseman & England (1986) obtained good agreement between the thickening of the crust computed from thin-sheet models and the known crustal thickness of the Tibetan Plateau. Jiménez-Munt & Platt (2006) modelled the flow of crustal material around the indenting corner of the Indian Plate in the area of the eastern Himalayan syntaxis and were able to produce a model with horizontal surface velocities which were in good agreement with the observed GPS velocity data. Flesch et al. (2001) used the thin-sheet approximation in combination with geological and geophysical data to estimate the vertically averaged effective viscosity distribution and deviatoric stress field of the India–Asia collision.

In addition to explaining several large-scale features of continental collision, thin-sheet models have the advantage that they yield a mathematically 2-D problem, which describes a 3-D deformation field. The advantages of reducing the dimensionality of the problem are, first, that analytical solutions are easier to obtain and, second, that numerical approximations are computationally inexpensive to obtain. Reducing the 3-D problem to a 2-D problem is achieved by assuming that the lithosphere can be represented as a continuous viscous sheet with effective, depth-averaged parameters. Vertical gradients of horizontal velocities are assumed to be zero (i.e. strain rate is constant with depth), flexural deformations (i.e. bending and buckling) are ignored and as a result the lithosphere deforms only by diffuse thickening (or thinning). A thin-sheet model seems intuitive for modelling lithospheric deformation during large-scale continental collision, because the thickness of the lithosphere (∼100 km) is considerably smaller than the horizontal dimensions (∼1000′s of km) of the lithospheric area affected by the collision. However, the particular sheet-like or plate-like geometry of the lithosphere, and especially its internal mechanical stratification (e.g. upper and lower crust and upper and lower lithospheric mantle), can alternatively be used as an argument for the application of so-called thin-plate models to study lithospheric deformation, where, in contrast to thin-sheet models, buckling and bending are first-order processes during lithospheric shortening (e.g. Burg & Podladchikov 1999; Cloetingh et al. 1999; Schmalholz et al. 2002).

The restrictions in the deformation style of thin-sheet models led Medvedev & Podladchikov (1999a,b) to enhance the original thin-sheet model and to apply their ‘extended thin-sheet’ model to a variety of geodynamic problems such as, for example, lithospheric folding. From a theoretical point of view, folding is expected to take place during shortening of a mechanically stratified lithosphere because it requires less strain energy than diffuse thickening, and for some parameters is the preferred deformation mode of mechanically stratified bodies. Such is likely the case with the continental lithosphere, which is shortened subparallel to its stratification during continental collision (Biot 1961; Burg & Podladchikov 1999; Burg & Schmalholz 2008; Cloetingh et al. 1999; Kaus & Schmalholz 2006). The main force inhibiting lithospheric folding is gravity and whether folding or thickening dominates during lithospheric shortening is controlled by the Argand number (i.e. the ratio of buoyancy stress to viscous compressive stress; England & McKenzie 1982), the effective viscosity ratio between lithospheric sublayers (e.g. upper and lower crust and upper mantle) and the power-law stress exponent of each sublayer (Schmalholz et al. 2002).

Furthermore, other deformation processes, such as lower crustal flow, thrusting along large-scale shear zones and strike slip faulting are also neglected in the standard thin-sheet models. Consequently, a variety of other models has been applied to study continental collision zones, in particular the India–Asia collision. For example, lower crustal flow relative to the upper crust and upper mantle is considered in channel flow models (e.g. Beaumont et al. 2001) and lower crustal flow models (e.g. Royden et al. 2008). Thrusting along plastic shear zones is modelled with crustal wedge models (e.g. Koons et al. 2002), and strike slip faulting has been studied with theoretical slip line models and analogue experiments (e.g. Tapponnier & Molnar 1976).

In general, a drawback of many models applied to continental collision is that the chosen model setup, rheology and boundary conditions typically impose a certain deformation style a priori. For example, thickening is promoted in thin viscous sheet models, whereas significant lower crustal flow is promoted in some channel flow models (e.g. Clark & Royden 2000; Medvedev & Beaumont 2006). Therefore, it is difficult to assess if the deformation style observed in a model is a good approximation of the natural deformation style, or simply the inevitable result of the a priori assumptions made in the model. It is thus desirable to use models that are able to simultaneously simulate several deformation styles and to determine under what conditions each deformation style is dominant.

In this study, we use a 3-D multilayer numerical model based on the finite element method to study different ductile deformation styles during continental indentation. The 3-D model consists of four layers representing a simplified mechanically stratified lithosphere: strong upper crust, weak lower crust, strong upper mantle and weak lower mantle. Three major viscous deformation styles can be distinguished in the model: (i) thickening and thinning, (ii) buckling (i.e. folding) and (iii) lower crustal flow. Our 3-D model utilizes a similar setup to that employed by the thin-sheet models of England & Houseman (1986) and Houseman & England (1986) because these thin-sheet models have been successful in explaining several large-scale features of the India–Asia collision. This also implies that our study is restricted solely to viscous (linear and power-law) rheologies, as used in the studies of England & Houseman (1986) and Houseman & England (1986). As such we can directly compare the results of thin-sheet models with the results of our 3-D multilayer models.

The objectives of this study are (i) to investigate under what conditions thin-sheet models are acceptable approximations to mechanically stratified 3-D models for simulations of continental indentation and (ii) to study the relative importance of the viscous deformation processes thickening, folding and crustal flow during continental indentation. Recently, Garthwaite & Houseman (2011) presented a comparison of thin-sheet models with 3-D models consisting of a single homogeneous viscous layer (i.e. no buckling and no lower crustal flow). They found, amongst other things, that the greatest differences between thin-sheet and 3-D models occur around the indenter corner and that differences increase with increasing power-law exponent.

3-D multilayer numerical models are still computationally challenging and require a careful choice of the algorithm used to solve the discrete system of equations. Therefore, we compared the performance of several solution algorithms to determine the most suitable algorithm for conducting the presented 3-D multilayer simulations.

2 Methods

2.1 Governing equations

The thin-sheet and 3-D models are based on the continuum mechanics conservation equations. The conservation equation representing the momentum balance is

1

with σ_ij being the total stress tensor, ρ the density, g the gravity and . Indices i,j run from 1 to 3 and refer to the three spatial coordinates x, y, z, respectively. According to the summation convention, repeated indices are summed. The superscript T indicates the transpose of a vector or matrix.

The incompressibility condition, with velocity vector and , is

2

with the symbol · representing a scalar product.

The deviatoric stress tensor is

3

where is the pressure and δ_ij the Kronecker delta.

We apply the constitutive relation

4

where is the strain rate tensor and is the second invariant of the strain rate tensor. H is the rheological coefficient and n the power-law exponent (H= 2μ when n= 1, in which case μ is the linear Newtonian viscosity).

Substituting eqs (3) and (4) into eq. (1) leads to

5

where

6

Eqs (2) and (5) represent the four governing equations in 3-D for the four unknowns v_x, v_y, v_z and p.

The thin-sheet equations are approximations of these 3-D equations, in which all parameters are averaged with depth and the vertical shear stress components are assumed to be negligible, that is, τ_xz=τ_yz= 0. Considering this assumption, inserting eq. (3) into eq. (1), leads to the following expressions

6

8

Here, the indices i, j now only run from 1 to 2 and represent the two horizontal coordinates x and y, respectively. Eq. (7) can be solved for pressure as

9

ρ_c and ρ_m are the densities of the crust and the mantle, respectively. The base of the lithosphere is located at z= 0 and in isostatic equilibrium, the top of the lithosphere is defined to be at z = L+h (England & McKenzie 1983). When the lithosphere is in isostatic equilibrium, h corresponds to h=S−b with (Fig. 1), where S is the crustal thickness.

Figure 1

Parameters of the lithosphere used to derive the thin-sheet equations. S is the crustal thickness and L the mantle thickness. h is the elevation and b the isostatic compensation depth. ρ_c and ρ_m are the densities of the crust and the mantle, respectively. The base of the lithosphere is at z= 0.

From eqs (2) and (4)τ_zz is defined as

10

Next, each parameter, which we denote by ϕ, is vertically averaged over the thickness of the lithosphere, by applying

11

Applying eq. (10) to eq. (6), and making use of eq. (8), gives

12

with

13

Eq. (12) can be rewritten using the lithospheric pressure p* which yields

14

p*_cand p*_m are the pressure in the crust and in the mantle, respectively. Evaluating the integrals and using eq. (9), leads to the following expression (the detailed evaluation is given in the appendix)

15

with , the isostatic amplification factor.

Substituting eqs (4) and (14) into eq. (11) leads to (where the bars denoting averages have been dropped for clarity)

16

with The first term on the right-hand-side of both equations in eq. (15) is a constant and thus vanishes when differentiated. By non-dimensionalizing using H, L, the right-hand-side of eq. (15) can be simplified to

17

where

18

and is the reference strain rate. Ar is the Argand number which is a measure of the relative contributions of the forces arising from variations of crustal thickness and the forces required to deform the fluid at a reference strain rate (England & McKenzie 1982; England et al. 1985). Other non-dimensional variables (denoted via stars) are defined as follows

19

After non-dimensionalization, eq. (15) can be re-written as (stars of non-dimensional parameters are dropped for clarity)

20

with

21

Using eq. (2) we can account for the change in crustal thickness S

20

for the thin-sheet model, eqs (19) and (20) represent a system of three equations, which can be solved for S, v_x and v_y, by applying the initial and boundary conditions. Eq. (20) can be solved after eq. (19) which is a 2-D mathematical problem providing the two horizontal velocities that can then be used to calculate (update) the crustal thickness.

2.2 Discretization

2.1 3-D multilayer model

The governing equations for the conservation of momentum, eq. (2) and the conservation of volume, eq. (5) are discretized using a standard mixed finite element method (FEM) (Brezzi & Fortin 1991; Hughes 1987). As the model domain possesses a layered geometry, and due to the geometric simplicity of this domain, we utilize hexahedral elements to discretize our domain. In the mixed finite element approach, the discrete velocity and pressure solutions are represented by different basis functions. The choice of these basis functions (or element types) requires careful consideration. In this work, for the simulations in Section 3, we utilize the element type Q₂P₋₁ which consists of a triquadratic basis function for velocity and a piecewise discontinuous linear basis function for pressure (Brezzi & Fortin 1991). This element is suitable for solving incompressible problems as it is inf-sup stable. Following the discretization of the Stokes equations, we obtain the following discrete system

23

with the discrete gradient of the stress tensor, the discrete gradient, Q^T_3-D the discrete divergence and the discrete body force. Matrices and contain spatial derivatives of the velocity shape functions. Vectors N_p, N_v contain the pressure and the velocity shape functions.

The mixed finite element method described above has been implemented within the code LaMEM [Lithosphere and Mantle Evolution Model (e.g. Schmeling et al. 2008)]. LaMEM is a parallel (distributed memory) FEM code developed using the PETSc library (Balay et al. 2009), which provides a flexible mechanism to change and configure the solution strategy used to compute the discrete velocity and pressure from eq. (21). The incompressible flow problem is solved using either a fully coupled (v, p) solver, a decoupled Schur complement approach, or a specialized, ‘hand-rolled’ conjugate gradient algorithm which iterates on the Schur complement, but performs simultaneous updates of the v and p solution at each iteration (Verfürth 1984). The Schur complement methods require the solution of a velocity subproblem. Picard iterations (i.e. successive solution substitution) are used to solve the non-linear problem when the power-law rheology is introduced. LaMEM has been benchmarked for 3-D power-law viscous folding, where the numerical solution was compared with the analytical solution of Fletcher (1991) for the growth rate of a 3-D single-layer fold. Furthermore, it has been benchmarked for the linear stages of the Rayleigh Taylor instability (e.g. Kaus & Podladchikov 2001) and for a free surface subduction benchmark (Schmeling et al. 2008).

Following the solution of the velocity and pressure, the nodes of the finite element mesh are advected in a fully Lagrangian manner with a method which is first order accurate in space and time. The material properties and rheology are carried with the mesh, and the new mesh geometry thus defines the updated material configuration.

2.2.2 Thin-sheet model

The discretization of the 2-D, thin-sheet model is performed using a standard finite element method (Hughes 1987). Eq. (19) is discretized in space using biquadratic (nine node) elements resulting in a system of equations of the form:

24

with and Matrices and contain spatial derivatives of the velocity shape functions. Crustal thickness S in the right hand side of eq. (22) is updated on every integration point in time via explicit, first-order finite differences. From eq. (20) we can thus write,

25

Eq. (23) then can be solved for S_new, leading to

26

2.3 Model setup and boundary conditions

For the thin-sheet simulations we use the setup of Houseman & England (1986) with a rigid indenter, which is represented by a cosine-shaped velocity function, applied at the front boundary. The same setup is used for the 3-D model, which is extended in the vertical direction (z) by means of four layers (high viscosity upper and low viscosity lower crust; high viscosity upper and low viscosity lower mantle) (Fig. 2), which represent a simplified continental lithosphere (e.g. Burov 2010). Both models either utilize a linear viscous or a power-law viscous rheology.

Figure 2

Viscosity profile with depth used for the linear viscous 3-D model (solid black line) and depth-averaged viscosity profile used for the linear viscous thin-sheet model (dashed black line). Indicated in grey letters are the resulting differential stresses for a strain rate of 10⁻¹⁵ s⁻¹. The vertically integrated stress is 2 × 10¹³ N m⁻¹.

The thin-sheet model consists of a single layer, representing the lithosphere, which possesses depth-averaged parameters (Fig. 2, Table 1). The 3-D model possesses four layers (upper and lower crust, upper and lower mantle), which differ in thickness and material parameters (Fig. 2, Table 1). Model dimensions are 500 km × 500 km in the x- and y-direction, respectively, and 120 km in depth for the 3-D model. Thin-sheet models were calculated with a resolution of 201 × 201 nodes on a desktop PC, whilst the 3-D models were calculated using 48 CPUs with a grid resolution of 129 × 129 × 65 nodes (i.e. a resolution of around 100³ nodes).

Table 1

Densities of crust and mantle are 2700 kg m⁻³ and 3200 kg m⁻³, respectively. All the parameters applied in the thin-sheet model correspond to the depth-averaged parameters of the 3-D model (Table 1) to enable a direct comparison of results. From a chosen differential stress profile (Δσ), together with a thickness (Δz) and a density profile (ρ), a given background strain rate ( s⁻¹), power-law exponents (n= 1, 3) and gravitational constant (g= 9.81 m s⁻¹), the viscosities

27

and the vertically integrated stress of the model lithosphere

28

can be calculated for the 3-D model. Differential stresses (Fig. 2) are chosen in such a way that the vertically integrated stress is in the range of 1.5–2 × 10¹³ N m⁻¹ (2 × 10¹³ N m⁻¹ in the linear viscous setup and 1.64 × 10¹³ N m⁻¹ in the power-law viscous setup), in agreement with estimates for the India–Asia collision system (England & Houseman 1986; England & Jackson 1989).

From the parameters defined for the 3-D model, the averaged quantities for the thin-sheet model are calculated using eq. (11) via using the midpoint rule, thus we have

29

along with the Argand number (eq. 17).

The boundary conditions applied to all simulations are shown in Fig. 3: Indentation is taking place at the front of the model in the y-direction with a depth-constant horizontal velocity of initially 1.5 cm yr⁻¹ (after half of the simulation the indentation velocity is reduced by a factor of 5). v_y is a cosine-function of the x-coordinates and v_x= 0. The lateral (to the left and right of the indenting boundary, v_x= 0) and back (opposite to the indenting boundary, v_y= 0) sides are kept planar and are free slip surfaces, that is, zero normal velocity and zero tangential traction. In the 3-D model the upper boundary is treated as a free surface, (i.e. the total stress normal to the surface is zero), whilst the bottom boundary is a horizontal, free slip surface. This fixed bottom boundary with free slip does not exactly agree with the assumptions of the thin-sheet model, where the bottom boundary is defined as the level of isostatic compensation. However, the bottom layer (lower mantle) of the 3-D model is weak and is therefore able to account for isostatic compensating movements. The maximal variation of σ_zz at the bottom of the entire model is in the order of 900 MPa and with an average density of 3050 kg m⁻³ in the model (Table 1), this maximal stress variation corresponds to 3 km of depth error (i.e. 2.5 per cent of the total thickness of the lithosphere), which makes the comparison of the two models valid.

Figure 3

Model setup and boundary conditions of (a) the thin-sheet model and (b) the 3-D multilayer model. Model dimensions are 500 km × 500 km in x- and y-direction and 120 km in z-direction in the 3-D model. Shortening is applied in the y-direction with a cosine-shaped velocity profile. Stars in (a) indicate the locations for which the evolution of crustal thickness is investigated in Section 3.3. and Fig. 6.

2.4 Performance test

The ability to perform high resolution, 3-D finite element simulations is frequently limited by the choice of the matrix solver used. Often solvers are limited by either excessive memory requirements, or excessive CPU requirements, or both. Here we examine the performance characteristics of a direct factorisation (LU) method compared to several preconditioned iterative methods. The results were evaluated on the basis of maximum memory usage, time required to setup the solver, number of iterations and overall CPU time required to solve the discrete Stokes problem for velocity and pressure. To address this, a number of tests were conducted using the 3-D model setup employed for the 3-D simulations in this study. The tests were conducted on ‘Brutus’, the central high-performance cluster of ETH Zurich and on ‘Achilles’, the cluster of the Computational Tectonics Group at University of Lausanne. ‘Brutus’ consists of AMD Opteron CPUs and compute nodes have 8–128 GB RAM. Also ‘Achilles’ possesses AMD Opteron CPUs and compute nodes have 128 GB or 256 GB RAM.

In our tests, we used three different nodal resolutions (33 × 33 × 17, 65 × 65 × 33, 129 × 129 × 65). The simulations associated with each mesh were executed on 1, 6 & 48 CPUs, respectively (Table 2). In our work we consider a Schur complement reduction (SCR) based solver implemented in LaMEM. This method requires, for each SCR iteration, the solution of a velocity subproblem. The performance of four different solver configurations, used with two element types, Q₁P₀ and Q₂P_–1 was examined. We identify the solvers used for the velocity subproblem as ‘direct’, ‘iterative-ASM’, ‘iterative-AMG’ and ‘iterative-GMG’. The details of each solver configuration are provided below.

Table 2

2.4.1 ‘Direct’

The velocity subproblem was solved using the parallel, multifrontal direct solver MUMPS (Amestoy et al. 2000).

2.4 ‘Iterative-ASM’

The preconditioned Krylov subspace method FGMRES (Saad 2003) was used to solve the subproblem. Here we utilised the domain decomposition preconditioner, ASM (additive Schwarz method; Toselli & Widlund 2005). In our experiments, the ASM preconditioner had an overlap of ten. The application of the inverse operator on each block of ASM, used FGMRES with an incomplete LU (ILU) preconditioner. The iterative scheme on the block problems employed a relaxed stopping condition with a tolerance of 10⁻².

2.4.3 ‘Iterative-AMG’

The same iterative method was used as in the ‘iterative-ASM’ solver configuration, however here we employed an algebraic multigrid (AMG) preconditioner. We used the AMG implementation provided via ML (Gee et al. 2006). AMG is known to be less robust when applied to systems of coupled PDE′s. Consequently, we applied ML to subblocks within K, which were defined via the coupling between like velocity components. The application of the inverse of each subblock consisted of FGMRES with a relaxed stopping condition tolerance of 10⁻². Five iterations of the conjugate gradient (CG) Krylov method preconditioned with block ILU were used as the smoother on each multigrid level. On the coarse grid, a parallel direct solve was performed using MUMPS.

2.4.4 ‘Iterative-GMG’

As in the other iterative approaches, FGMRES was used as the Krylov method. In this solver configuration, we utilised a geometric multigrid (GMG; Briggs et al. 2000). GMG requires a hierarchy of finite element meshes, and a coarse grid operator on each. The coarse grid operators were defined using Galerkin coarsening (Trottenberg et al. 2001). Five iterations of CG preconditioned with block ILU were used as the smoother on each multigrid level, however on the coarsest grid, the direct solver MUMPS was used. For the simulations, the coarse grid contained 17 × 17 × 9 nodes. We used either 2, 3 or 4 grid levels to obtain the three different resolutions used. Each fine grid was obtained by subdividing the elements in each spatial direction by a factor of two.

2.4.5 Performance test results

To conduct our tests, we utilised the boundary conditions, together with the model parameters of our linear viscous model setup (Table 1). Each simulation was run for 10 time steps, over which the average of the evaluated parameters was calculated. The averaged results of the tests performed on ‘Achilles’ are displayed in Table 2, and can be summarized as follows (Results of ‘Achilles’ and ‘Brutus’ are similar, however times reported by ‘Achilles’ are more reliable. Results reported in Chapter 3 were achieved on ‘Brutus’).

1. 1

   Memory usage is always larger for the Q₂P₋₁ elements than for the Q₁P₀ elements and significantly larger when applying the direct solver MUMPS, compared to all other solvers. GMG requires the least amount of memory.

2. 2

   Similarly, MUMPS shows much larger values for the solver setup time, while AMG and GMG perform quite uniformly and fast for the two lower resolutions. The solver setup time for Q₂P₋₁ elements in general is much larger than for the Q₁P₀ elements.

3. 3

   In terms of the CPU time required to solve for velocity and pressure, GMG out performs all other solvers.

4. 4

   The number of iterations needed to solve for pressure was the same for all solvers for a given resolution, which is due to the SCR applied in combination with all solvers.

5. 5

   The solution of velocity shows a less uniform picture in terms of iterations performed. It can be seen that using AMG as a preconditioner gives very robust results, as the number of iterations required for solving for velocity remains constant with increasing resolution.

For the range of mesh resolutions considered here, GMG performs the best with respect to memory usage and CPU time for both the Q₁P₀ and Q₂P_–1 elements. We note that the number of iterations required to solve the velocity subproblem is increasing using GMG, as the mesh resolution and the number of processors increase. We believe the reason GMG does not behave in the same manner as AMG, is because the model contains viscosity jumps on the order of 10² Pas. The interpolation and restriction operators used in GMG are derived directly from the mesh geometry. In this study, we did not align the coarse grid mesh geometry with the jumps in material properties, and thus the interfaces are slightly ‘blurred’ on the coarser grids. This leads to a less robust convergence behaviour compared to AMG. On the other hand, the AMG approach only uses the non-zero structure and coefficients of the matrix to define the coarse grid operators. As expected, this type of algebraic coarsening preserves information related to the jumps in material parameters and thus produces robust convergence. However, this robustness comes at the price of increased CPU time. GMG in combination with the Q₁P₀ elements yields the fastest results, however, experience with the use of LaMEM has shown that Q₂P_–1 elements produce stable pressure fields (non-oscillatory) and therefore only these elements were employed in the 3-D simulations presented in the following sections.

3 Simulations

3.1 Introduction

In the following sections, we present a comparison of deformation styles during continent indentation resulting from a 3-D multilayer model and a thin-sheet model. Two different rheologies were applied (eq. (25)): (1) linear (Newtonian) viscous (n= 1) and (2) power-law viscous (n= 3) rheology.

We performed a number of thin-sheet simulations (not shown here) to confirm that our results agree with those of Houseman & England (1986) and England & McKenzie (1983) concerning crustal thickness evolution for similar combinations of power-law exponents and Argand numbers. For both models and rheologies, two scenarios are investigated: (1) indentation of a rigid indenter with constant indentation velocity and (2) the effect of gravitational equilibration after decreasing the indentation velocity by a factor of 5. Other simulations were performed with a decrease of the indenter velocity by factors of 2 and 10, and the effect of completely stopping indentation was also investigated. We found that a factor 5 decrease in the indentation velocity best represents the outcome of all simulations and therefore these results are discussed in detail here. The effect of slowdown of the indenter is investigated because it is a feature observed for the India–Asia collision (e.g. Molnar & Stock 2009). Furthermore, slowdown of the indenter makes the effect of gravitational collapse, a process taking place throughout the simulation, better visible, as it is less overprinted by a faster indenting velocity. The numerical time step was 0.35 Ma and indentation velocity was reduced after 10 Ma, after which gravitational equilibration was acting for another 10 Ma.

3.2 Overall model behaviour

For the linear viscous setup of the thin-sheet model Ar= 2.9 and the initial crustal thickness is 36 km (Table 1). The indentation of a rigid indenter, simulated by a cosine-shaped velocity profile, leads to an overall thickening of the model domain (Fig. 4a). This thickening is most pronounced at the indenter corner (55 km crustal thickness). At the frontal model boundary, where the indentation velocity reduces to zero, thinning takes place and causes the formation of a depression (30 km crustal thickness). Gravity is acting as a counterforce to thickening during the entire simulation, however, its effect only becomes evident when the indentation velocity is decreased. After this, gravitational collapse visibly equilibrates the differences in elevation and crustal thickness.

Figure 4

Crustal thickness of (a, b) the thin-sheet and (c, d) the 3-D model after 10 Ma of indentation, for (a, c) a linear viscous and (b, d) a power-law viscous (n= 3) rheology. Contour lines are separated by 2.5 km.

The thin-sheet model with power-law viscous rheology uses n= 3 and Ar= 3.6 (Table 1). As in the linear model, overall thickening is observed during indentation (Fig. 4b), although it is more pronounced than in the linear viscous case. Especially the indenter corner exhibits a localized peak (85 km crustal thickness). After 10 Ma indentation, the velocity is decreased, which leads to equilibration of differences in topographic elevation.

During indentation, the results of the linear viscous 3-D simulation (Fig. 4c) are similar to that of the linear viscous thin-sheet model. The 3-D model starts with an initial lithospheric thickness of 120 km, of which the upper 36 km correspond to the crust. As in the thin-sheet model, indentation is simulated by a cosine-shaped velocity profile, which results in shortening in the left half of the modelled domain and to overall thickening. The indenter corner is strongly thickened (55 km crustal thickness) and a depression is formed where the indentation velocity is reduced to zero (30 km crustal thickness, Fig. 4c). In addition to this thickening process, folding of the upper mantle and upper crust is taking place (Fig. 5a). The formation of a fold, exhibiting a long wavelength (200 km), is visible both from the deformed grid in Fig. 5(a) as well as from the distorted (compared to the ones of the thin-sheet model) crustal thickness contour lines (Fig. 4c). The fold amplitude after 10 Ma of indentation is about 1 km in the upper crust and 250 m in the upper mantle. Folding is due to the layered viscosity profile of the 3-D model setup, which responds to shortening with the formation of a folding instability in the more viscous layers. The viscosity ratio between upper crust and lower crust is 80, between upper mantle and lower crust 200 and between upper mantle and lower mantle 100 (Table 1). As soon as the indentation velocity is decreased, the effect of gravitational collapse becomes evident and equilibrates the differences in elevation.

Figure 5

Deformed 3-D grid of (a) the linear and (b) the power-law viscous (n= 3) model after 10 Ma of indentation. Overall thickening and multilayer buckling is taking place. Grey shades on the surface indicate surface elevation.

Models in which n is increased from 1 to 3 result in more pronounced folding and two buckle folds with a smaller wavelength of 100–150 km, although the effective viscosity ratios (for the reference strain rate of 10⁻¹⁵ s⁻¹) are smaller than for the linear viscous simulation (Figs 4d, 5b). The effective viscosity ratio between upper crust and lower crust is 60, between upper mantle and lower crust 170 and between upper mantle and lower mantle 85 (Table 1). Maximum amplitudes of the folds after 10 Ma of indentation are about 3.2 km in the upper crust and 6.8 km in the upper mantle. The maximum dip of the fold limbs for these amplitudes and wavelengths is between 10 and 15° and can likely occur in nature. The overall amount of thickening, on the other hand, is in the same range as in the linear viscous case (maximum crustal thickness was again 55 km; Fig. 4d). After indentation gravitational collapse is taking place.

Wavelengths of folds in the 3-D model depend on the effective viscosity ratio between the layers and, for the values chosen here, are in agreement with typical fold wavelengths published for the crust and lithosphere, which are in the range of a few hundreds of kilometres (e.g. Burov et al. 1993; Nikishin et al. 1993; Burov & Molnar 1998; Caporali 2000). In addition, these wavelengths also agree with scaling laws published for 2-D folding (Biot 1961; Fletcher 1974; Schmalholz et al. 2002).

In synthesis, three deformation processes can be distinguished when observing the overall evolution of crustal thickness in the 3-D models: (i) overall thickening and thinning, (ii) folding of competent layers (upper crust and upper mantle) and (iii) compensating movements, which take place in the less competent layers of lower crust and lower mantle. Of these processes, thickening and compensating movements, that is, gravitational collapse, are also present in the thin-sheet model. Thickening develops accordingly in the 3-D and thin-sheet model for a linear viscous setup, however, differences in crustal thickness between the two models are more pronounced when a power-law rheology is utilised.

3.3 Evolution of crustal thickness at four selected points

The evolution of crustal thickness was investigated at four selected points (Fig. 3a). Point 1 is situated at the indenter corner, point 2 lies in the forming depression, point 3 in the area that develops buckle folds and point 4 is positioned roughly in the centre of the model. As expected, the thin-sheet model is a good approximation to the 3-D model in terms of overall crustal thickening—particularly during indentation. Crustal thicknesses of the linear viscous 3-D and thin-sheet models develop similarly (Fig. 6a). Crustal thickness increases (at points 1, 3 and 4) or decreases (at point 2) in an approximately linear manner. In the 3-D model after 10 Ma, that is, after slowdown of indentation, crustal thicknesses at the selected points develop in the same manner (i.e. linearly), albeit more slowly, that is, they keep on thickening (at points 1, 3 and 4) or thinning (at point 2). In the thin-sheet model, these large crustal thicknesses are not stable and collapse. While points 2 and 4 are showing a linear increase in crustal thickness, point 1 shows decreasing thickness, whilst the thickness at point 3 remains approximately stable. The crustal thicknesses at the different points slowly converge to an average crustal thickness of around 45–50 km, as the averaged viscosity of the thin-sheet is too weak to sustain gravitational forces. In the rheologically layered 3-D model, the crustal thickness is not decreasing because the competent layers can support the gravitational load and crustal thickening continues. The viscosity of the upper crust and upper mantle in the 3-D model is larger than the depth-averaged viscosity in the thin-sheet model (Fig. 2) and, therefore, the 3-D model can support higher gravitational stresses than the thin-sheet model because the lithosphere in the 3-D model includes thin but strong layers.

Figure 6

Development of crustal thickness for four selected locations (Fig. 3a) for (a) the linear viscous and (b) the power-law viscous (n= 3) thin-sheet (TS) and 3-D models.

With a power-law viscous rheology, the thin-sheet and 3-D model again behave similarly during indentation, with an approximately linear increase in crustal thickness at points 2 and 4 and an approximately exponential growth at points 1 and 3. The latter are situated in the area of the model, which experiences the strongest deformation, thus explaining the difference to points 2 and 4. At the corner (point 1), the thin-sheet model is extremely thickened (Fig. 6b), while thickening in the power-law viscous 3-D model is in the same range as in the linear viscous case. After slowdown of indentation, the thin-sheet model keeps thickening (except at point 2, where the thickness continues to decrease) and does not collapse in contrast to the result for the linear viscous rheology. This can be explained by the effective power-law viscosity, which increases when the strain rate is reduced (and decreases when the strain rate is increased). As pointed out by England & McKenzie (1983), there is an upper limit of crustal thickness which can be supported in a thin-sheet model for Argand numbers >0. The Argand numbers applied in our simulations are 2.9 and 3.6. In this case, large crustal thicknesses can only be supported if the lithosphere has a power-law rheology, and a linear viscous setup results in fast collapse. Houseman & England (1986) also identified that power-law exponents should be larger than 3 for Argand numbers between 1 and 10, to generate crustal thicknesses similar to those of the Tibetan Plateau.

3.4 Quantification of folding along a selected profile

As already mentioned in Section 3.2, folding is one of the processes taking place during indentation in the mechanically layered 3-D models. To quantify this process, fibre strain rates along vertical lines are investigated. In a folded layer fibre strain rates deviate from the average strain rate, calculated along the vertical line, due to relative extension, respectively compression in the hinge zone of the fold (Fig. 7a–d, e.g. Schmalholz & Podladchikov 2000). In our model, these vertical lines lie in a section (y, z) through the middle of the indenter, one line perforating the synform, the other the antiform (Fig. 7a–b). The strain rates considered for this quantification are , as direction of indentation is along the y-axis.

Figure 7

Fibre strain rates at the antiform of the linear viscous 3-D model and at the synform of the power-law viscous 3-D model: (a) Profile in y-direction through the centre of the indenter with locations of lines through the syn- and antiform of the linear viscous 3-D model. Arrow denotes antiform for which fibre strain rates are shown in (c). (b) Profile in y-direction through the centre of the indenter with locations of lines through the syn- and antiform of the power-law viscous 3-D model. Arrow denotes synform for which fibre strain rates are shown in (d). (c), (d) Fibre strain rates in the antiform of the linear viscous model and in the synform of the power-law viscous model. Dashed line indicates flip between extensional and compressional fibre strain rates relative to the depth-averaged strain rate. Strong upper crust (UC) and upper mantle (UM) show the linear increase, respectively decrease of fibre strain rates in a folded layer. In (d) is shown how the maximum difference between fibre strain rates in the more viscous upper crust and upper mantle (max) is calculated for (e) and (f). (e), (f) Evolution of difference between fibre strain rates over time, for the first 10 Ma of indentation.

Fig. 7(c)–(d) shows fibre strain rates—the deviation of strain rates from the average value—along the vertical line through the antiform of the linear viscous 3-D model (Fig. 7a) and through the synform of the power-law viscous 3-D model (Fig. 7b) after 10 Ma of indentation. Fibre strain rates are illustrated in terms of relative differences in percent with respect to the average strain rate (). The more viscous upper crust and upper mantle show a linear increase (antiform) and decrease (synform), respectively, of the fibre strain rates from the top of the layer to the bottom.

Fibre strain rates in the synform of the upper crust of the power-law viscous model are larger than the average due to compression taking place at the inner arc of the developing synform (Fig. 7d). Fibre strain rates decrease with depth and equal the average strain rate at the neutral line (zero strain, e.g. Ramsay & Huber 1983). Also, it is clearly visible that the neutral line does not necessarily have to be located in the middle of a layer. Fibre strain rates reach a minimum at the outer arc of the synform where extension is taking place. In the lower crust fibre strain rates show a parabolic pattern, indicating enhanced flow in the lower crust. The same pattern as in the upper crust can be observed in the upper mantle. The weak lower mantle shows increasing fibre strain rates with depth. At the antiform in the linear viscous model, the same behaviour can be observed, but reversed, that is, extension at the outer arc and compression at the inner arc of the folded layers. However, in the power-law viscous model fibre strain rates cover a much larger range. Deviations from the average strain rate range between –150 per cent and 200 per cent, whereas in the linear viscous case differences range from –50 per cent to 20 per cent.

In Fig. 7(e)–(f) the evolution of the maximum difference (max) between fibre strain rates within the upper crust and upper mantle are shown for the first 10 Ma of indentation. The latter half of the simulation, when indentation velocity is reduced by a factor of 5, was not considered to be able to separate the effects of buckling and compensating movements. The behaviour of both an antiform and a synform is investigated for both the linear and power-law viscous model. In the linear viscous model, for both layers (upper crust and upper mantle) and fold types (antiform and synform), the maximum difference between fibre strain rates increases with continued indentation. Upper crust and upper mantle behave very similarly. After an initial phase with very little differences (0–5 per cent) between fibre strain rates, when presumably thickening is predominant, differences in the fibre strain rate increase after 6 Ma of indentation, with a stronger increase (up to 30 per cent) in the antiform. In the power-law viscous model, differences between the maximum and minimum fibre strain rates in upper crust and upper mantle are much larger (up to 280 per cent) than in the linear viscous case, that is, buckling is much stronger. As in the linear viscous case, differences between fibre strain rates increase significantly after 6 Ma of indentation. Only the antiform of the upper mantle shows increased values at an earlier stage and reaches the largest value with 280 per cent difference between fibre strain rates within the layer. This also agrees with the larger amplitude generated in the upper mantle of the power-law viscous model (Section 3.2). Unlike in the linear viscous model, upper crust and upper mantle show a different behaviour, that is, differences between fibre strain rates are larger in the upper mantle. Thus in the power-law viscous model, the upper mantle is the dominant layer in terms of buckling.

3.5 Horizontal velocities in the thin-sheet model and layers of the -D model

The horizontal velocity fields of different layers of the 3-D model and the horizontal velocities of the thin-sheet model were compared to identify differing flow patterns which allow identification of the differences in the deformation style. Horizontal velocity fields are displayed (Figs 8 and 9) with white velocity vectors, composed of the x- and y-components of the velocity field, and with the corresponding overall magnitude () as filled contour plots. Furthermore, vertical profiles illustrating the velocities in upper and lower crust and upper mantle of the 3-D models are shown (Fig. 10) to compare velocities in the different layers. As before, only the results of simulations with a slowdown of indentation velocity by a factor 5 are discussed.

Figure 8

Comparison of horizontal velocities of the linear viscous thin-sheet model (top left-hand side) with the surface (top right-hand side), the lower crust (bottom left-hand side) and the upper mantle (bottom right-hand side) velocities of the 3-D model after (a) 10 Ma of indentation and (b) shortly after slowdown of indentation velocity by a factor of 5. White vectors indicate the horizontal velocity field, while the filled contours indicate the overall magnitude of the horizontal velocities ().

Figure 9

Comparison of horizontal velocities of the power-law viscous (n= 3) thin-sheet model (top left-hand side) with the surface (top right-hand side), the lower crust (bottom left-hand side) and the upper mantle (bottom right-hand side) velocities of the 3-D model after (a) 10 Ma of indentation and (b) shortly after slowdown of indentation velocity by a factor of 5. White vectors indicate the horizontal velocity field, while the filled contours indicate the overall magnitude of the horizontal velocities ().

Figure 10

Velocities of the linear and power-law viscous 3-D models after 10 Ma of indentation. The three uppermost layers [upper crust (UC), lower crust (LC) and upper mantle (UM)] are shown: (a) north–south profile of the linear viscous 3-D model, (b) west–east profile of the linear viscous 3-D model with lower crustal flow, (c) north–south profile of the power-law viscous 3-D model with pronounced folding of the upper mantle, (d) west–east profile of power-law viscous 3-D model with strong lower crustal flow. Position of profiles is shown in inset in (d). Downward pointing velocity vectors close to the indenter, respectively upward pointing velocity vectors away from the indenter are due to compensating movements invoked by gravitational forces equilibrating differences in topographic elevation.

In Fig. 8(a) the results of models with a linear viscous rheology after 10 Ma of indentation, with a velocity of 1.5 cm yr⁻¹, are displayed. The horizontal velocity fields at the surface and upper mantle of the 3-D model (Fig. 8a, top right-hand side and bottom right-hand side) and the thin-sheet model (Fig. 8a, top left-hand side) are very similar. Velocity vectors indicate a rotation around the indenter corner. Only minor differences occur within the lower crust (Fig. 8a, bottom left-hand side).

Once the indentation velocity is reduced by a factor of 5 (from 1.5 cm yr⁻¹ to 3 mm yr⁻¹), differences between thin-sheet and 3-D model are more distinct. The velocities from thin-sheet models are generally larger than those from the 3-D models, except at the indenter front where the same boundary condition was applied (Fig. 8b, top left-hand side). The lower crust shows differently shaped contours (Fig. 8b, bottom left-hand side) in the area of the buckle fold, that is, velocities in the order of 2–3 mm yr⁻¹ cover a larger area in front of the indenter. Also, the velocity vectors indicate a stronger rotation around the indenter corner than on the surface and in the upper mantle. In the lower crust, velocity vectors in the area of the indenter corner are larger and oriented more strongly to the East (i.e. the right side of the model). Thus, the lower crust flows relative to the upper crust and upper mantle, which is also visible in West–East (W–E) vertical profiles (Fig. 10a and b).

Compared to simulations with a linear viscous rheology, differences between thin-sheet and 3-D model are more pronounced if a power-law rheology is applied (Fig. 9a). Surface and upper mantle velocity contours are further apart than in the thin-sheet model (Fig. 9a, top right-hand side and bottom right-hand side), which can be attributed to the formation of two buckle folds in the competent layers (upper crust and upper mantle) of the 3-D model during indentation. Especially within the lower crust, the effect of these buckles on the velocity field becomes apparent (Fig. 9a, bottom left-hand side), and velocities are larger in the area of the antiforms of the buckle folds and smaller in the area of synforms. The orientation of velocity vectors is very similar in the thin-sheet and 3-D models and in comparison to the linear viscous simulation, flow around the indenter corner is less pronounced.

After a reduction in indentation velocity, the free surface and upper mantle again show distorted contours in the area of the buckle folds (Fig. 9b, top right-hand side and bottom right-hand side). In the lower crust the pattern of the buckle folds with two antiforms and a synform is clearly visible and velocity vectors show a stronger clockwise rotation around the indenter corner (Fig. 9b, bottom left-hand side). The pattern of the velocity field from the thin-sheet model does not significantly alter after slowdown of indentation (Fig. 9b, top left-hand side). We also note that the flow direction does not greatly differ between the free surface and the upper mantle obtained from the 3-D and the thin-sheet models. However, in the lower crust, velocity vectors are not aligned with those obtained from the thin-sheet model, particularly in the area of the indenter corner. This is also visible in vertical profiles (Figs 10c and d). In the area of the indenter corner (at position 150–200 km in the W–E profile) the horizontal component of the velocity vector in the lower crust is significantly larger than in the upper crust and upper mantle.

Comparing the results from the linear (Fig. 8) and power-law viscous rheologies (Fig. 9), one can observe larger differences between the 3-D and thin-sheet models when a power-law viscous rheology is applied. This is due to stronger multilayer folding and the following stronger equilibrating movements, after the reduction of the indentation velocity. Yet, our results also indicate that the thin-sheet model is a good approximation of a full 3-D model, particularly during indentation with constant velocity. The large-scale patterns and magnitudes of the horizontal velocity fields agree well. The free surface and upper mantle of the 3-D model basically show the same flow pattern and are very similar to the averaged thin-sheet velocities.

While in Section 3.2 the lack of buckling in thin-sheet models became apparent by comparing crustal thicknesses, the comparison of horizontal velocity fields of the thin-sheet model and different layers of the 3-D model, reveals the mechanically decoupled character of these layers with different strength, which can deform relative to each other. In addition, different processes are dominant within the different layers: The competent upper crust and upper mantle buckle, whereas the less viscous lower crust is deforming by thickening and lower crustal flow relative to its bounding layers.

3.6 Quantification of crustal flow

To quantify lateral displacement of material in the crustal layers of the different models, an approach introduced by Houseman & England (1993) is used in an adapted form. If zero strain in x-direction (perpendicular to the direction of indentation) and conservation of crustal volume is assumed, the amount of crustal thickening in the model directly in front of the indenter should correspond to the reduction of crustal volume by the indenter

30

Here A represents the indented area and S₀ is the initial thickness, which is 36 km for the entire crust and 18 km for upper and lower crust. S(x,y) corresponds to the final crustal thickness. The indented volume (A·S₀) should be equivalent to the volume added by crustal thickening (). However, as there is flow in the x-direction, the area reduced by the indenter is usually not equal to the area added directly in front of it.

To quantify the amount of crustal flow using eq. (28), two regions of the model are investigated (Fig. 11): (i) the indentation front, where the indentation velocity is constant (A1) and (ii) the transitional part of the frontal boundary, where the velocity decreases to zero (A2). The reduction of area by the two parts of the indenter (A1-I and A2-I) is compared with the area added inside the model (A1-C and A2-C for the entire crust, A1-UC and A2-UC for the upper crust, A1-LC and A2-LC for the lower crust, Fig. 11) after 20 Ma of indentation, that is, at the end of the simulation. Differences between the area reduced by the indenter and the area added in the model are reported in per cent in Tables 3 and 4 (dA1-C, dA2-C, etc.).

Figure 11

Areas distinguished for estimating crustal flow. The indented volume (A1-I and A2-I) should be equivalent to the volume added by crustal thickening (A1-C/A1-UC/A1-LC and A2-C/A2-UC/A2-LC). However, as there is flow in the x-direction, the area reduced by the indenter is usually not equal to the area added directly in front of it.

Table 3

Table 4

Comparing areas calculated for the entire crust (A1-C, A2-C) in the thin-sheet and the 3-D models with the actual reduction of area by the indenter (A1-I, A2-I, Table 3), the crustal areas added in front of the indenter are smaller. This means that lateral crustal flow away from the indenter front is taking place. Comparing linear and power-law viscous thin-sheet models, differences between area reduction by the indenter (A1-I, A2-I) and the area added in front of it are smaller with a power-law viscous setup, that is, dA1-C and dA2-C are larger in the power-law viscous model. This implies that lateral crustal flow is less strong with a power-law viscous rheology. For both rheologies, dA2-C is larger than dA1-C, indicating stronger lateral displacement of crustal material away from A1-C. The same observations as for the thin-sheet model can be made for the crust of the 3-D model. In general, the respective thin-sheet and 3-D models correlate well and show similar values for the differences between reduced and added area (dA1-C, dA2-C). dA1-C corresponds very well for the linear viscous thin-sheet and 3-D models and as well for the power-law viscous thin-sheet and 3-D model. dA2-C only has similar values for the linear thin-sheet and 3-D model. However, in the power-law viscous thin-sheet and 3-D model the values for dA2-C differ strongly, with less lateral flow in the 3-D model.

Comparing the indented and added areas of upper and lower crust separately in the 3-D model (Table 4), reveals that differences between indented and added area generally are considerably smaller in the upper crust (dA1-UC > dA1-LC, dA2-UC > dA2-LC), that is, more material is laterally displaced in the lower crust. Differences between the indented and added A1 (dA1-UC, dA1-LC) area are twice as large in the lower crust, that is, twice as much material is laterally displaced away from A1 in the lower crust than in the upper crust. Differences between lateral flow in the upper and lower crust get even larger for A2, where in the case of the power-law viscous 3-D model the area added in front of the indenter exceeds the amount of area indented. For power-law viscous rheology around six times more lower crustal material is moved laterally than upper crustal material (Table 4). As for the entire crust, crustal flow in the upper and lower crust of the power-law viscous model is less than in the linear viscous model.

In summary, comparing flow in the upper and lower crust separately revealed, that flow in the weaker lower crust is significantly larger than in the stronger upper crust and differences in material flow between upper crust and lower crust are larger for power-law viscous rheologies.

3.7 Second strain rate invariant in -D model

The second strain rate invariant (i.e. second invariant of the strain rate tensor) is a representative measure for the 3-D strain rate field and is displayed in North–South and East–West slices through the model. This calculation is performed after 10 Ma of indentation and shortly after the onset of slowdown by a factor of 5, for both the linear viscous and the power-law viscous setup.

During indentation, strain rates are largest at the indenter corner and generally in the less competent lower mantle and lower crust (Figs 12a–b, 13a–b). This can be attributed to the strong thickening, taking place in these layers. In addition, isostatic compensating movements take place in these parts of the model. Strain rates generally decrease with increasing distance from the indenter front and corner. In the lower mantle high strain rates reach furthest to the East (i.e. to the right), which indicates that the lower mantle deforms due to isostatic compensating movements and material evades from below the strongly thickening western part of the model. The lowest strain rates occur in the competent upper crust and upper mantle, which deform by both thickening and buckling during the overall simulation. Strain rates vary between one and two orders of magnitude.

Figure 12

The second strain rate invariant of the 3-D model is shown on (a, c) north–south and (b, d) east–west slices for the linear viscous model setup for (a, b) the situation after 10 Ma of indentation and (c, d) shortly after slowdown of indentation by a factor of 5.

Figure 13

The second strain rate invariant of the 3-D model is shown on (a, c) north–south and (b, d) east–west slices for the power-law viscous (n= 3) model setup for (a, b) the situation after 10 Ma of indentation and (c, d) shortly after slowdown of indentation by a factor of 5.

After reducing the velocity of indentation, the differences in strain rate between the competent and less competent layers increase (Figs 12c–d, 13c–d), with the competent layers showing very low strain rates. Lower crust and lower mantle accommodate the gravitational equilibrating movements resulting from the differences in elevation, which is indicated by larger strain rate values. High strain rates in the lower crust indicate lower crustal flow, taking place relative to the above and below lying competent layers. Note, however, that these high strain rates are not equally distributed (e.g. Figs 12c–d), with lower crustal flow being especially pronounced at the indenter corner and in the synforms and antiforms of folds. In the power-law viscous model, high strain rates are more smoothly distributed (Figs 13c–d), but concentrated on the lower half of the lower crust. Conversely, high strain rates in the lower mantle are situated in the upper parts of this layer and do not cover an area as large as during faster indentation.

Differences in the second strain rate invariant between the different layers reveal their dominant deformation process. Less competent layers, like lower mantle and lower crust, are characterised by larger strain rates and mainly deform by thickening and compensating (flow) movements. Competent layers on the other hand exhibit lower strain rates and deform by thickening and buckling. The main difference between the power-law (Fig. 13) and linear viscous (Fig. 12) models, apart from the amount of buckling, is that strain rates vary over two orders of magnitude in the power-law viscous model, whilst in the linear viscous model strain rates vary between 1 and 1.5 orders of magnitude.

4 Discussion

Thin-sheet models assume that the lithosphere is mechanically coupled and deforms as a continuous, uniform sheet. Some studies (e.g. Flesch et al. 2005; Wang et al. 2009) compared surface GPS velocities with SKS shear wave splitting data assigned to the lithospheric mantle and suggested that the data indicate a coherent deformation in the crust and lithospheric mantle. However, while surface GPS velocities represent an instantaneous deformation field, the anisotropy revealed by shear wave splitting data is attributed to the lattice preferred orientation of olivine crystals in the mantle, which in turn correlates with the major axis of the finite strain ellipsoid, and thus the integrated deformation field. This means, GPS and SKS shear wave splitting data reflect deformation fields of different time frames and are as such a questionable choice for comparison. Furthermore, the question under what terms the anisotropy of shear waves can be considered as an indicator for the direction of flow of the mantle is still debated (Kaminski & Ribe 2002). This question could be resolved in future, by combining 3-D models like the one presented with texturing models (e.g. D-Rex, Kaminski et al. 2004). Nevertheless, the similarity between surface GPS and SKS shear wave splitting data is used by, for example, Flesch et al. (2005) and Wang et al. (2009) as evidence for a mechanically coupled lithosphere in a large area of the India–Asia collision. In contrast, our 3-D simulation results indicate that the upper crust and upper lithospheric mantle can exhibit similar horizontal velocity fields and, hence, a coherent deformation, even though upper crust and lithospheric mantle are mechanically decoupled by an intermediate weak lower crust (Figs 8 and 9). The mechanical decoupling causes, for example, the significant lower crustal flow relative to the upper crust and upper mantle. In the presented simulations, the effective viscosity of the lower crust is about two orders of magnitude smaller than the effective viscosities of the upper crust and upper lithospheric mantle (Fig. 2). The upper crust and the mechanically decoupled upper lithospheric mantle deform coherently during both deformation phases, that is, indentation and gravitational equilibration (Figs 8 and 9). Therefore, similar patterns in surface GPS and mantle SKS shear wave splitting data do not necessarily indicate a mechanically coupled lithosphere and are unsuitable as significant evidence for it. Consequently, the fundamental assumption of the thin-sheet model, that lithospheric deformation can be considered as deformation of a continuous viscous sheet, seems questionable. Nevertheless, for predicting the overall thickness evolution (Figs 4 and 6) and the horizontal velocities of the competent layers in the 3-D multilayer model (Figs 8 and 9), the thin-sheet model is a good first-order approximation and results largely coincide with the 3-D model. However, for a more accurate quantitative prediction of, for example, strain rates (and hence stresses) and velocities, 3-D multilayer models are required because they consider processes such as buckling and lower crustal flow and these processes cause strain rate and velocity variations in the same order (and even more) than the strain rates and velocities resulting from diffuse thickening.

The intensity of buckling can be quantified by investigating fibre strain rates (Fig. 7). In the thin-sheet model strain rates are per definition constant with depth, however, in the mechanically layered 3-D model, deviations from an average value can be quantified and indicate the intensity of buckling. With a systematic investigation of fibre strain rates in the 3-D model for different rheological parameters, it would be possible to define a threshold marking the transition between thickening and buckling. Observations made for the performed 3-D simulations indicate, as expected form the analytical solutions (e.g. Schmalholz et al. 2002), that folding is stronger with a power-law viscous rheology, that is, differences between fibre strain rates cover a much larger range (ca. 300 per cent) than in the linear viscous model (ca. 30 per cent). In the power-law viscous model, where differences between fibre strain rates are much larger, significant buckling is taking place. Here, the fold amplitude and the difference between fibre strain rates are largest in the upper mantle, which is the strongest layer in the applied model setup and is dominating the buckling-deformation (Table 1). Also, the thickness of the upper mantle and upper crust remains nearly constant whereas the thickness of the weak lower crust varies significantly because it must accommodate the deformation resulting from the different size of the fold amplitudes in the upper mantle and upper crust (Fig. 7b). The quantification of buckling during lithospheric shortening is important because buckling causes structural (geometric) softening and can have an important impact on the localization of strain (e.g. Schmalholz ., 2005). If thermomechanical coupling is considered then the strain localization around buckle-folds can lead to thermally activated shear zone formation which was studied in 2-D numerical simulations of lithospheric shortening (e.g. Burg and Schmalholz, 2008; Schmalholz ., 2009).

Common patterns between the 3-D and thin-sheet model can be identified from the comparison of horizontal velocities of the two models. Thin-sheet models are visually an acceptable approximation of the 3-D models during indentation. Nevertheless, differences in the velocity field between the 3-D and thin-sheet model become large as soon as gravitational collapse becomes significant and a power-law rheology is applied. Whilst the velocity fields in the upper crust and upper mantle of the thin-sheet and 3-D model are in close agreement, the lower crust shows a different pattern due to lower crustal flow. Lower crustal flow is induced by the slowdown of indentation and the enhanced effect of gravitational collapse, leading to an equilibration of differences in elevation, taking place in the less competent layers, such as the lower crust (Fig. 10). The 3-D results therefore clearly reveal additional processes, which are important during continent indentation, such as the decoupled behaviour of the different layers of the lithosphere with magnitudes of velocities and directions of the horizontal velocity field showing gradients with depth. Corresponding flow fields on the surface and in the upper mantle, therefore do not rule out the existence of a decoupling lower crust.

The crustal thicknesses of 3-D and thin-sheet model are comparable (Figs 4–6). The temporal evolution of crustal thickness for both models coincides, particularly during indentation and with a linear viscous rheology, generating comparable crustal thicknesses. However, around the indenter corner in the power-law viscous thin-sheet model thickening is unrealistically strong and deviates from thickening in the 3-D model. After slowdown of the indentation velocity, thicknesses in the linear viscous thin-sheet model decrease, whilst both the 3-D models and the power-law thin-sheet model continue to thicken or to thin. Linear viscous thin-sheet models are not able to support large topographic differences and quickly collapse when indentation velocity decreases, due to lacking stabilizing high viscosity layers. And even though the linear viscous models exhibit a slightly higher effective lithospheric strength than the power-law viscous models, strong thickening in thin-sheet models is only feasible with a power-law rheology.

With regard to lateral flow of the entire crust (Table 3), the results from the thin-sheet and 3-D models agree well. When using a linear viscous rheology, the thin-sheet model approximates the results of the 3-D model rather well. However, the results of the 3-D model show that considerably more material of the lower crust is moved than of the upper crust (Table 4). Therefore, if a relatively weaker lower crust exists, the thin-sheet models significantly overestimate flow of the upper crust while significantly underestimating flow of the lower crust.

Considering the aforementioned similarities between thin-sheet and 3-D model, we consider the thin-sheet model an acceptable first-order approximation to investigate continent indentation, especially when considering large scale problems of thousands of kilometres (the impact of indentation on a large scale orogenic system was also pointed out by Ellis (1996)). Thin-sheet models agree best with their corresponding 3-D multilayer model when a linear viscous rheology is applied and during indentation with a constant velocity. Crustal thickening and the horizontal velocity field are well reproduced. However, a mechanically layered 3-D setup does make a difference and processes like buckling and lower crustal flow should not be ignored when discussing continent–continent collision, as they can significantly influence geometry, strain rates and stresses, especially on smaller scales in the order of hundreds of kilometres around the collision zone.

Garthwaite & Houseman (2011) reach conclusions qualitatively similar to ours with their recent comparison of a 2-D thin-sheet and a 3-D model. However, they use a 3-D model, which is mechanically homogeneous and not layered, hence preserving many of the fundamental simplifying assumptions of the original thin-sheet model. As in this study, they found that differences between thin-sheet and 3-D models are large at the indenter corner and with increasing power-law exponent. Furthermore, they accordingly conclude that thin-sheet approximations are well suited for large-scale problems. In comparison to the 3-D model of Garthwaite & Houseman (2011) and to all the thin-sheet models applied to the India–Asia collision, the main advantage of our 3-D model is that it does not exclude buckling and lower crustal flow a priori. From our study, we observe that even from a very simple rheological setup with only four layers with different material parameters, a range of deformation mechanisms can be observed. We showed that apart from large strain indentation and differential thickening—also observed in the thin-sheet model—additionally significant multilayer buckling and lower crustal flow are taking place in the 3-D model. Both multilayer folding and lower crustal flow have been suggested based on geophysical data and geological observations (e.g. Royden 1996; Cloetingh et al. 1999; Beaumont et al. 2001; Shin et al. 2009) and 3-D models such as the one presented allow quantification of these mechanisms. As some of these mechanisms are often taking place simultaneously, models focusing on one specific process lack information related to several deformation processes, and thus their relative importance during continental indentation cannot be assessed.

We only performed simulations for a power-law exponent of three and the simulations show that 3-D results for power-law viscous rheology deviate stronger from the corresponding thin-sheet results than 3-D results for a linear viscous rheology. Also, Garthwaite & Houseman (2011) found that differences between thin-sheet and 3-D models increase with increasing power-law exponent. Effective power-law exponents of flow laws describing the deformation of the crust and mantle are likely larger than three (e.g. Sonder & England 1986), at least in some depth range such as the upper lithospheric mantle. For example, the deformation of the upper mantle can be described by a high stress flow law (based on Peierl′s mechanism) for olivine (e.g. Goetze 1978) and this high stress flow law can be transformed into the standard form of a power-law flow law (Schmalholz & Fletcher 2011). The corresponding effective power-law exponents are likely between 10 and 20 for natural deformation conditions of the upper lithospheric mantle, especially at the onset of collision. This means that for natural conditions processes such as buckling and lower crustal flow are likely first-order processes around the collision zone within a spatial scale corresponding to the order of the thickness of the lithosphere, that is, a few hundreds of kilometres. Consequently, 3-D multilayer models are preferable to thin-sheet models when studying the tectonic evolution of the continental collision zone.

The thickness evolution and the horizontal velocity field in our simulations is strongly controlled by the applied indentation boundary conditions and fundamental differences between thin-sheet and 3-D multilayer models are not really expected. For other model setups, closer to natural conditions, such as for example with a deforming indenter or with laterally varying strength within the model domain or a more complex rheology, larger differences between thin-sheet and 3-D models can be expected.

This study focuses on viscous deformation processes, which implies that many other potentially important effects during continental collision are ignored. Additional features, such as plasticity, or a more sophisticated geometry, such as including a slab, are deliberately dispensed with, to be able to compare our 3-D model with the original thin-sheet models of England & Houseman (1986) and Houseman & England (1986). The aim is to understand the first-order effects of including the third dimension, that is, the layered nature of the lithosphere, and to study which additional processes arise. The main scope of this study is a theoretical comparison of thin-sheet models and their equivalent 3-D multilayer models in terms of viscous deformation processes. We were able to show that 3-D models provide a more accurate and complete picture of continental collision than thin-sheet models. 3-D models allow studying the timing, locality and relative importance of different processes, which take place simultaneously on a scale of hundreds of kilometres.

5 Conclusions

We performed a comparison of a numerical 3-D multilayer model with its corresponding thin-sheet equivalent for the scenario of continent indentation. The 3-D model employs four layers with different rheological properties representing strong upper and weak lower crust and strong upper and weak lower lithospheric mantle. Two deformation phases (indentation and gravitational collapse) and two rheologies (linear and power-law viscous) were considered. Differences between the 3-D and thin-sheet models are significant (i) during gravitational collapse, after the indentation velocity slowed down by a factor 5, for both rheologies and (ii) for both deformation phases when a power-law viscous rheology was used. Buckling is significant and the related spatial variations in strain rate are up to two orders of magnitude around the collision zone and indenter corner within a distance on the order of the lithospheric thickness (i.e. hundreds of kilometres) for a power-law viscous rheology. Lateral mass flow of weak lower crustal material is also significant and can be up to six times larger than lateral flow of upper crustal material. Both buckling and lower crustal flow are ignored in thin-sheet models. Nevertheless, thin-sheet models approximate the horizontal velocity field, the amount of overall crustal thickening and overall crustal flow quite well, especially during indentation with a linear viscous rheology. Thin-sheet models are, therefore, a good first-order approximation to model lithospheric deformation on the large-scale, that is, the scale considerably larger than the thickness of the lithosphere (thousands of kilometres), and away from the collision zone for linear viscous rheologies.

3-D multilayer models are necessary to better understand lithospheric deformation and to better interpret geophysical data. For example, our 3-D multilayer model shows that upper crust and upper lithospheric mantle can deform coherently although being mechanically decoupled by an intermediate weak lower crust (with significant lower crustal flow) which questions the usage of similarities between surface GPS and mantle SKS shear wave splitting data as evidence for a mechanically coupled lithosphere. 3-D multilayer models not imposing a certain viscous deformation style a priori are needed to better quantify the deformation processes during continent indentation and to assess their relative importance especially around the collision zone and indenter corners.

Acknowledgements

We thank Sergei Medvedev, Susan Ellis and Stephane Labrosse for their constructive review of the manuscript. We thank ETH Zurich for the use of ‘Brutus’ and the University of Lausanne for the use of ‘Achilles’. SML and SMS acknowledge support from ETH grant 0–20497-08. DAM was financially supported by the ETH Zurich Postdoctoral Fellowship Program and BK by ERC Starting Grant MODEL-258830.

Appendix

Here, the detailed operations performed to get from eq. (13) to eq. (14) are described. Eq. (13) is given as

31

where p*_c and p*_m are the pressure in the crust and in the mantle, respectively:

32

Inserting eq. (A2) into (A1) leads to

33

To facilitate the evaluation of the integrals in eq. (A3), variables are simplified and the integral boundaries are changed accordingly. For the integration of mantle pressure, the mantle thickness h_m= L−b−z is introduced, the variable of integration thus becomes dh_m=−dz′ and the limits of integration go from L−b to 0. Integration of the crustal pressure is done accordingly: The crustal thickness is h_c= L+h−z and the integration variable changes to dh_c=−dz′. The limits of integration are from h+b to 0. These transformations lead to the following expression:

34

Evaluating the integrals leads to

35

For the normalization it is assumed that L≫h, and for simplifying eq. (A5) the relations h=S−b and and eq. (9) are used. This leads to the final equation for , corresponding to eq. (14):

36