Control of 3-D tectonic inheritance on fold-and-thrust belts: insights from 3-D numerical models and application to the Helvetic nappe system

We apply three-dimensional (3-D) thermo-mechanical numerical simulations of the shortening of the upper crustal region of a passive margin in order to investigate the control of 3-D laterally variable inherited structures on fold-and-thrust belt evolution and associated nappe formation. We consider tectonic inheritance by employing an initial model configuration with basement horst and graben structures having laterally variable geometry and with sedimentary layers having different mechanical strength. We use a visco-plastic rheology with a temperature-dependent flow law and a Drucker–Prager yield criterion. The models show the folding, detachment (shearing off) and horizontal transport of sedimentary units, which resemble structures of fold and thrust nappes. The models further show the stacking of nappes. The detachment of nappe-like structures is controlled by the initial basement and sedimentary layer geometry. Significant horizontal transport is facilitated by weak sedimentary units below these nappes. The initial half-graben geometry has a strong impact on the basement and sediment deformation. Generally, deeper half-grabens generate thicker nappes and stronger deformation of the neighbouring horst, while shallower half-grabens generate thinner nappes and less deformation in the horst. Horizontally continuous strong sediment layers, which are not restricted to initial graben structures, cause detachment (décollement) folding and not overthrusting. The amplitude of the detachment folds is controlled by the underlying graben geometry. A mechanically weaker basement favours the formation of fold nappes, while stronger basement favours thrust sheets. The model configuration is motivated by applying the 3-D model to the Helvetic nappe system of the Central Alps of France and Switzerland. Our model reproduces several first-order features of this nappe system, namely (1) closure of a half-graben and associated formation of the Morcles and Doldenhorn nappes, (2) overthrusting of a nappe resembling the Wildhorn and Glarus nappes, and (3) formation of a nappe pile resembling the Helvetic nappes resting above the Infrahelvetic complex. Furthermore, the finite strain pattern, temperature distribution and timing of the 3-D model is in broad agreement with data from the Helvetic nappe system. Our model, hence, provides a 3-D reconstruction of the first-order tectonic evolution of the Helvetic nappe system. Moreover, we do not apply any strain softening mechanisms. Strain localization, folding and nappe transport are controlled by initial geometrical and mechanical heterogeneities showing the fundamental importance of tectonic inheritance on fold-and-thrust belt evolution.

1 Introduction

Fold-and-thrust belts are common in nature and typically associated with orogenic belts, such as the Himalayas or the European Alps. The structural interpretation of fold-and-thrust belts is based on the interaction between the crystalline basement and the overlying sedimentary cover. Two endmember deformation styles are commonly distinguished: thin-skinned deformation, without significant basement involvement, and thick-skinned deformation, with significant basement involvement (Rodgers, 1949; Pfiffner, 2006). Due to their importance for the fundamental understanding of mountain building processes and for natural resources exploration, the formation of fold-and-thrust belts has been studied since several decades with field and modelling studies (e.g. Davis et al., 1983; Dahlen, 1984; Dahlen and Suppe, 1988; Beutner, 1977; Price and McClay, 1981; Gillcrist et al., 1987; Butler, 1989; Ramsay, 1989; Buchanan and Buchanan, 1995; Dunn et al., 1995; Mitra, 1997; Lacombe and Mouthereau, 2002; Wissing and Pfiffner, 2003; Simpson, 2011; Yamato et al., 2011; Ruh et al., 2012; Fernandez and Kaus, 2014; Bellahsen et al., 2012; Bauville and Schmalholz, 2015; Lacombe and Bellahsen, 2016; Bauville and Schmalholz, 2017). However, the mechanical deformation processes controlling fold-and-thrust belt evolution are still incompletely understood. One challenge for understanding fold-and-thrust belt evolution is that the formation, spacing, orientation and time-sequence of thrusts, shear zones and folds are controlled by two different factors. First, the mechanical deformation behaviour of rocks, which can be dominated either by brittle frictional sliding or by ductile creep and can further be strongly affected by various mechanical softening mechanisms, such as frictional strain softening, reduction of effective friction by fluid overpressure, grain size reduction or thermal softening. Second, the geometrical configuration, such as half-graben structures or orientation of sedimentary layers, and variations in rock strength, e.g. between basement and cover or within the cover by alternation of strong, such as carbonates, and weak, such as shales, sediments. Therefore, a main challenge for understanding the thermo-mechanical evolution of fold-and-thrust belts is to determine whether the major thrusts and shear zones have been controlled by a particular rheological softening mechanism or by pre-existing geometrical and mechanical heterogeneities, referred to here as tectonic inheritance. In nature, there is most likely a continuous transition between these two controlling factors.

Many studies employing numerical models have been performed with a focus on the impact of different rheological models. Studies investigated the impact of brittle, brittle–ductile, visco-plastic and visco-elasto-plastic rheological models on fold-and-thrust belt evolution, and they also studied, for example, the impact of fluid pressure and associated reduction of the effective friction coefficient (e.g. Simpson, 2011; Ruh et al., 2014; Poulet et al., 2014; Bauville and Schmalholz, 2015; Granado and Ruh, 2019). Other numerical studies focused more on the impact of tectonic inheritance, in the form of mechanical heterogeneities, on fold-and-thrust belt evolution (Wissing and Pfiffner, 2003; Bauville and Schmalholz, 2017). To illustrate the two major control factors, namely, rheology and tectonic inheritance, in a simple way, let us consider the deformation of a linear viscous material under homogeneous pure shear. Adding a circular inclusion with a smaller viscosity to the viscous material will not generate a significant shear zone inside the linear viscous material for this deformation configuration. The only possibility to generate a shear zone in the viscous material is to add a softening mechanism, such as thermal softening (Jaquet et al., 2015; Kiss et al., 2019) or grain size reduction with grain-size-sensitive flow laws (Bercovici and Ricard, 2003; Austin et al., 2008). In contrast, if the linear viscous material is sheared over a non-planar interface, resembling a half-graben, then a shear zone can develop inside the linear viscous material even without any softening mechanism (Bauville and Schmalholz, 2017). Hence, strain localization can be caused by material softening and by mechanical heterogeneity.

The majority of numerical modelling studies use two-dimensional (2-D) models, which are often able to produce results that are to a first-order approximation comparable with geological reconstructions and cross sections. However, it is well known that the style of fold-and-thrust belts can vary considerably along the strike of the belt (e.g. Hamilton, 1988; Mitra and Fisher, 1992; Mitra, 1997; Mouthereau et al., 2002; Fitz Diaz et al., 2011; Nemčok et al., 2013). Inherited, laterally varying pre-existing structures are important, for example, in the Alps where pre-Alpine laterally varying passive margin structures presumably exert a strong control on the deformation style (Pfiffner, 1993, 2011). A recent study of the Iberian passive margin by Lymer et al. (2019) highlights the complex 3-D architecture of such margins. Their results imply that fault systems disappear laterally or link together in lateral direction along the margin, consequently creating discontinuities and geometrical asymmetries. Therefore, it is important to consider the 3-D inherited heterogeneities of passive margins when studying fold-and-thrust belts that resulted from the deformation of passive margins, as is the case for the Helvetic fault-and-thrust belt (Pfiffner, 2011).

Here, we apply a 3-D thermo-mechanical numerical model to investigate the fundamental impacts of mechanical heterogeneities, representing graben structures and sedimentary layering, on the deformation style during fold-and-thrust belt formation. A particular aim is to apply our model results to the Helvetic nappe system in the Central Alps of France and Switzerland (see next section). We employ an initial model configuration that mimics a simplified upper crustal region of a rifted margin and is composed of a basement and several sedimentary units. The passive margin contains a half-graben system that has a variable depth and width in the lateral direction. Moreover, we apply laboratory-derived temperature-dependent dislocation creep flow laws for all our model units and consider a brittle frictional Drucker–Prager yield stress. We use a typical velocity boundary condition to simulate the large-scale deformation conditions during tectonic wedge formation (Simpson, 2011; Ruh et al., 2014). In order to focus on the large-scale processes, we keep our model relatively simple and concentrate on the thermo-mechanical processes on the macroscale and the impact on kilometre-scale structures. Hence, we do not consider microscale processes such as grain size reduction involving secondary mineral phases and pinning. Furthermore, our model does not include hydro-chemical coupling. Hence, we do not model processes such as fluid release and decarbonatization (Poulet et al., 2014).

The aims of our study are to (1) understand the impact of lateral variations in half-graben geometry on the deformation style, (2) investigate the importance of the spatial distribution of a competent layer and half-graben geometry during fold-and-thrust belt formation, and (3) discuss the application of our numerical models to the formation of the Helvetic nappe system.

2 Overview of the Helvetic nappe system

The Helvetic nappe system is a fold-and-thrust belt complex that was formed during the Alpine orogeny when the European passive margin collided with the Adriatic margin (e.g. Trümpy et al., 1980; Pfiffner, 2015). The system consists of a pile of tectonic nappes, which mainly comprise Mesozoic and Cenozoic sediments derived from the former European continental margin (Masson et al., 1980; Ramsay, 1981; Pfiffner, 1993; Escher et al., 1993; Steck, 1999; Pfiffner et al., 2011). Commonly, the nappe system is subdivided into the structurally upper Helvetic nappes, considered mainly as allochthonous thrust nappes, and the structurally lower Infrahelvetic complex, which can involve par-autochthonous fold nappes (Pfiffner, 1993; Pfiffner et al., 2011) (Fig. 1).

