At what stress level is the central Indian Ocean lithosphere buckling?

Internal contrasts in strength are responsible for lithospheric buckling. These are quantified by comparing the Indian Ocean data to two-dimensional visco-elasto-plastic numerical models where the material properties depend on temperature and pressure. The central Indian Basin is known for its intraplate seismicity and long wavelength undulations of the sea floor and associated gravity signatures. To simulate the amplitudes of undulations that reach 1 km within about 11 Ma of compression and 60 km of shortening, the required mean yield strength of the lithosphere is 400 MPa. If either a hydrostatic fluid pressure is considered in the crust or a mechanical decoupling at the depth of the Moho, small crustal wavelengths are superimposed on the long wavelength deformation, in agreement with observations. It is then possible to match the alternative indication that buckling commenced only 4 Ma ago, with a total amount of shortening of 30 km, and with a required yield strength of 200 MPa. About 10% of homogeneous thickening accompanies buckling. Taking into account variable thermal conductivity demonstrates that a change in the geotherm is sufficient to increase buckling amplitudes by 10% within 5 Ma

Control of Fault Weakening on the Structural Styles of Underthrusting-Dominated Non-Cohesive Accretionary Wedges

Underthrusting is a typical process at compressive margins responsible for nappe stacking and sediment subduction. In nature, underthrusting is often associated with weak basal faults, although static mechanical analysis (critical taper theory) suggests that weak basal faults promote accretion while strong basal faults promote underthrusting. We perform mathematical analyses and numerical simulations to determine whether permanent fault weakening promotes or inhibits underthrusting. We investigate the control of permanent fault weakening on the dynamics of a strong-based ( (1 − λ ∗ b )μ b ≈ (1 − λ ∗ )μ ) non-cohesive wedge ( μ and μ b are internal and basal friction, respectively). We control the wedge material strength by a spatially constant fluid overpressure factor ( λ ∗ ), and fault strength by a plastic strain weakening factor ( χ ). First, we use the critical taper theory to determine a mechanical mode diagram that predicts structural styles. Then, we perform numerical simulations of accretionary wedge formation to establish their dynamical structural characteristics. We determine a continuum of structural styles between three end-members which correspond to the theoretical mechanical mode transitions. Style 1 is characterized by thin tectonic slices and little to no underthrusting. Style 2 shows thick slices, nappe stacking, and shallow gravity-driven tectonics. Style 3 displays the complete underthrusting of the incoming sediments, that are exhumed when they reach the backstop. We conclude that in the condition of an initially strong wedge base, permanent fault weakening promotes underthrusting. Thus, this contribution enlightens the control of the dynamic evolution of material properties on the formation of subduction channels, slope instabilities, and antiformal nappe stacks

Prediction of faulting from the theories of elasticity and plasticity: what are the limits?

Elasticity, rigid-plasticity and elasto-plasticity are the simplest constitutive models used to describe the initiation and evolution of faulting. However, in practice, the limits of their application are not always clear. In this paper, we test the behaviour of these different models using as examples tectonic problems of indentation of a die, compression with basal shear, bending of a plate and normal faulting around a dike. By comparing the results of these tests, we formulate some guidelines that may be useful for the selection of an appropriate constitutive model of faulting. The theory of elasticity reasonably predicts the initiation of the fault pattern but gives erroneous results for large strains. The theory of rigid-plasticity is more appropriate for large deformations, where the geometry of faults can be found by the method of characteristics. This method works well for zones of failure that are not severely constrained by elastic material outside e.g. when faults are connected to the free-surface, a viscous substratum or a zone of weakness. Non-associated elasto-plasticity is the most complete theory among those considered in this paper. It describes the evolution of faults from the initiation of localized deformations to the formation of a complicated fault network

3-D numerical modelling of crustal polydiapirs with volume-of-fluid methods

Gravitational instabilities exert a crucial role on the Earth dynamics and in particular on its differentiation. The Earth’s crust can be considered as a multilayered fluid with different densities and viscosities, which may become unstable in particular with variations in temperature. With the specific aim to quantify crustal scale polydiapiric instabilities, we test here two codes, JADIM and OpenFOAM, which use a volume-of-fluid (VOF) method without interface reconstruction, and compare them with the geodynamics community code ASPECT, which uses a tracking algorithm based on compositional fields. The VOF method is well-known to preserve strongly deforming interfaces. Both JADIM and OpenFOAM are first tested against documented two and three-layer Rayleigh–Taylor instability configurations in 2-D and 3-D. 2-D and 3-D results show diapiric growth rates that fit the analytical theory and are found to be slightly more accurate than those obtained with ASPECT. We subsequently compare the results from VOF simulations with previously published Rayleigh–B´enard analogue and numerical experiments.We show that the VOF method is a robust method adapted to the study of diapirism and convection in the Earth’s crust, although it is not computationally as fast as ASPECT. OpenFOAM is found to run faster than, and conserve mass as well as JADIM. Finally, we provide a preliminary application to the polydiapiric dynamics of the orogenic crust of Naxos Island (Greece) at about 16 Myr, and propose a two-stages scenario of convection and diapirism. The timing and dimensions of the modelled gravitational instabilities not only corroborate previous estimates of timing and dimensions associated to the dynamics of this hot crustal domain, but also bring preliminary insight on its rheological and tectonic contexts

