The Role of Isostasy in the Evolution and Architecture of Fold and Thrust Belts

Fold-thrust belts (FTB) form at plate convergent boundaries, often adjacent to a cordillera or a collisional orogen or adjacent to a subduction zone as part of an accretionary wedge. FTBs develop within a dynamic framework in which tectonic processes, surface processes, and isostatic processes interact under rock mechanical properties that evolve in response to temperature, strain, and fluid-rock interactions. Their development has been well-studied through analog and numerical experiments (e.g., [1–7]). These studies show that the architecture of an FTB is influenced by many processes including convergence rate, preexisting basement structures [8], the mechanical stratigraphy of the sedimentary pile [9], and surface processes [10–12]. During the build-up of an FTB, tectonic and sedimentary loading along with the critical taper angle [13] controls the internal deformation of the wedge [14]. The increasing load further influences FTB architecture by forcing local subsidence and/or more regional plate bending, which provides accommodation space for further sedimentation, lowers the surface slope, and increases the height of the wedge. Hence, the bending and subsidence of the surface supporting the load are of critical importance.

When numerical experiments aim to study high-resolution small-scale FTB development, the model domain is generally computationally restricted to the upper crust. Therefore, modellers forego processes that occur underneath the model such as isostasy (Figure 1), that may influence the evolution of an FTB. As a result, numerical experiments are often run with a rigid model base. An alternative approach is to employ a boundary condition at the base of the model that simulates processes occurring underneath the model domain.

For a line load $qx$⁠, the full flexural response $wx$ is given by the summation of these solutions. For a given surface load, the wavelength and amplitude of the bending vary with the plate’s elastic thickness (⁠$Te$⁠). This methodology is described in detail in [18].

Though the elastic formalism presented above explains reasonably well the bending of an oceanic lithosphere entering a subduction zone, its applicability to FTBs developing adjacent to orogenic domains is debatable. Indeed, in these domains where warmer geotherms are expected, the isostatic response to loading is—at least in part—accommodated by viscous mass redistributions inside and/or below the plate (e.g., [19, 20]). Isostatic glacial rebound demonstrates the importance of viscous flow in the sublithospheric mantle. The rate of isostatic postglacial uplift, from a few mm yr^-1 to a few cm yr^-1 [21–23], is comparable to tectonic loading rates. To properly capture the magnitude and wavelength of isostatic subsidence, numerical models of FTBs must include the entire lithosphere as well as a layer of asthenosphere of significant thickness (i.e., a total depth of 200-250 km). Typically, the isostatic boundary condition attached to the base of the computational grid allows the asthenosphere to flow in-and-out to maintain the pressure at the base of the model constant. However, to properly capture in sufficient details the structural architecture of FTBs, and resolve vertical motion of few 100 s of meter, our numerical experiments require a computational grid of $spatialresolution<100m$⁠. The combination of these two requirements (model 200-250 km deep and $gridresolution<100m$⁠) is out of reach of most advanced numerical codes. Therefore, in our numerical experiments of FTBs, we implement a boundary condition directly underneath the FTB that mimic Airy-like isostasy. For this, we use a high-viscosity virtual layer at the base of the model and a boundary condition at the base of the computational grid allowing that layer to move down and out of the computational domain to maintain the basal pressure constant. The subsidence of the virtual viscous layer captures the integrated viscous flow anywhere in the lithosphere and the subjacent asthenosphere that accommodates isostasy (Figure 1). As the magnitude of the subsidence and its wavelength is a function of the density of the virtual high-viscosity layer, we can explore a range of an Airy-like isostasy by simply varying the density of the high-viscosity layer.

In what follows, we explore how the interaction between horizontal mass transfer, related to tectonic processes, and vertical mass transfer, related to either elastic flexure or viscous flow, impacts the overall architecture of FTBs. We present a series of numerical experiments where the progressive stacking of thrust sheets occurs under contrasting isostatic mechanisms. To simulate isostasy, we use two distinct basal boundary conditions to mimic the elastic bending of the lithosphere or its subsidence in response to the ductile flow. Our experiments illustrate how different modes of isostasy influence the architecture of FTBs.

Our numerical experiments solve the conservation of energy, mass, and momentum equations on a cartesian grid, using Underworld, a coupled thermomechanical code [24–26]. The governing equations and numerical method are described in detail in [24]. In our experiments, the numerical domain is 64 km long and 16 km deep, with a grid resolution of 80 m. Model materials are elastic-viscoplastic, and we disregard the role of temperature and strain rate on the viscosity as FTBs form under moderate temperatures. In the formulation of plastic rheology, cohesion and coefficient of friction evolve to account for strain-weakening. Material properties are listed in Table 1. Our rheological parameters deliver experimental outcomes reasonably similar to natural examples and consistent with other studies [8, 27–29].