Figure 1 Simplified geological cross sections along strike of the Helvetic nappe system. (a) Glarus nappe complex (modified after Pfiffner, 2015). (b) Doldenhorn nappe (modified after Kirschner et al., 1999). (c) Morcles nappe (modified after Escher et al., 1993). (d) Morcles nappe at the Belledonne massif (modified after Epard, 1990). UH/NH represents Ultrahelvetics/North Helvetics. CS represents cover sediments.

The onset of Alpine burial of the proto-nappe system is constrained by the last deposited sediments with ages of ca. 28 to 34 Ma (Kirschner et al., 1995; Nibourel et al., 2018). Peak metamorphic conditions in the Helvetic nappe system occurred between ca. 25 and 17 Ma, indicating the end of the main phase of nappe stacking (Kirschner et al., 1995; Nibourel et al., 2018). The main phase of nappe formation and stacking occurred presumably during a period of ca. 10 to 15 Ma (Masson et al., 1980; Milnes and Pfiffner, 1980; Burkhard, 1988; Kirschner et al., 1995; Nibourel et al., 2018). Uplift and exhumation of the Helvetic nappe system occurred between ca. 20 Ma and today (Kirschner et al., 1995; Nibourel et al., 2018). We focus here on the main phase of nappe formation and stacking, and we do not consider the subsequent uplift and exhumation of the Helvetic nappe system.

The Helvetic nappe system exhibits a wide range of nappe geometries, including two commonly considered endmember nappe styles, namely, fold nappes and thrust nappes or thrust sheets (Termier, 1906; Epard and Escher, 1996). Fold nappes are recumbent folds with fold amplitudes of several kilometres and with a stratigraphic inversion in a prominent overturned limb. Thrust nappes are coherent allochthonous rock sheets that are displaced along a basal shear zone and lack a prominent overturned limb.

We consider here four simplified geological sections across the Helvetic nappe system in the Central Alps of France and Switzerland, and we focus on several prominent nappes within these sections (Fig. 1). The first-order tectonic features of these four cross sections and the associated nappes will be compared with our 3-D modelling results.

The first section includes the Glarus nappe complex of eastern Switzerland (Fig. 1a). Geological reconstructions suggest a displacement of approximately 50 km from its original location along a thin basal thrust zone that is composed of Mesozoic sediments (e.g. Pfiffner, 2015). Studies suggest that the Glarus basal thrust originates inside Carboniferous strata, allowing for the transport of the Glarus nappe consisting of Permian Verrucano units at its base (e.g. Schmid, 1975; Pfiffner, 1993, 2015). Observations on the thrust zone suggest earlier viscous-dominated deformation followed by dominantly brittle deformation (Herwegh et al., 2008). A number of studies investigated the complex deformation behaviour of the thrust zone and suggest the involvement of pressurized fluids that resulted in hydro-fracture networks and the reduction of friction at the base (e.g. Burkhard et al., 1992; Badertscher and Burkhard, 2000; Badertscher et al., 2002; Herwegh et al., 2008; Hürzeler and Abart, 2008). Poulet et al. (2014) suggested a superposition of viscous and brittle deformation mechanisms due to ductile shear heating resulting in decarbonatization and the release of overpressurized fluids causing brittle fracturing.

The second section includes the Doldenhorn fold nappe, belonging to the Infrahelvetic complex, which has been overthrusted by the Wildhorn nappe, belonging to the Helvetic nappes (Fig. 1b). The Doldenhorn nappe consists of Mesozoic and Cenozoic par-autochthonous sediments that have been squeezed and sheared out of a pre-Alpine half-graben, referred to here as the North Helvetic basin (Fig. 2). The term North Helvetic basin refers here to a Jurassic extensional basin restricted to the western region of the Helvetic nappe system and should not be confused with the Tertiary North Helvetic Flysch basin that is part of the North Alpine foreland basin. The Doldenhorn nappe is derived from the inversion of the North Helvetic basin, which is restricted to the area west and south of the Aar massif. The Gellihorn and Jägerchrütz nappes are minor nappes and their sediments are considered deposits on a basement high, likely a horst, which separated the half-graben including the Doldenhorn sediments from the more distal marginal basin, referred to here as Helvetic basin, on which the Wildhorn sediments have been deposited (Masson et al., 1980) (Fig. 2). The Doldenhorn nappe roots in the Aar basement massif. Studies indicate metamorphic peak temperatures in the Doldenhorn nappe of up to 380 ^∘C (Herwegh and Pfiffner, 2005; Ebert et al., 2007a; Girault et al., 2020). These temperatures allowed for ductile deformation and folding of the Doldenhorn nappe during nappe formation. Colder temperatures around ca. 250 ^∘C in the structurally higher Wildhorn nappe were likely responsible for a deformation style resembling more a thrust nappe.

Figure 2 Simplified palaeogeographic map of the Lower Cretaceous showing the assumed distribution of the basement massifs and sedimentary units forming the discussed tectonic nappes of the Helvetic nappe system (modified after Epard, 1990, and Pfiffner, 2015). The cross section on the top left shows a representative section across the Helvetic nappe system in the SW region of the nappe system (modified after Jaquet et al., 2018).

The third section includes the Morcles fold nappe (Fig. 1c), belonging to the Infrahelvetic complex. It is overlain by a major thrust nappe, which is termed in this region the Wildhorn super-nappe (Steck, 1999). The term super-nappe is used, because the nappe can be subdivided, from bottom to top, into the Diablerets, Mont Gond and Sublage nappes. Similarly to the Doldenhorn nappe, the Morcles nappe is considered to be the result of the closure of the North Helvetic basin and the subsequent extrusion of sediments during compressional Alpine tectonics (e.g. Ramsay, 1981; Pfiffner, 1993). This North Helvetic basin comprised kilometre-thick sequences of shale-rich units and competent carbonate units. The strongly thinned and overturned limb is in contact with the crystalline Aiguilles Rouges massif below and its autochthonous sediments. Between the Morcles and Wildhorn nappe is a minor sedimentary nappe, the Ardon nappe, which is considered to be originating from the horst region, from a similar palaeogeographic position as the Jägerchrütz nappe (Fig. 2). Estimates of metamorphic peak temperatures range between 250 and 380 ^∘C (Leloup et al., 2005; Boutoux et al., 2016; Girault et al., 2020) and therefore support a dominantly ductile deformation regime. Furthermore, the deformation of the Morcles nappe is constrained by finite strain measurements. The data highlight a pattern of increasing strain from the front of the nappe towards the root zone and also from the top to the bottom. Strain ellipses show X∕Y ratios > 400 at the contact between the overturned limb and the basement-cover (Ramsay and Huber, 1987; Casey and Dietrich, 1997). Microstructural observations of the basal mylonitic shear zone in the overturned limb of the fold nappe indicate ductile creep in the calcite-rich lithologies (Austin et al., 2008). The Wildhorn super-nappe, as a whole, resembles more a thrust sheet but exhibits significant internal deformation. For example, the Diablerets nappe is separated from the Mont Gond nappe by an isoclinal fold indicating significant ductile deformation inside the super-nappe.

Due to the structural Rawil depression, there is no continuous outcrop from the Doldenhorn towards the Morcles nappe. However, several geological reconstructions suggest, that the Doldenhorn and Morcles nappes originate from the same laterally continuous North Helvetic basin and that the Wildhorn nappe and super-nappe originate from the same laterally continuous Helvetic basin of the Mesozoic passive European margin (Burkhard, 1988; Epard, 1990). The North Helvetic basin is considered absent in the eastern region of the Helvetic nappe system, which explains the absence of significant recumbent fold nappes in the Infrahelvetic complex below the Glarus thrust. To a first-order approximation, the Doldenhorn nappe is structurally similar to the Morcles nappe in that both nappes exhibit long inverted limbs, which are absent in the Glarus nappe complex.

The fourth section (Fig. 1d) includes the Morcles nappe and the Chaînes Subalpines and is located in the French Alps. There are different interpretations concerning the tectonic relationship between the Morcles nappe and the Chaînes Subalpines in this region (Epard, 1990; Pfiffner et al., 2011). Here, we follow the interpretations of Epard (1990) and Collet (1943), who consider this region of the Chaînes Subalpines as belonging to the same tectonic unit as the Morcles nappe. In other words, we assume that in the considered region the Chaînes Subalpines and Morcles nappe derived from the North Helvetic basin (Fig. 2). In this section, no Helvetic nappes are outcropping. In contrast to the Morcles fold nappe in the third section, the nappe is here not a fold nappe but rather a thrust nappe (compare Fig. 1c and d, respectively), because there is no prominent overturned limb.