Elastic flexure controls magma trajectories and explains the offset of primary volcanic activity upstream of mantle plume axis at la Réunion and Hawaii hotspot islands

Surface volcanism at la Réunion and Hawaii occurs with an offset of 150–180 km upstream to the plume axis with respect to the plate motion. This striking observation raises questions about the forcing of plume-lithosphere thermo-mechanical interactions on melt trajectories beneath these islands. Based on visco-elasto-plastic numerical models handled at kilometric resolution, we propose to explain this offset by the development of compressional stresses at the base of the lithosphere, that result from elastic plate bending above the upward load exerted by the plume head. This horizontal compression adopts a disc shape centered around the plume axis: (i) it is 20 km thick, (ii) it has a 150 km radius, (iii) it lays at the base of the elastic part of the lithosphere, i.e., around ∼50–70 km depth where the temperature varies from ∼600 °C to ∼750 °C, (iv) it lasts for 5 to 10 My in an oceanic plate of age greater than 70 My, and (vi) it is controlled by the visco-elastic relaxation time at ∼50–70 km depth. This period of time exceeds the time during which both the Somalian/East-African and Pacific plates drift over the Reunion and Hawaii plumes, respectively. This indicates that this basal compression is actually a persistent feature. It is inferred that the buoyant melts percolating in the plume head pond below this zone of compression and eventually spread laterally until the most compressive principal elastic stresses reverse to the vertical, i.e., ∼150 km away from the plume head. There, melts propagate through dikes upwards to ∼35 km depth, where the plate curvature reverses and ambient compression diminishes. This 30–35 km depth may thus host a magmatic reservoir where melts transported by dykes pond. Only after further magmatic differentiation can dykes resume their ascension up to the surface and begin forming a volcanic edifice. As the volcano grows because of melt accumulation at the top of the plate, the lithosphere is flexed downwards, inducing extra tensile stress at 30–35 km depth and compression at ∼15 km depth (induced by the edifice load). It implies that now the melts pond at ∼15 km and form another magmatic reservoir lying just underneath the crust. These processes explain the ponding of primary (shield) melts at ∼35 km and ∼15 km depths as recorded below La Reunion, Mauritius or Hawaii volcanoes, all shifted by ∼150 km with respect to the plume axis

Do faults trigger folding in the lithosphere?

A number of observations reveal large periodic undulations within the oceanic and continental lithospheres. The question if these observations are the result of large- scale compressive instabilities, i.e. buckling, remains open. In this study, we support the buckling hypothesis by di- rect numerical modeling. We compare our results with the data on three most proeminent cases of the oceanic and continental folding-like deformation (Indian Ocean, Western Gobi (Central Asia) and Central Australia). We demon- strate that under reasonable tectonic stresses, folds can de- velop from brittle faults cutting through the brittle parts of a lithosphere. The predicted wavelengths and finite growth rates are in agreement with observations. We also show that within a continental lithosphere with thermal age greater than 400 My, either a bi-harmonic mode (two superimposed wavelengths, crustal and mantle one) or a coupled mode (mono-layer deformation) of inelastic folding can develop, depending on the strength and thickness of the lower crust

Interseismic deformation at subduction zones investigated by 2D numerical modeling: case study before the 2010 Maule earthquake

We study interseismic deformation preceding the Mw8.8 2010 Maule earthquake by means of two-dimensional finite-element modeling. Our goal is to gain insight into the fundamental factors controlling elastic strain build-up and release in subduction zones, and to evaluate different modeling approaches of surface displacement as observed by GPS. We developed a linear elasticity solver that allows us to implement a realistic subducting plate geometry constrained by geophysical data. We test the influence of subducting plate thickness, variations in the updip and downdip limit of a 100% locked interplate zone, elastic parameters, and velocity reduction at the base of the subducted slab. We compared our modeled predictions with interseismic GPS observations along an EW profile crossing the Maule earthquake rupture area, in order to determine best fitting parameters. Our results indicate little influence of the subducting plate thickness at a given downdip limit, which itself has a strong influence on surface deformation. However, the fit to observations is achieved only after reducing the velocity at the base of the subducted slab below the trench region to 10% of the far-field convergence rate. We link this novel result to complementary numerical models that gradually evolve toward considering longer time-scales and complex rheology in order to evaluate the mechanical meaning of the above mentioned inferred kinematic conditions. This allowed us to link the velocity reduction at the base of subducting slabs with a long-term state of high flexural stress resulting from the mechanical interaction of the slab with the underlying mantle. Even a small amount of theses high deviatoric stresses may transfer towards the upper portion of the slab as strain energy that could participate into the mechanical loading of the megathrust and therefore in triggering large earthquakes there