From top to bottom, the model consists of “sticky air” [30], 200 m of “loose sediment” (cohesion of 0 Pa), 4 km of sedimentary rock made up of eight layers of equal thickness but contrasting competence, a 3 km thick basement, and an underlying 4 km thick basal layer (Figure 2). At the top of the model, we impose a free slip condition. Material enters from the right wall at a rate of 1 cm yr^-1 with no slip in the vertical direction. On the left wall, material below the basement is evacuated horizontally at a rate of 1 cm yr^-1. Above the basement, the left wall acts as a rigid backstop. We apply our isostasy boundary conditions (flexural or ductile flow) to the base of the model.

We run two suites of models using either flexural isostasy or ductile flow isostasy. In alignment with previous studies (e.g., [6, 16, 17]), we simulate isostatic flexure with an elastic beam formulation at the base of the model (Figure 2(b)). We use GFlex, a widely used flexural response model [18], to calculate the amount of flexure at each model timestep and use that to displace the model vertically. In line with other studies, Youngs’ modulus (⁠$7×1010Pa$⁠) and Poisson’s ratio (0.25) are held constant, and we run experiments with varying $Te$⁠.

For models where isostasy is accommodated by viscous flow, we use an Airy-like boundary condition at the base of the model. This condition maintains constant basal pressure by kinematically evacuating material from the isostatic scaling layer and through the base of the computational grid, allowing the model to subside in response to surface loading. The magnitude and wavelength of the subsidence are controlled by the density of the basal layer. The higher the density is, the smaller the ductile outward flow needed to maintain the basal pressure, and therefore, the lesser the subsidence and the longer its wavelength (Figure 2(c)). A limitation of this approach is that there is no isostatic compensation postthickening. Erosional processes are accounted for via slope instabilities and a hillslope diffusive law, with a diffusion coefficient of 1 m² yr^-1. We present six key models from two series of numerical experiments (see complete set in the supplementary section). Models LF1–LF3 use the elastic flexure isostasy condition, and DF1–DF3 use theAiry-like ductile flow condition. Experiments DF1 and DF3 deliver smaller/higher subsidence over longer/shorter wavelength, respectively. Experiments LF1 and LF3 have smaller $Te$ and larger $Te$ and therefore a larger and smaller flexural response, respectively. Experiments DF2 and LF2 are intermediate scenarios. Experiments are compared after reaching 21 km of shortening (approx. 33% shortening). Experiments are evaluated in terms of their topography, internal structure, faults development, stress, and strain rate.

This numerical experiment (Figure 3(a)) approximates the situation of a stronger, more viscous lithosphere. Applying a scaling density of 10,000 kg m^-3 to the basal layer delivers a smaller isostatic subsidence and an enhanced subsidence wavelength. In this experiment, the FTB forms above a basement, which remains subhorizontal. Faults form in sequence from the hinterland (to the left) to foreland (to the right). Individual faults are short-lived and accommodate little shortening before becoming inactive (Figures 4(a) and 5(a)). The basement is minimally deformed and is detached from the overlying sedimentary sequence, which accommodates most of the shortening. The stacking sediments into an FTB generate approximately 2 km of topography, which tapers down towards the leading fault. Towards the hinterland, the fold belt develops a gravitationally unstable high topography, which results in gravity-driven faulting and reactivation of older faults.

This experiment (Figure 3(c)) approximates the situation of a weaker, less viscous lithosphere. A scaling density of 3200 kg m^-3 to the basal layer enables a larger subsidence over a shorter wavelength. The outcome is an FTB which is narrower compared to experiment DF1, as deformation is localized to a narrower region of intense shortening and deformation. While the sedimentary stack is thickened by about 2 km, its topography is no more than a few hundred meters due to the strong subsidence of the basement. A set of near vertical faults accommodate the differential vertical motions associated with isostatic adjustment. These steeply dipping faults interfere with the set of structures that accommodate horizontal shortening. In addition, the sequence of thrusts that develop in the cover sequence consists of back thrusts and out-of-order thrusts, which complicates the FTB’s architecture (Figures 4(c) and 5(c)). This contrasts with the simpler in-sequence thrusting documented in experiment DF1 and suggests that Airy-like ductile isostasy may have a distinct structural signature.

In this experiment, the scaling density is 4000 kg m^-3 (Figure 3(b)). The magnitude and the wavelength of the subsidence fall in between the two previously discussed end-member experiments. Experiment DF2 produces up to 500 m in topography and a relatively narrow FTB, compared to experiment DF1 but wider compared to experiment DF3. The basement subsides as thrust sheets stack. Gravitational instabilities and gravity-related faulting add to the structural complexity. Subsidence is maximized underneath the load created by the overturned basement nappes.