Some studies suggested mechanical explanations for the variation in nappe style within the Helvetic nappe system. The formation of fold nappes in the southwestern part of Switzerland and the lack of such fold nappes in the northeastern part are explained by lateral variations of the mechanical stratigraphy, which is the alteration of mechanically strong, such as carbonates, and weak, such as shales, sedimentary units (Pfiffner, 1993, 2011). Different thickness ratios, n, of weak to strong sedimentary units cause a different mechanical response during shortening. Low ratios, n<0.5, favour imbricate thrusting and harmonic folding, while higher ratios of n favour fold nappes and detachment folding (Pfiffner, 1993). This impact of thickness ratios on deformation style, which was derived by field observations, is also supported by 2-D numerical simulations (Jaquet et al., 2014). Moreover, Von Tscharner et al. (2016) quantified, with 3-D numerical models of viscous deformation, the impact of laterally varying half-graben depth on the folding of sedimentary layers in the half-graben. Their 3-D model results also confirm that laterally varying sediment thickness has a strong impact on fold amplification and nappe formation. However, their models do not generate thrust nappes nor the stacking of nappes, as observed in the Helvetic nappe system. A comparison of observed finite strain gradients across the Morcles fold nappe with finite strain gradients, resulting from a theoretical thermo-mechanical shear zone model utilizing calcite flow laws, suggests that the Morcles fold nappe was generated by heterogeneous shearing during Alpine shortening (Bauville and Schmalholz, 2013), as was already suggested by kinematic models (Ramsay et al., 1983; Dietrich and Casey, 1989; Casey and Dietrich, 1997).

3 Methods

3.1 Numerical method

We apply the concept of continuum mechanics to describe the deformation of rocks with a system of partial differential equations (e.g. Mase, 1970). To solve the resulting system of equations numerically, we apply the 3-D thermo-mechanical parallel code LaMEM (Kaus et al., 2016; https://bitbucket.org/bkaus/lamem, last access: June 2020) for our simulations. The equations describing the conservation of mass, linear momentum and energy are the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial v_i}{\partial x_i}}=\mathrm{0},\text{(2)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial P}{\partial x_i}}+{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}=\mathit{\rho }{g}_i,\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\rho }{C}_{\mathrm{p}}{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\lambda }{\displaystyle \frac{\partial T}{\partial x_i}}\right)+H_{\mathrm{R}}+H_{\mathrm{S}},\end{array}$$

where x_i (i=1, 2, 3) refers to Cartesian coordinates in the three spatial directions (i=1 indicates x direction, i=2 y direction and i=3 z direction), v_i represents the components of the velocity vector, P is pressure (negative mean stress), ${\mathit{\tau }}_{ij}={\mathit{\sigma }}_{ij}+P{\mathit{\delta }}_{ij}$ represents components of the deviatoric Cauchy stress tensor (with δ_ij being the Kronecker delta), ρ is density, g_i=[0 0 g] is the gravity acceleration vector with g being the gravitational acceleration, C_p is the specific heat, T the temperature and λ the thermal conductivity. The source term H_R refers to the radiogenic heat production and $H_{\mathrm{S}}={\mathit{\tau }}_{ij}{\dot{\mathit{ϵ}}}_{ij}$ for shear heating. The components of the deviatoric strain rate tensor ${\dot{\mathit{ϵ}}}_{ij}$ are defined by the visco-plastic constitutive equations:

$$\begin{array}{}\text{(4)}& {\dot{\mathit{ϵ}}}_{ij}={\dot{\mathit{ϵ}}}_{ij}^{\mathrm{vs}}+{\dot{\mathit{ϵ}}}_{ij}^{\mathrm{pl}}={\displaystyle \frac{{\mathit{\tau }}_{ij}}{\mathrm{2}{\mathit{\eta }}_{\mathrm{eff}}}}+{\dot{\mathit{ϵ}}}_{\mathrm{II}}^{\mathrm{pl}}{\displaystyle \frac{{\mathit{\tau }}_{ij}}{{\mathit{\tau }}_{\mathrm{II}}}},\end{array}$$

where ${\dot{\mathit{ϵ}}}_{ij}^{\mathrm{vs}}$ is the viscous strain rate tensor, ${\dot{\mathit{ϵ}}}_{ij}^{\mathrm{pl}}$ is the plastic strain rate tensor, η_eff is the effective viscosity, τ_ij represents components of the deviatoric stress tensor and ${\mathit{\tau }}_{\mathrm{II}}={\left(\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\tau }}_{ij}{\mathit{\tau }}_{ij}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}$ is the square root of the second invariant of the deviatoric stress tensor.

The temperature-dependent viscosity, η_eff, for the considered dislocation creep is

$$\begin{array}{}\text{(5)}& {\mathit{\eta }}_{\mathrm{eff}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(B_n\right)}^{-\frac{\mathrm{1}}{n}}{\left({\dot{\mathit{ϵ}}}_{\mathrm{II}}\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}}\exp \left({\displaystyle \frac{E_n}{nRT}}\right),\end{array}$$

where n is the stress exponent, ${\dot{\mathit{ϵ}}}_{\mathrm{II}}={\left(\frac{\mathrm{1}}{\mathrm{2}}\dot{{\mathit{ϵ}}_{ij}}{\dot{\mathit{ϵ}}}_{ij}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}$ is the square root of the second invariant of the strain rate tensor, B_n the creep constant, E_n the activation energy and R the universal gas constant. The components of the plastic strain rate tensor ${\dot{\mathit{ϵ}}}_{ij}^{\mathrm{pl}}$ are determined by enforcing the Drucker–Prager yield criterion given by

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\mathrm{II}}\le {\mathit{\tau }}_{\mathrm{Y}}=\sin \left(\mathit{\varphi }\right)P+\cos \left(\mathit{\varphi }\right)C.\end{array}$$

Here, τ_Y is the yield stress, ϕ the friction angle and C the cohesion.

We do not model erosion and sedimentation in our model, mainly because the main phase of basement uplift and exhumation started during the end of the main nappe-forming phase, after ca. 20 Ma, and because our models do not generate significant lateral topographic variations (see Sect. 4).