Flat versus normal subduction zones: a comparison based on 3-D regional traveltime tomography and petrological modelling of central Chile and western Argentina (29°–35°S)

Our study compares the seismic properties between the flat and normal subduction regions in central Chile, to better understand the links between the slab geometry, surface deformation and the deeper structures. In comparison with previous studies, we show the most complete 3-D regional seismic tomography images for this region, in which we use (1) a larger seismic data set compiled from several short-term seismic catalogues, (2) a denser seismic array allowing a better resolution of the subduction zone from the trench to the backarc and into the upper ∼30 km of the slab and (3) a starting 1-D background velocity model specifically calculated for this region and refined over the years. We assess and discuss our tomography results using regional seismic attenuation models and estimating rock types on the basis of pressure and temperature conditions computed from thermomechanical models. Our results show significant seismic differences between the flat and normal subduction zones. As expected, the faster seismic velocities and increased seismicity within the flat slab and overriding lithosphere are generally consistent with a cooler thermal state. Our results are also consistent with dehydration of the mantle above the subducted Juan Fernandez Ridge at the eastern tip of the flat slab segment, indicating that the latter retains some fluids during subduction. However, fluids in the upper portion of the flat slab segment are not seismically detected, since we report instead fast slab seismic velocities which contradict the argument of its buoyancy being the cause of horizontal subduction. The forearc region, above the flat slab, exhibits high Vs and very low Vp/Vs ratios, uncorrelated with typical rock compositions, increased density or reduced temperature; this feature is possibly linked with the aftershock effects of the Mw7.1 1997 Punitaqui earthquake, the flat slab geometry and/or seismic anisotropy. At the surface, the seismic variations correlate with the geological terranes. The Andean crust is strongly reduced in seismic velocities along the La Ramada–Aconcagua deformation belt, suggesting structural damage. Slow seismic velocities along the Andean Moho match non-eclogitized hydrated rocks, consistent with a previous delamination event or a felsic composition, which in turn supports the extent of the Chilenia terrane at these depths. We confirm previous studies that suggest that the Cuyania terrane in the backarc region is mafic and contains an eclogitized lower crust below 50-km depth. We also hypothesize major Andean basement detachment faults (or shear zones) to extend towards the plate interface and canalize slab-derived fluids into the continental crust

Field observations and numerical models of a Pleistocene-Holocene feeder dyke swarm associated with a fissure complex to the east of the Tatara-San Pedro-Pellado complex, Southern Volcanic Zone, Chile

Magma is transported through the lithosphere as dykes which, during periods of unrest, may feed eruptions at the surface. The propagation path of dykes is influenced by the crustal stress field and can be disturbed by crustal heterogeneities such as contrasting rock units or faults. Moreover, as dykes propagate, they themselves influence the surrounding stress field through processes of stress transfer, crustal deformation and seismic failure. The result is the formation of arrested dykes, as well as contrasting strike and dip angles and dyke segmentation. Here, we study the mechanisms of dyke injection and the role played in modifying the stress field and potential propagation paths of later dyke injections. To do this we combine field data from an eroded and well-exposed shallow feeder dyke swarm with a suite of two-dimensional FEM numerical models. We mapped 35 dyke segments over a ~1 km long dyke swarm exposed ~5 km to the East of Pellado Volcano, in the Tatara-San Pedro-Pellado (TSPP) volcanic complex, Southern Volcanic Zone of the Andes. Detailed mapping of the swarm elucidates two preferential strike orientations, one ~N80°E and the other ~N60°E. Our numerical models simulate both the TSPP volcanic complex and the studied dyke swarm as zones of either magmatic excess pressure or as a rigid inclusion. The crustal segment hosting the volcanic complex and dykes is modelled using an elastic domain subjected to regional compression in select model cases. Model outputs provide the stress and strain fields resulting from the different geometries and applied boundary loads. The model results indicate that individual dyke injections can locally rotate the principal stresses such as to influence the range of orientations over which later dykes will form. The orientation of σ1 at the dyke tip ranges over 60° (±30° either side of the dyke tip) indicating that the strike orientation of later dykes will fall within this range. The effect of adding a bulk regional compression is to locally increase the magnitude of favorably oriented tensile stresses in the bedrock but to reduce the range of σ1 orientations to 40° (±20°). This implies that under a far-field transpressive stress regime, as is common in Andean settings, regional dyke swarms will tend to maintain their strike orientation parallel to the regional bulk stress. These results should be accounted for when studying periods of volcanic unrest in order to discern the location and orientation of potential fissure eruptions in active volcanic areas such as the Southern Volcanic Zone of the Andes