Model LF1 (Figure 3(d)) has a $Te$ of 90 km (quarter-wavelength of 242 km). This model represents the end-member scenario where the tectonic loading exceeds the vertical accommodation space created by flexure. The thrust system generally forms in sequence. Due to the high elastic thickness, the basement retains its initial horizontal profile. This experiment produces approximately 3.5 km of topography, with minimal basement flexure. The high topography towards the backstop produces a gravitational instability and collapse, which contributes to the structural complexity of this model towards the hinterland.

Model LF3 (Figure 3(f)) has a $Te$ of 10 km (quarter-wavelength of 47 km). Here, the small elastic thickness of the lithosphere allows the basement to easily flex as the FTB load is emplaced. Faults form in sequence from hinterland to foreland with relatively little complexity. This model generates approximately 2 km of topography and undergoes a maximum subsidence of approximately 1.3 km near the left wall. The low topography of the model inhibits the formation of gravitational instabilities to the extent seen in model LF1.

This experiment (Figure 3(e)) represents an intermediate scenario between the two previous endmembers. LF2 has a $Te$ of 20 km (quarter-wavelength of 78 km), allowing the model to flex as a response to the FTB load, resulting in approximately 2.4 km of topography. The deformation style of this experiment becomes simpler later in its evolution and towards the foreland and consists of structural repeats. Large horizontal faults cut across many preexisting faults. Toward the model’s foreland, younger faults detach along these decollement zones (Figures 4(e) and 5(e)).

In the ductile flow, Airy-like, and suite of models, while all experiments develop an in-sequence thrust system verging towards the foreland, the isostatic adjustment tends to impede the outward propagation of FTBs. High strain rate regions, documenting the active deformation front, are further towards the foreland as the isostatic rate decreases (Figures 5(a)–5(c)). Consequently, when FTBs are associated to a strong subsidence driven by localized ductile flow isostasy, FTBs tend to be narrower, less elevated, and thicker. A stronger subsidence also gave rise to more complex FTB architecture as thrust faults are active longer, accommodating more shortening over a smaller region, and do not conform to a simple in-sequence kinematic (Figure 4). In these fold belts, a set of vertical faults accommodates isostatic subsidence.

Conversely, in the flexural isostasy suite of models, we found that an increase in flexural subsidence does not correspond to an increase in structural complexity. However, similarly to the Airy-isostasy models, when high topographies develop (for high $Te$⁠), gravitational instability and collapse near the hinterland of the models contribute to a higher structural complexity. A key difference between the two modes of isostasy is the relatively planar basement profile that results from the flexural isostasy condition, compared to the Airy-isostasy. This simpler basement profile likely encourages the propagation of bedding-parallel faults (i.e., décollement) that are absent in the Airy-like suite of models with a high isostatic rate. This suggests that décollement may be signature of FTBs that form under high $Te$ or low isostatic subsidence rate conditions, whereas vertical faulting may be signature of a relatively high isostatic rate due to a more ductile basement.

While 2-D numerical experiments cannot capture the diversity and range of complexity that exists in nature, insights from our numerical experiments can nevertheless tentatively be applied to natural examples. Here, we compare our results with the Subandean FTB (Figure 6). Along the South American Cordillera, the Neogene Subandean zone is the transitional region between the Cordillera’s elevated hinterland and the undeformed platform of the South American continent. Along the Subandean zone, the $Te$ varies from a maximum of ~90 km in the bend region of the Andes within Bolivia (latitude of 20°S), to a minimum of ~10 km away from the bend region towards Peru to the north and Argentina to the south (approx. latitudes 13°S and 25°S, respectively) [31, 32]. This change in $Te$ correlates with strong variations in surface heat flow [32]. Hence, this region is ideal to assess how the strength of a plate impacts the development of FTBs, via the modulation of isostasy. In a pioneer study, [31] noted that in parts of the Andes where a high $Te$ minimizes isostatic flexure, fold belts are wider and structurally simpler compared to regions where $Te$ is lesser. More recently, [32] reemphasize this point, finding a strong correlation between structural style and $Te$⁠. Additionally, ductile flow of the lower crust [33] and mantle [34] influence how shortening is accommodated across the Andes. Ouimet and Cook [35] have suggested that there is redistribution of material in the lower crust through ductile flow, related to crustal thickening and surface uplift. Postseismic uplift reveals mantle and lower crust viscosities as low as 10¹⁸ Pa s in parts of the Andes [36, 37].