The system of equations is discretized with a staggered grid finite difference approach and solved with LaMEM. Material properties are advected by employing a marker-and-cell method. In order to maintain computational stability for large time steps, we employ a stabilized free surface boundary condition. For detailed information on the LaMEM code, see Kaus et al. (2016) and the documentation on the official website (https://bitbucket.org/bkaus/lamem, last access: June 2020).

3.2 Model configuration

The model configuration (Fig. 3a) has a lateral width of X=70 km, a length of Y=210 km and a height of Z=40 km. Our numerical resolution for most simulations is $\mathrm{128}\times \mathrm{256}\times \mathrm{128}$ ($X\times Y\times Z$) grid points with a mesh refinement in the Z direction using 96 grid points between Z=8 and 32 km and 16 grid points each for Z<8 and Z>32 km. Our model consists of five units which are, from top to bottom, a sticky air, a cover, a strong layer, a weak unit and a basement unit. Each unit has distinctive mechanical properties, such as, for example, a flow law, friction angle or cohesion that corresponds to natural materials. In addition to the standard parameters of the flow laws, we add a prefactor f to the creep constant B_n to facilitate a brittle–ductile transition zone at a depth of 8–10 km. The details of the mechanical properties for each unit are listed in Table 1. The initial model geometry is based on a simplified and idealized cross section of the European crustal continental passive margin. The basement constitutes the bottom of the model domain and has a maximum height of 15 km. It involves a half-graben (HG) that is separated by a horst from a larger distal basin (HB) to the right (Fig. 3b). HG represents the half-graben region in which the sediments of the Morcles and Doldenhorn nappes have been deposited, HB represents the Helvetic basin in which the sediments of the Wildhorn nappe have been deposited, and the horst between HG and HB represents the domain on which the sediments of the Ardon and Gellihorn nappes have been deposited (Fig. 2). In the reference configuration HG is becoming shallower and finally disappears towards the NE direction (positive X direction; Fig. 3a). The half-graben, between 110 km $>Y>\mathrm{90}$ km, is subjected to different geometries in our different model configurations with respect to its lateral extent (X direction). The total dimensions of the graben system are fixed with a maximum depth of 7 km, a length of 25 km (Y direction) and a width of 50 km (X direction). In contrast, the geometry of the distal basin (Y<90 km) is constant and has no variations in the X direction for all model configurations. Both half-graben and basin are filled with a weak unit, mimicking shale-rich sediments, that is overlain by a strong layer with a thickness of 1.5 km, mimicking strong carbonates. Additionally, we apply a 10 km thick sediment cover, mimicking the North Helvetic Flysch and Ultrahelvetic units. We assume that these cover units have been thrust on top of the proto-Helvetic nappe system before the onset of nappe formation. Finally, we use a low-density low-viscosity sticky air layer (η=10¹⁹ Pa s, ρ=1 kg m⁻³) with a thickness of 15 km on top to simulate a free surface boundary condition. The chosen parameters for the sticky air layer follow the recommendations by Crameri et al. (2012) for our model height of 40 km and a characteristic effective viscosity around the brittle–ductile transition of ca. 2.5×10²² Pa s.

Figure 3 (a) 3-D geometry of our reference model configuration using a laterally disappearing half-graben and highlighting the distribution of our numerical phases and their material properties. Here, the colour green refers to weak materials and dark blue to strong materials. The configuration is based on the palaeogeographic reconstructions as indicated by the abbreviations on the model: AR represents Aiguilles Rouges massif, MB represents Mont Blanc massif, HB represents Helvetic basin, NH/UH represents North Helvetic Flysch and Ultrahelvetic units, and HG represents North Helvetic graben. (b) Display of our three different half-graben systems used as initial setup. (c) Cross-sectional view at X=0 km for three models utilizing half-graben G-1 in combination with three different stratigraphic distributions. (d) Sketch displaying our boundary conditions. (e) Three cross-sectional slices along the X direction to explain our model evaluation. (f) Subfigure highlighting our passive marker overlay, which is used to distinguish the different structural units, HG (green) and HB (blue), from each other. Temperature isotherms are shown in degrees Celsius (^∘C).

Table 1 Table listing the material properties of our model units for most models, where ρ is the density, A is the pre-exponential factor, f is custom prefactor, n is the stress exponent, Q is the activation energy, λ is the thermal conductivity, ϕ is the friction angle, and C is the cohesion. Some additional parameters are constant: the thermal expansion coefficient $\mathit{\alpha }=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ K⁻¹, the heat capacity C_p=1050 J K⁻¹ and the radiogenic heat production $Q_{\mathrm{r}}={\mathrm{10}}^{-\mathrm{6}}$ W m⁻³. We use creep flow laws from Schmid et al. (1977),¹ Hansen and Carter (1983)² and Kronenberg et al. (1990).³

In order to investigate the impact of the geometry of HG on the structural evolution of our model, we use two different series of simulations. In the first one, we vary the geometry of HG in the lateral direction using three different geometries (Fig. 3b), namely, G-1 to G-3. Here, we also define our reference model A.G-1. Model A.G-1 comprises a geometry of HG, which linearly shallows and narrows out in the lateral X direction. The next configuration, G-2, is HG without lateral geometrical variations, which is bounded by a 90^∘ wall to the adjacent basement. Our third configuration, G-3, is HG with a full graben that shallows out in lateral direction while the length in the Y direction is constant.

In the second series (Fig. 3c) we employ the same HG as in the reference configuration, G-1, but modify the spatial connection between the strong layer and the underlying basement. In model B.G-1 we decrease the height of the horst in such a way that the strong layer and a part of the weak unit rests on the top of the horst. Additionally, we also thicken the left part of the strong layer which is connected to the basement in a wedge-like shape (Fig. 3c). In the last model, C.G-1, we extend the strong layer to the left of the HG (Fig. 3c). For C.G-1, the layer is not connected with the basement but underlain by an additional 0.5 km thick layer of weak units (Fig. 3c).

To test the impact of the vertical strength distribution in the basement, we performed an additional simulation, D.G-1, with the same configuration as the reference model, A.G-1, but we cut off the deviatoric stresses in the basement at 40 MPa. Numerically, this is done by setting the cohesion in the basement to 40 MPa and setting the friction angle to zero. Such a yield criterion corresponds essentially to a pressure-insensitive von Mises yield strength and mimics a semi-brittle deformation or a low-temperature plasticity.

We apply free slip boundary conditions on all sides of the model except on the top where we model a free surface boundary condition with the sticky air method (Fig. 3d). In order to mimic the kinematic conditions during tectonic wedge formation, we apply a constant velocity boundary condition on the bottom face and the left XZ face of our model (Fig. 3d). We induce bulk shortening of the model by moving the left boundary and the bottom boundary in positive Y direction with a constant velocity of v_y=1 cm yr⁻¹. This boundary condition is similar to typical sandbox analogue models of accretionary wedge formation. The bulk shortening strain rate ${\dot{\mathit{ϵ}}}_{\mathrm{bg}}$ is recalculated from the velocity boundary condition for every time step. The shortening of the model in the Y direction is balanced by an elongation in the Z direction, with no changes of the width in the X direction.

Initially, we apply a linear, vertical temperature gradient of 16.6 ^∘C km⁻¹ with a fixed temperature of 20 ^∘C at the surface and 435 ^∘C at the bottom of the model.

Furthermore, we use passive marker lines and patterns to improve the visualization of 2-D cross sections of our 3-D model (Fig. 3f). Also, we slightly change the colour scheme to highlight the two major basement structures HG and HB, together with the sediments on top of the basement horst (yellow). The darker colours inside the HG and HB indicate the strong layer, while the lighter colours correspond to the weak units. The layering and patterns are passively advected with the numerical velocity field. Hence, they do not influence the material properties and deformation and are simply there for visualization purposes.

4 Results

4.1 3-D model evolution

We first provide an overview of the general model evolution of all six performed 3-D simulations. Figure 4 shows the structural evolution for three different bulk shortenings for each model configuration. All models, except model C.G-1 (Fig. 3c), show the formation of nappe-type structures in the strong layers and stacking of the strong layer from HB on top of the strong layer of HG. Model C.G-1 does not generate any nappe-like structure or overthrusting but generates detachment folds.

Figure 4 3-D model evolution for our six simulations, using a graphical threshold to highlight the deformation of the strong layer, weak units and the basement. Rows correspond to the configuration and columns to the total bulk shortening γ_b.

In the first stage of our reference model A.G-1 (Fig. 4a), the strong layer of HB is detached from the horst by the formation of a shear zone. Both the strong layer and the weak units below are thrusted on top of the horst, culminating in a horizontal displacement of about 10 km. Our model generated a structure resembling a thrust sheet or thrust nappe. The HB is progressively closed. In the region of HB, the weak sediments and the basement are thickened with ongoing bulk shortening. Initially, the basement experiences an uplift at the right boundary due to the imposed velocity discontinuity. The basement uplift increases the topography at the right boundary of the model (i.e. the backstop). The HG shows minor signs of deformation at 23 %. The strong layer is partially sheared and squeezed out in front of HG. Additionally, we observe a slight depression of the strong layer in the rear of HG. In the next shortening stage (Fig. 4b), the closure and inversion of HG is in progress. The strong layer has been squeezed out of HG. The horizontal displacement of the strong layer of HG decreases in positive X direction associated with a decreasing depth of HG. Furthermore, the thrust nappe from HB has advanced on top of HG and the horst exhibits a dome-like shape. At 47 % bulk shortening (Fig. 4c), HG is almost completely closed, resulting in the formation of a elongated cusp dipping towards the backstop. The infill of HG has been extruded and compromises a nappe that is overthrust by a nappe from HB. The front of the nappe from HB is essentially straight along the lateral X direction, indicating that the displacement of the nappe is essentially unaffected by the formation of the nappe below, whose front varies significantly along the X direction. We only observe a slightly higher topography above the nappe from HG, causing a slight tilting of the nappe from HB towards the positive X direction.

Model A.G-2 (Fig. 4b) deviates in the structural evolution from model A.G-1. At 23 % bulk shortening (Fig. 4b), the strong layer of HG is displaced more out of HG than for model A.G-1. The basement horst is slightly sheared towards the left, which is against the shortening direction. This shearing of the horst increases with progressive bulk shortening. At 47 % bulk shortening (Fig. 4b), the HG is nearly closed at its top, but the closure of the HG did not form a cusp, as for model A.G-1. In contrast, the sediments still in the HG become thicker with depth, because for the considered geometry of the HG it is more difficult to squeeze out all the sediments (Fig. 4b). The strong layer of HG has been completely detached from the horst, showing a higher slope in front of the extruded sediments. Overall, this model shows a significantly higher degree of basement involvement, with the updoming horst separating the nappe from HG from the thrust nappe from HB. Nevertheless, the front of the nappe from HB was displaced in the same uniform manner as in the reference model, A.G-1.

Model A.G-3 displays a similar evolution as model A.G-1 (Fig. 4c). The geometry of HG results in an almost straight front of the squeezed out strong layer. There is a slight curvature at the lateral boundary between the lateral end of HG and basement (Fig. 4c). The final stage (Fig. 4c) shows a nappe originating from HG that is almost homogeneous in the lateral X direction. However, even though the depth of the half-graben is constant up to X=50 km, there is a notable decrease of nappe height from X=30 to 50 km. This tilting of the nappe front towards the positive X direction can be explained by the adjacent basement, which might affect the degree of half-graben closure due to the strength contrast between the sediment and basement units. At 47 % bulk shortening (Fig. 4c), the nappe from HB did not fully thrust over the fold nappe.

Model B.G-1 (Fig. 4d) exhibits a continuous strong layer from HG to HB across the horst, but it has the same geometry of HG than reference model A.G-1. At 23 % bulk shortening (Fig. 4d), we observe a similar deformation of the strong layer at HG than in model A.G-1. However, a larger amount of weak sediments from HB is displaced across the horst against the shortening direction. At 35 % bulk shortening, the HG is less closed compared to the same bulk shortening of the reference model (Fig. 4d). Inside the initially continuous strong layer, a shear zone develops that forms a nappe of sediments from the region of the HB. At 47 % bulk shortening (Fig. 4d), two nappe structures have formed: one consisting of sediments from HG and a structurally higher one consisting of sediments from HB. The two nappes are also stacked. However, in this model the strong layer has no connection anymore to the horst.

Model C.G-1 (Fig. 4e), with a strong layer resting above the entire basement, shows detachment folding. The lateral variations in basement geometry are not significant enough to generate a shear zone in the strong layer, which eventually would form a nappe structure. The deformation of the basement and HG is similar to model B.G-1. The detachment fold with the largest amplitude originates from the region of HG (Fig. 4e). With progressive shortening, this detachment fold is displaced along the top of the basement against the direction of shortening (Fig. 4e).

In model D.G-1 (Fig. 4f), with a stress cut-off in the basement at 40 MPa, the shortening is more homogeneously distributed in the basement, so that during the initial stages of shortening (23 %; Fig. 4f) the basement at large distance from the backstop is already thickened. Consequently, the basement uplift around the backstop is significantly lower compared to all other models. Furthermore, during nappe formation some parts of the uppermost basement are also sheared off from the basement. The nappe forming from the strong layer of HG resembles a fold nappe.

4.2 2-D numerical cross sections

We discuss the deformation in the 3-D models by analysing six 2-D cross sections parallel to the shortening direction (Y direction) but at six different locations along the lateral X direction (Fig. 5). The six cross sections are located along the lateral direction from X=0 to X=50 km with 10 km spacing.

Figure 5 Graph showing lateral cross sections for two different states of bulk shortening for the geometrical evolution of model A.G-1. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position.

For the reference model A.G-1 (Fig. 5) at 23 % shortening (Fig. 5, left column), the sediments from HB have been detached from their original position and thrusted across the basement horst all along the lateral direction (everywhere in the left column of Fig. 5). The displaced sediments from HB resemble a thrust nappe. Around the backstop (right side) there is significant basement uplift. The cross sections show, from top to bottom, the shallowing and disappearance of HG. The strong layer is already sheared out of HG, whereby the horizontal displacement is larger for deeper HG. Both nappes originating from HB and HG are deforming at the same time. The temperature around the top of the horst is ca. 300 ^∘C.

The passive marker lines in the weak sediments of HB indicate ductile flow generating passive shear folds. In the basement the passive marker symbols (crosses) indicate an increase in shear strain from the top to the bottom of the basement due to decrease in basement viscosity associated with a temperature increase. At 47 % bulk shortening (Fig. 5g–l, right column), the nappe from HB is thrust completely above the nappe from HG. The HG is almost closed with the weak sediments residing now inside a cusp between the neighbouring basement units, which dips towards the backstop. In the deeper region of HG the passive markers in the weak sediments indicate a fold nappe structure. The frontal part of the nappe from HG resembles a thrust nappe, formed exclusively by the strong layer, whereas the main part of the nappe, including also the weak sediments, resembles a fold nappe (Fig. 5g–i). Between the upper nappe from HB and the lower nappe from HG are cover sediments that have been dragged between the two nappes during overthrusting of the nappe from HB (Fig. 5g–k). The nappe from HG disappears in the lateral direction with the disappearance of HG. For the cross section without HG (Fig. 5l), there is only a nappe from HB resembling a thrust nappe. Despite the significant lateral variation of the depth of HG and the associated lower nappe, the horizontal displacement of the upper thrust nappe is essentially the same along the lateral direction. Consequently, nappe formation of sediments from HG does not affect the horizontal displacement of the overthrusting nappe from HB. The horst exhibits a significant internal deformation in the regions with a deep HG, indicated by the passive marker symbols. The temperature isotherms are affected by the deformation and thickening of the model. Generally, the sediments are getting hotter during the deformation. At 23 % bulk shortening the strong layer of HG has a temperature of ca. 300 ^∘C, while the temperature of the strong layer of HB is less than 300 ^∘C. At 47 % bulk shortening the isotherms indicate that both strong layers were heated by ca. 50 ^∘C.

In addition to the cross sections showing the structural and thermal evolution of model A.G-1, we display the same cross sections but indicating the magnitude of the deviatoric stress invariant, τ_II, to quantify the state of stress (Fig. 6). The largest stress is ca. 140 MPa and occurs at the brittle–ductile transition in the cover, at a depth of ca. 5 km. At 23 % bulk shortening (Fig. 6a–f), the top of the basement exhibits τ_II values between 40 and 80 MPa. After 47 % bulk shortening, τ_II values in the basement are strongly decreased down to 10 to 20 MPa, due to the increase in basement temperature. Overall, there are no considerable lateral variations in τ_II.

Figure 6 Graph showing lateral cross sections for two different states of bulk shortening for the second invariant of the deviatoric stress τ_II of model A.G-1. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position (colour map from Crameri (2019)).

Figure 7 Graph showing lateral cross sections for two different states of bulk shortening for the geometrical evolution of model A.G-2. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position.

The cross sections for model A.G-2 indicate overall a similar structural evolution as our reference model (Fig. 7) in the sense that the sediments are sheared out of HG and HB, form nappe-like structures and are piled at 46 % bulk shortening. However, the different initial geometry of HG also generates differences in the structural evolution. At 23 % bulk shortening the nappe from HB is less displaced across the horst for a deeper HG. The reason is that, due to the different geometry of HG, the basement horst is also significantly sheared towards the left model side and represents, hence, a less stiff obstacle for a deep HG. Therefore, both the nappe from HB and the horst are sheared together towards the left model side. For a shallower HG, the horst represents a mechanically stiffer obstacle for the nappe from HB, and hence the nappe is overthrust more to the left. At 46 % bulk shortening (Fig. 7g–l), we also observe extrusion of sediments from HG, but for the deepest regions of HG a significant amount of the sediments remains trapped between the basement at depth. Around the region where HG is deepest, the horst is sheared significantly and even sheared slightly on top of the basement that was initially to the left of HG (Fig. 7g–i). The strong layer from HG is essentially disconnected from the horst and a significant part of this strong layer resembles more a thrust nappe. The temperature evolution of model A.G-2 is similar to the one of the reference model. The results indicate that the geometry of HG has a strong impact on the structural evolution of the nappes and the basement, although the first-order structural evolution is similar to the reference model.

Figure 8 Graph showing lateral cross sections for two different states of bulk shortening for the geometrical evolution of model A.G-3. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position.

Model A.G-3 (Fig. 8) shows a similar structural evolution as the reference model. The main difference can be observed in the cross section located at the lateral boundary of HG at an X position of 50 km (Fig. 8l). There, sediments of HG have been extruded laterally out of HG, which can be seen by a pocket of sediments inside the basement.

Figure 9 Graph showing lateral cross sections for two different states of bulk shortening for the geometrical evolution of model B.G-1. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position.

Figure 10 Graph showing lateral cross sections for two different states of bulk shortening for the geometrical evolution of model C.G-1. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position.

Model B.G-1 (Fig. 9) shows a different nappe evolution than model A.G-1. At 23 % bulk shortening (Fig. 9a–f) in the region where the horst is overlain by weak sediments and a strong layer (yellow), the strong layer is sheared off the horst and pushed across HG by the strong layer from HB (Fig. 9a). In the region where the horst is overlain only by the strong layer (yellow), the strong layer is continuously sheared and dragged by the sediments from HB (Fig. 9b). The degree of shearing decreases with decreasing thickness of the yellow layer (Fig. 9b–d). At 47 % bulk shortening, HG is also almost closed and nearly all sediments have been squeezed out (Fig. 9g–f). At X=0 km (Fig. 9f) the strong layer is completely detached from the basement. The layer initially resting on the horst has been sheared above the sediments from HG and is itself overthrust by the sediments from HB. This structure resembles the vertical stacking of three nappes whose sediments were originally horizontally next to each other. With decreasing thickness of the strong layer from the horst, its shearing across the sediments from HG is also decreasing (Fig. 9i–k). Despite a horizontally continuous strong layer from HG to HB, a major shear zone developed in the strong layer and caused the generation of a nappe-like structure. The shear zone development inside the strong layer is only due to the geometrical variation of the underlying basement, because the strong layer has homogeneous material properties.

Model C.G-1, with a continuous strong layer across the entire model domain (Fig. 10), shows a very different evolution compared to the reference model, because no prominent shear zones form in the strong layer, which could develop a nappe-like structure and significant overthrusting. The strong layer initially above HG develops a detachment fold. At 23 % bulk shortening (Fig. 10a–f), this detachment fold already shows a great variability in amplitude in lateral X direction of the model. During the initial stages of folding, the core of this fold was filled with half-graben sediments and progressively displaced towards the left where the basement top is horizontal (Fig. 10a–e). In the region without half-graben, only a small amount of weak sediments were available to fill the fold core so that considerable amplification was inhibited (e.g. compare Fig; 10f and a), as is the case for detachment, or décollement, folding (e.g. Epard and Groshong, 1993; Schmalholz et al., 2002; Butler et al., 2019. During bulk shortening, the detachment fold is displaced across the basement towards the left. The HG is also closed after 47 % bulk shortening (Fig. 10g–j). Similarly to the other models the weak half-graben sediments have been extruded and some of these sediments have filled the core of the detachment fold. The green marker lines in Fig. 10g indicate that the fold was continuously fed by the half-graben sediments during ongoing deformation. Hence, the amplitude of the largest individual fold decreases towards the model side without a half-graben (Fig. 10g and l). The strong layer initially from HB forms a thick nappe-like structure with a slightly overturned layer at the front. However, no overthrusting above the sediments from HG occurred.

Compared to the reference model, model D.G-1 (Fig. 11) shows only a minor basement uplift at the right boundary at 21% bulk shortening (Fig. 11a–f). Instead, the bulk shortening is distributed more homogeneously throughout the basement, resulting in significant thickening also in the basement on the left model side. Because the top of the basement is weaker, the strong layers can shear off and displace slices of the basement (e.g. Fig. 11a). Further bulk shortening results in the extrusion from sediments from HG, resulting in a structure resembling a fold nappe. This fold nappe is overthrust from the sediments from HB, resembling a thrust nappe (Fig. 11g–l). The last deformation stage is shown for 53 % bulk shortening, because 6 %–7 % more bulk shortening is required to overthrust the sediments from HB compared to the reference model. The strong layer from HG forming the fold nappe exhibits an overturned limb that is still in contact with the basement. Thus, this model does not form any structure similar to a thrust nappe for the sediments from HG. The passive marker symbols show gentle folding of the basement, also on the model side without half-graben (Fig. 11l). In contrast, the reference model (Fig. 5l) exhibits updoming of the basement without well-developed folding.

Figure 11 Graph showing lateral cross sections for two different states of bulk shortening for the geometrical evolution of model D.G-1. Cross sections are taken in 10 km steps from X=0 to X=50 km. Additionally, we display the isothermal lines in degrees Celsius (^∘C). Columns correspond to the bulk shortening state and rows to cross-section X position.

4.3 Nádai strain and Lode's ratio

We compute the Nádai strain ε_s (Nádai and Hodge, 1963) and Lode's ratio ν (Lode, 1926) for our simulations, which are two quantities to quantify 3-D finite strain. Both parameters are computed from the natural logarithm of the finite strain principal axes ϵ₁, ϵ₂ and ϵ₃. The Nádai strain ε_s is the octahedral shear strain and provides a non-dimensional value for the strain magnitude:

$$\begin{array}{}\text{(7)}& {\mathit{\epsilon }}_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{3}}}}\sqrt{{\left({\mathit{\epsilon }}_{\mathrm{1}}-{\mathit{\epsilon }}_{\mathrm{2}}\right)}^{\mathrm{2}}+{\left({\mathit{\epsilon }}_{\mathrm{2}}-{\mathit{\epsilon }}_{\mathrm{3}}\right)}^{\mathrm{2}}+{\left({\mathit{\epsilon }}_{\mathrm{3}}-{\mathit{\epsilon }}_{\mathrm{1}}\right)}^{\mathrm{2}}}.\end{array}$$

The Lode's ratio provides information on the strain symmetry and strain regime, and it can exhibit values in the interval of [−1, 1]. In particular, ν<0 indicates a prolate strain ellipsoid (constrictional strain), ν=0 a prolate–oblate strain ellipsoid (plane strain), and ν>0 an oblate strain ellipsoid (flattening strain):

$$\begin{array}{}\text{(8)}& \mathit{\nu }={\displaystyle \frac{\mathrm{2}{\mathit{\epsilon }}_{\mathrm{2}}-{\mathit{\epsilon }}_{\mathrm{1}}-{\mathit{\epsilon }}_{\mathrm{3}}}{{\mathit{\epsilon }}_{\mathrm{1}}-{\mathit{\epsilon }}_{\mathrm{3}}}}.\end{array}$$

We display both quantities for model A.G-1 (Fig. 12) for three different cross sections, which provides a representative image of the finite strain distribution also for the other models. The first section is taken at the maximum depth of HG (X=0 km), the second at half of the depth (X=25 km) and the last at the section where the half-graben disappears (X=50 km). In general, ε_s values are highest on the right side on top of the basement horst, ranging between 3.5 to 4. The extent of this shear zone, for the same intensity, changes laterally with changing depth of HG. For example, the shear zone roots much deeper into the basin on the side without HG (compare Fig. 12a–c). Moreover, ε_s outlines the fold nappe structure, indicating that the fold nappe is surrounded by highly strained material. Here, maximal values of ε_s increase laterally with decreasing depth of HG. Another lateral difference is the intensity of the shear zone located on the left side of the model in front of the extruded sediments. This shear zone also gains intensity in lateral direction with decreasing half-graben depth. Similarly to ε_s, the Lode's ratios ν (Fig. 12d–f) pinpoint the shear zone on top of the horst and outline the fold nappe.