The Peruvian Subandean zone overlies a relatively weak lithosphere with a $Te$ of 10 to 20 km [38, 39]. Moreover, parts of the Peruvian Andes are in isostatic equilibrium at the Moho [40]. The weak lithosphere paired with a slow convergence rate makes this area an ideal case study for the scenario in which isostatic subsidence keeps up with the tectonic loading. This region has accommodated up to 40 km of shortening and is characterized by thin-skinned imbricate thrusting [41], disharmonic folding, back thrusts, and triangle zones [42]. In this region, the basement has recorded a strong subsidence localized where the imbricate stack of thrust sheets is at its thickest (Figure 6(b)). Our numerical experiment DF3 suggests that this strong localized subsidence could be accommodated by ductile flow in the basement, rather than elastic flexure. The geometry of the basement-FTB contact contrasts with that observed further south in the Bolivian Subandean zone.

In the Bolivian Subandean zone (Figure 6(c)), $Te$ increases between 40 and 90 km [32], in a region that was shortened by approximately 40% over the last 12 m.y. [43]. The FTB architecture consists of lateral structural repetition of decreasing age in the direction of propagation [43]. Locally, structural complexities derive from the activation of multiple decollement levels [44]. Unlike the Peruvian Subandean zone, the cover is entirely detached from the basement, the surface of which is planar and gently tilted towards the hinterland. The Bolivian Subandean zone is approximately 30% wider than the Peruvian Subandean zone. These attributes are consistent with numerical experiment DF1 and LF1, which suggests that, in the Bolivian Subandean zone, ductile flow in the basement may not be important.

To test this idea, we designed two numerical models (models A and B) that approximate the thickness and stratigraphy of the Subandean zone. Both models use the Airy-style isostasy condition to test the role of ductile flow in the structuration of the Subandean zone. In model A, we maximize ductile flow and isostatic subsidence to reflect the Peruvian Subandean zone, and in model B, we minimize ductile flow and subsidence to reflect the Bolivian Subandean zone. The material and rheological setups between models A and B were held constant to isolate the effect of ductile flow on FTB evolution. The material setup and boundary condition of these models are summarized in Figure 7 and Table 2.

Model A (Figure 8(a)), in which isostasy subsidence keeps up with the tectonic loading, shares key features with its Peruvian analogue. In particular, the basement records a strong subsidence underneath a relatively narrower FTB. Shortening in the FTB is accommodated by imbricated thrust, steepening faults, pop-up and pop-down structures, and back thrust. In comparison, model B (Figure 8(b)), which reflects a less ductile Bolivian zone, displays gently sloping top basement profile, over which a wider, structurally simpler FTB evolves through the progressive lateral propagation of structural repetition towards the foreland.

Our models do not attempt to capture the details of the structural architecture of the FTBs in the Bolivian and Peruvian Subandean zones. However, the geometry of the basement-basin interfaces is reproduced to the first order and agrees with the proposition that strong subsidence underneath an FTB can be explained by ductile flow in the basement.

Our numerical experiments confirm that the structural architecture of FTBs is sensitive to the interplay between the vertical motions induced by isostasy and the horizontal motions due to tectonics. They also show that the mode of isostasy is a critical to the development of FTB. Our experiments show that when isostasy is accommodated by viscous flow, subsidence beneath the surface load is localized, impeding the outward propagation of the FTB and limiting topography. This leads to narrower and structurally more complex FTBs, in which steeply dipping faults accommodate a stronger subsidence. When viscous flow is subdued and subsidence is less important, FTBs propagate farther outwards via the lateral repetition of simpler structures. Importantly, the smaller isostatic subsidence leads to higher topography. We conclude that in FTBs associated to orogenic systems, viscous flow plays an important role in their development. For more rigid lithospheres, where isostasy is better described by lithospheric flexure, we found that a change in $Te$ and therefore a change in flexure play a lesser role in the structural style of FTBs. Nevertheless, larger $T$_e delivers FTB with higher topography, with increased erosion and sediment flux, into adjacent basins.

Data supporting the results can be found in the manuscript text, figures, and discussed references. The underlying open-source framework Underworld used to run our models can be found here doi:10.5281/zenodo.1436039. Input scripts are available at doi:10.5281/zenodo.6629741.

An earlier iteration of this manuscript was presented at a conference here doi:10.5194/egusphere-egu21-13767.

The authors declare that there is no conflict of interest regarding the publication of this paper.

This research was supported by the National Computational Infrastructure (NCI), through the National Computational Merit Allocation Scheme, the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, and Australian Research Council grants ARC-IH130200012 and ARC-LP190100146. We thank J. Giordani and Dr. R. Beucher for their expert support with Underworld and two anonymous reviewers for contributing to greatly improve the manuscript.