Figure 12 Graph showing three different cross sections of model A.G-1, displaying Nádai strain (a–c) and Lode's ratio (d–f) (colour map from Crameri, 2019). Row-wise cross sections at X=0, 25 and 50 km, respectively.

Overall, ν values are close to zero, which indicates a plane strain deformation. The largest deviations from plane strain are located in the cross section of maximum depth of HG (Fig. 12d). Furthermore, values around the fold nappe in the range of 0.5 imply a flattening regime. In contrast, the shear zone at the rear of the basement horst displays negative values that are of the order of −0.2 to −0.5, indicating a constrictional deformation.

5 Discussion

5.1 Impact of lateral geometry variations and rheological layering

Our models show the detachment of sedimentary units and their subsequent horizontal transport. Depending on the model configuration, the sedimentary units resemble fold or thrust nappes. The strain localization necessary to detach the sedimentary units and to transport them horizontally without significant internal deformation occurred without any mechanical softening. In our models, the main cause of the strain localization is the initial geometrical configuration and the variations of mechanical strength between the model units. Strain localization and associated shear zone formation due to geometrical and strength variations has been shown with 2-D numerical simulations to occur even for linear viscous material (Bauville and Schmalholz, 2017). For such strain localization, no variation of effective material properties develops, and it is hence sometimes termed kinematic strain localization, which is fundamentally different to so-called dynamic softening mechanisms in which the effective material properties change, for example due to local heating, grain size reduction or fluid infiltration (Bauville and Schmalholz, 2017). Moreover, the formation of fold nappes by pushing ductile material against a rigid obstacle, without dynamic softening, was shown, for example, with laboratory deformation experiments (e.g. Bucher, 1956) and numerical simulations (e.g. Peña and Catalán, 2004). Dynamic softening mechanisms are most likely active in nature, but their intensity and the required deformation conditions are still contentious. Furthermore, dynamic softening mechanisms most likely intensify the strain localization shown in our models, but our models show that, in principle, such dynamic softening mechanism are not essential to generate nappe-like structures.

In cross-sectional view our results are similar to results by Kiss et al. (2020), who performed 2-D high-resolution numerical simulations for a model configuration close to the one of our reference model at the lateral side with the deep HG. The evolution, timing and final structures of the nappe stack are similar in both studies. The similarity of the 2-D and 3-D results, which were obtained by two different numerical algorithms, indicates that the modelled structural evolution of the nappe system is robust. Kiss et al. (2020) focused on the impact of (i) the viscosity ratio between model units, (ii) the yield stress in the basement, (iii) one versus two strong sedimentary layers and (iv) softening mechanisms on the evolution of the fold-and-thrust belt and associated nappes. They showed that all the above-mentioned points have a considerable impact on the structural evolution. In our study, we focus on the impact of the initial 3-D geometry on the 3-D spatial evolution of structures. Our simulations show, for example, the evolution of a curvature of the front of the underlying nappe (Fig. 4) and locally significant deviations from plane strain (Fig. 12), which can also be seen by the lateral indentation of the basement by sediments (Fig. 8l). These distinct 3-D effects are only due to the initial lateral variability of the geometry of HG since the applied lateral bulk deformation is zero. Ultimately, one goal of such 2-D and 3-D models is to reconstruct the tectonic evolution of fold-and-thrust belts and to assess the differences in bulk effective viscosity of major tectonic units, such as basement and sediments.

The emplacement of the sediments from HB in the style of a thrust nappe is observed in all models with exception of model C.G.-1. The thrust nappe exhibits a laterally uniform front despite the deformation and nappe formation in HG below (Fig. 4). Even in the case of model A.G-2, where the basement horst is strongly incorporated in the deformation, we do not observe large horizontal displacement gradients in the thrust nappe. Hence, deformation of such relatively small graben systems might play only a minor role on the displacement variations along strike of fold-and-thrust belts. In our models, the uniform thrust front is the result of the initially straight boundary between basement horst and HB. Consequently, we infer that the initial large-scale basin architecture plays a major role in the geometry of orogenic salients. Several studies (e.g. Thomas, 1977; Marshak et al., 1992; Boyer and Elliott, 1982) of different orogens also indicate that the sediment basin thickness is of particular importance in the expression of salients. In the context of thickness variation, Macedo and Marshak (1999) investigated the effect of variable basin geometries during bulk shortening using 3-D sandbox models. Their study implies that basin-controlled salients are strongly controlled by the basin topography, that is, variations of sedimentary thickness inside the basin. This observation is in accordance with our study, which shows that lateral changes in the lithology, such as thickness and rheology, produce different salients (Fig. 4). There are, however, additional conditions such as the indenter shape or the direction and orientation of convergence of the colliding plates that influence the overall shape of orogens (e.g. Laubscher, 1972; Ries, 1976; Mitra, 1997). Naturally, the investigation of such large-scale boundary effects would require a different model configuration than in our study.

Volume variations related to the lateral changes in the depth of HG are always expressed in the lateral variations of the thickness of the extruded nappe. The thinning out of the initial graben structure is reflected in a decrease of height and a decrease of length of the nappe. These lateral variations also affect the lateral amplitude variation of the major detachment fold of simulation C.G-1. The local geometry around the major detachment fold resembles a thin-skinned tectonic style where most deformation is concentrated in the cover sediment. However, much of the sediments filling the core of the detachment fold originate from HG, which is more than 10 km away from the detachment fold (Fig. 11g–l). Closure of HG resulted in movement of sediments from HG into the core of the detachment fold. This result shows that detachment fold initiation and progressive evolution can be controlled by the inherited basement structure. Additionally, this finding suggests that predominately thin-skinned tectonics can passively be influenced by underlying heterogeneities due to variations in the basement architecture. Nevertheless, the expression of such detachment folds in our models also depends largely on the thickness ratio between incompetent and competent units. Our findings are, therefore, in broad agreement with previous field observations and numerical studies (e.g. Pfiffner, 1993; Jaquet et al., 2014).

Figure 13 Comparison between the geological cross sections of the Morcles nappe, Doldenhorn nappe and Glarus nappe complex with sections taken from our numerical simulations. Columnwise (left to right): (a–c) geological-cross sections, (d–f) sections from model A.G-1, and (g–i) sections from model D.G-1. Black numbers in sections of (a–c) are estimated peak temperatures in degrees Celsius taken from Ebert et al. (2008). UH/NH represents Ultrahelvetics/North Helvetics. CS represents cover sediments.

5.2 Comparison with the Helvetic nappe system

Models A.G-1 and B.G-1 are able to reproduce several first-order structural features of the Helvetic nappe system. In addition, we are also able to connect the two different main kinematic deformation phases of the western Helvetics and the eastern Helvetic nappes. The formation of the main basal thrust originates in weak sedimentary units, comparable to the Cretaceous Palfris shales, of our Helvetic basin equivalent. Similarly to the Prabé deformation phase in the west and the contemporaneous Calanda phase in the east (e.g. Pfiffner, 2015), this basal thrust aids in the transport of a laterally uniform thrust nappe in our simulations. Given the simplifications of our model, we attribute the resulting major thrust front to the formation of a basal thrust of the Wildhorn super-nappe in the west and the Glarus thrust in the east. Continuous bulk shortening leads to shear localization at the contact between the strong layer and the horst in the basin. Subsequently, the layer is detached and translated with the help of the weak sediments across the basement horst and onto the half-graben, resembling a thrust sheet or nappe. We record a first noticeable deformation of the half-graben approximately 2 Ma in our simulations. This timing is in agreement with a study by Jaquet et al. (2018), which proposes a similar time interval between the onset of the Helvetic nappes and the onset of the basal thrusting at the Morcles half-graben. Vertical extrusion and squeezing out of the half-graben sediments takes place simultaneously with stacking of the major thrust nappe on top at about 6 to 7 Ma. The resulting model nappe stack shows laterally major structural differences. The overall along-strike variation is comparable to the lateral structural variations in the Helvetic nappe system, where the prominent fold nappes are only located in the southwest (Pfiffner, 1993). At the end of our simulations we record approximately 10 Ma for the complete process of nappe stack formation. This timing is in broad agreement with studies of the Morcles nappe complex (e.g. Kirschner et al., 1996; Boutoux et al., 2016) that suggest an emplacement duration of 10 to 15 Ma.

Figure 13 shows a comparison between geological cross sections taken along strike of the Helvetic nappe systems with selected enlarged numerical cross sections in lateral direction of our model. The first cross-section (a) of model A.G-1 (Fig. 13d–f) shows our equivalent of the Morcles nappe at the maximum depth of the half-graben. The result here is similar to 2-D numerical studies by Bauville and Schmalholz (2015), who investigated fold nappe formation and nappe stack formation also in application to the Helvetic nappe system, respectively. To a first-order approximation, we reproduce an extruded fold nappe with a shape comparable to the Morcles nappe.

The strong layer is still connected to the adjacent basement horst and the deformed internal weak units, which is highlighted by the green passive marker lines resembles a recumbent fold. The inner part of the fold nappe roots into a steep cusp, which is analogous to the Chamonix zone between the Aiguilles Rouges and Mont Blanc massif (Fig. 13a). The Chamonix zone is characterized by subvertical to overturned sediments representing the cover of the Mont Blanc massif and locally by subvertical sediments representing the cover of the Aiguilles Rouges massif (Paréjas, 1922; Oulianoff, 1924; Ayrton, 1980; Ayrton et al., 1987; Epard, 1986; Pfiffner et al., 2011). The overturned cover sediments of the Mont Blanc massif suggest that the Mont Blanc massif did not behave as an effectively rigid basement block but was internally deformed (e.g. Epard, 1986). The subvertical cover sediments of both massifs, having opposite younging directions, support our numerical model, which shows an overall ductile closure of the half-graben. Consequently, we argue that the half-graben was not closed by dominantly brittle inversion tectonics during which the Mont Blanc massif was thrust as an effectively rigid block above an effectively rigid Aiguilles Rouges massif (e.g. Granado and Ruh, 2019).

Figure 14 Comparison between finite strain from the Morcles nappe (taken from Bastida et al., 2014 showing a recoloured version of the synthetic cross section with strain ellipses from pressure shadows after Casey and Dietrich, 1997, including green strain ellipses from Ramsay and Huber, 1987, the cross section is modified from Dietrich and Casey, 1989) and selected finite strain ellipses from model A.G-1 and model D.G-1.

Moreover, we find a significant amount of material from the overburden between the thrust nappe and the fold nappe but also below the fold nappe. The material between the two nappes resembles the Ultrahelvetics in the geological cross sections, whereas the material below the fold nappe corresponds to undeformed cover-sequences of the basement. Turning now to the Morcles equivalent of model B.G-1 (Fig. 13g), we find a significantly different structure of the extruded fold nappe. There are several major differences that are connected to the initial mechanical layering inside the half-graben. First, the fold nappe shows a far greater length-to-height aspect ratio than the natural analogue. Secondly, the strong normal limb is not in contact with the basement horst. Due to the continuous strong layer, the structure shows greater shearing at the contact between the fold and thrust nappe.

Proceeding with cross-section (e) of our numerical model A.G-1 (Fig. 13e), we find similarities with the geological cross section of the Doldenhorn nappe (Fig. 13b). The extruded nappe shows similar aspect ratio and shape like the Doldenhorn nappe with the length significantly exceeding the thickness. However, cross-section (h) of model B.G-1 (Fig. 13h) provides an even better match with the structure of the Doldenhorn nappe. The resulting fold nappe shows a comparable shape and size. Due to the initial rheological layering, this model forms a thin thrust sheet on top of the fold nappe. We interpret this thrust sheet as an analogue to the Gellihorn and Jägerchrütz nappes, which are essentially thin thrust sheets on top of the Doldenhorn nappe (Pfiffner, 2011).

Finally the numerical cross sections without half-graben (Fig. 13f and i) exhibit a similar displacement and shape for both of our models. The finding suggests that underlying nappe formation in HG and minor variations in the vertical rheological layering has a negligible effect on the horizontal displacement parallel to the shortening direction along strike of the major thrust sheet. To a first-order approximation, this structure is comparable to the Glarus nappe complex displayed in Fig. 13c.

Additionally, we also find a first-order agreement when comparing the finite strain pattern between field measurements from the Morcles nappe and our numerically calculated finite strain (Fig. 14). Generally, the aspect ratio of the strain ellipses increases towards the bottom and root zone of the nappe, indicating significant shearing. Furthermore, we observe less deformed ellipses in the top, normal limb of our Morcles equivalent with ellipses of higher aspect ratios (Fig. 14b and d) near the contact between the fold nappe and the thrust nappe. Figure 14b shows a better match for the finite strain pattern in the upper limb, whereas Fig. 14c displays a better fit for the contact zone. We also observe a subhorizontal orientation of the major finite strain axis towards the root zone of the nappe (Fig. 14c). Numerical finite strain computation is a useful tool to compare numerical models with natural observations. The comparison can efficiently be used to recognize where the numerical model requires adaptation to the real data.

However, we could not reproduce certain first-order features of the Helvetic nappe system with our models. The first shortcoming is the protrusion and shearing out of the strong layer of the half-graben. In 3-D view this process is expressed as a thin carpet in front of the fold nappe (e.g. Fig. 4c and f), resembling a nose-like feature in the 2-D cross sections. Shearing out of the layer inhibits the formation of an outer recumbent limb as observed in the Morcles nappe (Ramsay et al., 1983). Here, our results are in contrast with previous 2-D thermo-mechanical numerical modelling results by Bauville and Schmalholz (2015), who successfully simulated the formation of a recumbent fold limb during half-graben inversion. One of the main reasons for this deviation lies in the viscosity ratio (η_R) between overlying strong layer (η_L) and the half-graben infill (η_I). In the simulations of Bauville and Schmalholz (2015), both the layer and the infill have the same viscosity (η_R=1), whereas our simulations, using different creep flow laws for each units, show that the infill viscosity η_I can be up to 3 orders of magnitude smaller than the strong layer viscosity η_L (η_R=1000). Secondly, temperature estimates from several authors (e.g. Kirschner et al., 1999; Herwegh and Pfiffner, 2005; Ebert et al., 2007b) indicate lower peak temperature conditions along the different basal nappe contacts in comparison with our simulations. For example, Ebert et al. (2008) reported increasing temperatures of 270 to 390 ^∘C from the frontal part of the Morcles nappe to its root. Here we find a close match with root zone temperatures, whereas our front part of the fold nappe exhibits temperatures in the range of 370 ^∘C (e.g. Fig. 13d). Similar temperature comparisons between our Doldenhorn and Glarus nappe equivalents show the same trend (see Fig. 13). Moreover, the cross sections of model B.G-1 exhibit slightly higher temperatures inside the core of the fold nappe, as a consequence of the significantly thicker thrust nappe above. In summary our models show a good match for the root zone temperatures, but higher temperatures at the front of the nappe complexes compared to natural observations. Ebert et al. (2008) suggested a horizontal temperature gradient of approximately 6 ^∘C km⁻¹ along the nappe interface for all three cross sections. We only register a horizontal gradient of 2 ^∘C km⁻¹ along the basal thrust surfaces. Additional constraints are given by Leloup et al. (2005) and Boutoux et al. (2016), who indicate peak metamorphic temperatures of 320 ^∘C for the Aiguilles Rouges massif and 400 ^∘C for the Mont Blanc massif. Despite these discrepancies, our model conforms to the temperature trends of natural observations, showing increasing temperatures from the front to the root zones of the nappes. Further improvements could be made by adjusting the initial geothermal gradient or by modifying the initial geometrical configuration. For the latter case, primarily, Nibourel et al. (2018) demonstrated that the Aar massif experienced a 10–15^∘ southwards dip in relation to the isotherms. Essentially, such a configuration would, for example, place the half-graben system further up in the isotherms. This adjustment would lead to cooler temperatures in the front of the nappe and greater horizontal temperature gradients during and after nappe emplacement.

In terms of large-scale structural components our simulations lack the formation of parasitic folds, smaller imbricate thrusts or the detachment between different levels inside the major thrust complex. Modelling of such smaller-scale features would require (i) a higher numerical resolution in combination with (ii) a drastically more complex rheologically layering. In addition, field observations indicate the existence of several large-scale shear zones in the Mont Blanc and Aar massif. Our models do not reproduce such features, mainly because of our mechanically homogeneous basement unit which represents a large simplification compared to the natural complexities. There are, for example, several studies regarding the Aar massif that correlate the formation of ductile and brittle shear zones to inherited pre-Alpine heterogeneities (e.g. Berger et al., 2017; Mair et al., 2018). In particular, structures such as foliations, mafic dykes or folds present mechanical anisotropies that can culminate in localized strain and shear zones (e.g. Bell, 1978; Carreras et al., 2010; Herwegh et al., 2017; Wehrens et al., 2017).

Finally, our model does not fully resolve the large-scale kinematics and intricate deformational history of the Helvetic nappe system. For example, reconstructions between the SW and the NE of the Helvetics show a differential horizontal displacement of the thrust nappes of up to 50 km (e.g. Pfiffner, 2015). Here, variations are likely due to the initial basement geometry and convergence of subducting plate, e.g. obliquity of the plate boundary which has not been taken into account in our model.

Figure 15 Comparison between selected numerical cross sections and simplified geological cross sections. Row-wise (top to bottom) (a) Simplified Morcles section modified after Pfiffner (2015) (indicated in the map of Fig. 1 with the dashed profile line) and (d) section from model D.G-1. (b) Simplified Morcles section modified after Escher et al. (1993) and (e) cross-section model F.G-1. (c) Simplified section modified after Epard (1990) and (f) section from model A.G-2. UH/NH = Ultrahelvetics/North Helvetics. CS = Cover sediments.

5.3 Comparing geological with modelled cross sections: the Morcles nappe

We show and compare three different, simplified geological cross sections of the Morcles nappe with our numerical results (Fig. 15). The cross sections are from different locations along strike of the Morcles nappe and highlight, together with our numerical cross sections, the lateral variability of the Morcles nappe. Differences in the presented geological cross sections most likely do not exclusively indicate true lateral variations in the tectonic structures, but they may also result from different geometrical projection methods and different geological interpretations. Such potential apparent differences should be kept in mind when comparing numerical models with constructed geological cross sections, which can be considered geometrical models that are also based on assumptions like the numerical models. We start with a comparison of a geological cross section of the Morcles nappe near the Sanetsch Pass and a section of model B.G-1 (Fig. 15a and d). Here, the geological cross section shows a relatively thin Morcles nappe with a sheared lower and upper limb. This observation is similar to our cross section which exhibits an elongated sheared fold nappe with sheared lower and upper limbs. Moreover, in this model we observe a significant deformation of the basement at the contact of the cover sediments, which is in agreement with microstructural observations that suggest a brittle–ductile emplacement of the Morcles nappe (Ebert et al., 2007a; Austin et al., 2008). Figure 15b shows a Morcles cross section by Escher et al. (1993) from further southwest in the nappe system. In this case, model D.G-1 (Fig. 14e) provides a better first-order fit than B.G-1. Our fold nappe exhibits a strongly sheared recumbent limb that still reaches into the root zone of the nappe. Furthermore, the strong layer is still connected to the basement horst and is in contact with the thrust nappe where it is sheared, displaying a similar saw-tooth shape as observed in the geological cross section. Also, the overall aspect ratio is closer to the geological cross section than in model B.G-1. In contrast, both of the previous numerical sections do not match the observations for the nappe structure even further southwest in France near Megève (Fig. 15c). We assume that in this region the Morcles nappe and the Chaînes Subalpines constitute one large tectonic nappe that originated from the North Helvetic basin. Here, the sediments are rather squeezed out of the half-graben due to a greater thickness of shales instead of forming a recumbent limb (Epard, 1990). Therefore, we suggest that a cross section of model A.G-2C (Fig. 15f) shows a closer match to the geological structure in this region. In summary, different cross sections of the Morcles nappe can be compared to different numerical models with different initial conditions regarding the geometry or rheology. Hence, it is a challenge to capture the evolution of 3-D geological structures with a single 3-D model, because of the geometrical and rheological uncertainties and the likely significant natural variability. Furthermore, our results imply that one has to be even more cautious in the application of a single 2-D model to the formation of 3-D fold-and-thrust belts and geological structures in general.

6 Conclusions

The presented 3-D numerical simulations show the formation of a fold-and-thrust belt resulting from the wedge-type deformation of the upper crustal region of a passive margin. The formation of the sedimentary nappes results from the shearing off and detachment of sedimentary units form the half-graben, horst and basin. Nappe detachment, transport and stacking occur for standard visco-plastic rheological models without any applied material, or dynamic, softening mechanisms. Nappe formation mechanisms and their geometry are controlled by the initial basement geometry and the strength contrast between basement and cover sediments and the strength contrast within the sediments. The results indicate the fundamental importance of tectonically inherited structures on the evolution of fold-and-thrust belts and tectonic nappes. Consequently, the results emphasize the importance of the initial geological situation before fold-and-thrust belt formation.

Modelled nappe-like structures generated from sediments from the half-graben with laterally varying depth show that the nappe geometry strongly depends on the amount of sediments available for nappe formation. The calculation of 3-D finite strain shows that the deformation during the formation of nappes originating from the half-graben largely deviates from a 2-D plane strain deformation and is of a flattening type (Lode's ratio > 0). Nappes originating from the laterally homogeneous basin show a more or less laterally straight displacement front, indicating that the laterally heterogeneous deformation during nappe formation in the half-graben region does not considerably affect the overthrusting nappe.

The applied strength of the basement has a strong control on the resulting nappe geometry. A relatively stronger basement favours the formation of nappes resembling a thrust nappe, while reduced basement strength (modelled here by a stress cut-off at 40 MPa) favours the formation of fold nappes.

Nappe-like structures and a pile of nappes form only in the models when the strong sedimentary layer is limited laterally by a basement block, i.e. when the layer is buttressed by the basement. If the layer is horizontally continuous above the basement, then detachment folding is favoured and a nappe pile does not form.

The initial model configuration is based on geological reconstructions of the Helvetic nappe system. The modelled 3-D structures represent a thermo-mechanical model for the tectonic evolution of the Helvetic nappe system, which can explain several first-order observations of this system:

1. Formation of a nappe resulting from the closure of a half-graben is observed. Depending on the model configuration, this nappe can be more similar to a fold or a thrust nappe. The modelled nappes are applicable to the Morcles and Doldenhorn nappes of the Helvetic nappe system.

2. Formation of a laterally homogeneous thrust nappe is observed, which overthrusts and is stacked above the underlying nappe resulting from half-graben closure. This nappe is applicable to the Wildhorn and Glarus nappes.

3. The detachment of minor sediment units originally located on the horst and their emplacement between the upper thrust nappe and the underlying fold nappe from the half-graben is observed. These minor nappes are applicable to the Ardon, Jägerchrütz and Gellihorn nappes.

4. The entrapment of weak sediments is observed, which were originally situated structurally above the sediments eventually forming nappes between the two major nappes. These entrapped sediments are applicable to the Ultrahelvetic units.

5. The modelled temperatures, temperature gradients, finite strain gradients and duration of deformation are in overall agreement with data from the Helvetic nappe system.