An Assessment of Stress States in Passive Margin Sediments: Iterative Hydro-Mechanical Simulations on Basin Models and Implications for Rock Failure Predictions

: Capturing the past and present hydro-mechanical behavior of passive margin sediments raises noticeable interest, notably in geo-hazard risk assessment and hydrocarbon exploration. In this work, we aim at assessing the stress states undergone by these sedimentary deposits through geological time. To do so, we use an iterative coupling between a basin simulator and a ﬁnite element mechanical solver. This method conciliates a computation of the full stress tensors with a dynamic and geologically detailed modelling of the sedimentation. It is carried out on a dedicated set of 2D synthetic basin models, designed to be representative of siliciclastic deposition in passive margins and integrating variations in their geological history. Contrary to common assumptions in operational basin modelling studies, our results imply that passive margin sedimentary wedges are multidimensional mechanical systems, which endure signiﬁcant non-vertical stress without external tectonic input. Our results also highlight the variability of the stress states through space and time, with a strong control from the geometry and lithological heterogeneities of the deposits. Lastly, we used the simulation results to predict a location and timing for the development of weakness zones in the sedimentary stacks, as privileged areas for rock failure. The outcome underlines the inﬂuence of the basal tilt angle, with a slight tilt impacting the wedges stability to a similar extent as a substantial increase in sedimentation rate. Altogether, this study emphasizes the need for careful consideration of non-vertical stresses in basin simulations, including in passive tectonic contexts. It also suggests that the iterative coupling method employed is a promising way to match industrial needs in this regard.


Introduction
Passive margins represent a cumulated length of 105,000 km of coastal areas around the world and are privileged locations for marine sediment deposition [1]. They are well-studied environments, notably in terms of morphology, geology, and underground resources [2,3]. On the latter topic, passive margins are historically key areas in the hydrocarbon industry, as they host more than 300 giant oil and gas fields [4].
Passive margin sediments are subject to a large range of hydraulic and mechanical phenomena, which creates specific interests. Remarkably, understanding their hydro-mechanical behavior is critical in civil engineering and geo-hazard prediction, as passive margin sediments are naturally likely to collapse at local or regional scales [5,6]. In petroleum exploration, analyzing slope instability and mechanical compaction through geological ages helps to accurately estimate the location of distal sand reservoirs and their volumes at the present time [7,8]. Furthermore, in all conventional petroleum systems, hydrocarbon migration is controlled by gradients of fluid overpressure and impacted by natural fracturing of cap rocks [9,10]. Lastly, assessment of present-day pore pressure appears as a crucial step to ensure drilling safety [11,12].
In the present study, our objective is to assess the stress states undergone by siliciclastic sediments in passive margins and the development of weakness zones as favorable locations for rock failure. To do so, our simulation strategy is rooted on an iterative hydro-mechanical coupling between a basin simulator, ArcTem [13], and a finite element mechanical solver, Code Aster [14]. The hydro-mechanical coupling was applied on a set of synthetic basin models, designed with the target to be representative of siliciclastic sedimentation in passive margins and of operational models in the hydrocarbon industry. This aims at depicting the stress states endured by the different parts of the sedimentary wedges through deposition, with an emphasis on non-vertical effects. Our goal is also to appraise how basal slope tilting and sedimentation rate variations affect the location and expansion of weakness areas, thus addressing the geological control on rock failure preconditioning through passive margins history.

Observations
Passive margins are locations of siliciclastic sediment accumulations, up to 12 km thick [15]. These accumulations usually consist of a stacking of several stratigraphic units, following well-described patterns in terms of geometry and lithology [16]. Most layers present a clinoform shape, with maximal slopes ranging from 1 to 10 • in siliciclastic environments [17].
Passive margin deposits are highly compactable and experience significant volume reduction, but important under-compaction and fluid overpressure are also observed [18]. This excess of pore pressure is notably induced by high sedimentation rates and low permeability layers, keeping the hydrodynamic system out of equilibrium [19,20]. Sequence stratigraphy patterns, controlling the relative distribution of sealing shaly facies and draining sandy ones, reflect on the pore pressure profiles and cause substantial lateral flows [21,22].
Sedimentary wedges in passive margins are dynamic and unstable systems, leading to failure events of different scales and amplitudes. Locally, buried and overpressured sediments are prone to natural hydraulic fracturing, known to damage the sealing capacity of shaly layers [23,24]. At a larger scale, instability can lead to the failure of the sedimentary slope, triggering important mass transport, from superficial sliding to wider slumping and even broad landslips [25][26][27][28][29][30]. Lastly, at a regional scale, gravity collapse of passive margin sediments set common structural features, with normal faulting observed in the proximal part of the wedge and fold-and-thrust belts in the most distal areas [31,32].
Although punctual and exceptional events are often needed to trigger the amplest incidents, failure at all scales is favored by preconditioning factors resulting in weakness areas in the sedimentary stack [5,6,33]. Pore pressure rise is noted as a key factor, in relation to sealing lithology or high deposition rates [34][35][36][37]. However, the geometrical evolution of the wedge is also reported to play a major role in its destabilisation, with an influence of the clinoform slopes [38,39], progradation rates [40], and basal angles [41,42]. Lastly, mechanical heterogeneities can be decisive, with regional collapses facilitated by deep low-friction layers, for instance evaporites [43] or overpressured shales [44].

Traditional Approaches
Several methods are described in the literature to model the hydro-mechanical behavior of passive margin sediments. They can be sorted into two main families.
A first approach relies on slope stability analysis. It covers a range of methods of increasing complexity, including critical Coulomb wedge theory, infinite slope and limit analyzes, their combinations, and their extensions [45][46][47][48][49][50][51]. These methods use slope and pore pressure estimates to compute the shear stress and the rock strength, leading to a likelihood of failure for the wedge. However, they are often based on a strong idealization of the sedimentary stack, notably in terms of geometry, lateral facies variations, and overpressured layers [6]. In addition, this kind of modelling is essentially static. Even when it is implemented in a sequential manner, each time step remains independent from the previous ones, with no consideration of fluid flows and past stress states [52].
A second approach relies on basin simulation technologies [53][54][55][56][57][58][59]. From a representation of the basin at the present day, the evolution of its geometry through geological ages is restored with a backward process. Then, the physical variables of interest are computed for each geological event with a forward simulation using this geometrical framework. This methodology can handle a high degree of realism in the description of the sedimentary deposits, notably in terms of slope breaks, truncated layering, lateral lithological heterogeneities, and complex overpressure distribution. The resulting models are also intrinsically dynamic: the simulation results at a given geological age are strongly dependent on the flow and stress history of the margin, and they are also significantly impacted by the evolutive geometry previously restored [60].
Nonetheless, using basin simulation methodologies, it is classical to reduce the stress state to its sole vertical component, to overcome numerical complexity or absence of calibration data [56,57,61,62]. Implicitly, the contribution of lateral efforts to porosity loss, fluid overpressure build-up and system instability is then estimated as marginal. Horizontal stresses are often considered as direct functions of the burial depth or simple ratios of the vertical stress. In the latter case, the ratios are homogenous in the sedimentary stack or solely variable with depth and lithology [57,63].
The limitations of the usual stress state simplification and their impact on hydro-mechanical simulation results have already been documented [64,65]. Notably, they can imply incorrect estimates for sediment compaction and rock failure likelihood. In a petroleum exploration context, this may produce misleading conclusions on cap rock integrity, migration flowpaths, and present-day overpressure. These observations motivated the design of several alternative methods aiming to better account for the 3D nature of the stress tensor [66,67]. Nevertheless, the focus is traditionally on basins overcoming important lateral tectonic solicitations, with large shortening or extension rates. As a result, the value of integrating non-vertical stresses in the modelling of passive geological settings appears more questionable, especially when considering the associated increase in computation costs. Basin modellers often deal with passive margin sediments as one-dimensional mechanical systems through geological ages [68][69][70], in line with simplistic models for horizontal stresses computation and the assumption that the maximal principal stress is perfectly aligned with the vertical axis [18].

Iterative Hydro-Mechanical Coupling
In this work, we aim at conciliating the geological realism and dynamic aspect of basin simulations with the integration of lateral mechanical effects. To do so, we use a coupling between two codes: ArcTem, a basin simulator [13], and Code Aster, a finite element mechanical solver [14]. ArcTem relies on the classical 1D stress simplification while Code Aster computes a full stress tensor.
The coupling is implemented in a prototype software named A 2 [71]. Following an iterative procedure, the two codes exchange inputs and outputs at each simulation step, until providing a consistent set of porosity, permeability, fluid pressure, and stress tensor in each cell of the basin model. This differs from the integrated numerical model of [72] but shares similar objectives.
The coupling algorithm uses an updated Lagrangian approach and is inspired from the work of [73] and [74] for reservoir simulations. The same grid is used by the two codes and the convergence criterion is based on a comparison between the porosity values they compute. A full description can be read in [71] and [75].

Rheology
The ArcTem simulator uses a direct relationship between porosity and vertical effective stress, known as the Schneider's law [76,77]. This law is widely used in the basin modelling community for its practicability, notably to reproduce porosity-depth data measured from cores or well logs [57].
In the iterative coupling, we model the rheology of passive margin sediments with a tensor extension of the Schneider's law. Concretely, the model describes an elasto-plastic behavior, with the rock stiffness varying with compaction. It accounts for the full stress tensor but is designed to trigger the same volume variations as the Schneider's law in the case of tabular sedimentation with no lateral deformation (i.e., oedometer test conditions). It is presented in detail in [71]. Other relevant models describing sediment compaction can be found in the literature [78,79]. However, we see the one used herein as a suitable compromise with several practical needs, notably applicability on geologically detailed passive margin models, easy comparison with classical basin simulators, and integration in the iterative coupling at a reasonable computation cost.

Set of Synthetic Basin Models
This study relies on a dedicated set of five synthetic 2D basin models, designed to be representative of actual siliciclastic sedimentation in passive margins. They share a common design, but differ in geological history, with variations in basal slope tilt and sedimentation rate.

Geological Setting
Our modelling work focuses on a few stratigraphic sequences of siliciclastic deposition. It aims at representing the sedimentary dynamics classically observed in a passive margin, notably the morphology of stratigraphic system tracks, with gradating deposition, discordant stratification, and lateral facies' variations. These geological features are commonly integrated in operational basin models of the petroleum industry, but more rarely in published hydro-mechanical simulations, as they often involve a higher level of idealization.
Among other inspirations, the modeled sequences can be related to the sediments of Eocene age in the Brazilian oil-bearing Santos basin, as described in [80]. In this example, they were deposited in a relatively short amount of space and time (30 km and 12 My). More generally, the modeled sequences represent common sedimentation patterns, but the accumulations in passive margin basins usually consists in their repetition through longer time periods and deposited volumes. Consequently, the structural evolution of passive margins at large time or space scales (i.e., hundreds of km and My) is beyond the scope of this work.

Common Features
All models in the set correspond to a 33 km long section and a succession of 40 gradating sedimentary layers, distributed in stratigraphic sequences of transgressive, highstand, lowstand, and turbiditic system tracks ( Figure 1). The sedimentation is globally prograding and most layers present a clinoform morphology, with a slope dipping up to 10 • . This angle refers to the highest value measured in clastic deposits on passive margins around the world [17]. Maximum sediment pile thickness is 3 km. Under the sediments, 2 km of a stiff geo-material are previously set up to create a hard basis to the model.

Variations between Models
Keeping a common design, we introduced different geological histories for the models in the set. Our objective is to appraise the geological control on failure preconditioning in passive margin sediments, and to assess the development of weakness areas in terms of timing and expansion. Among the numerous geological factors varying through a margin history, we specifically focused on the sedimentation rate and on the angle of the basal slope, which usually undergoes a progressive tilt during the subsidence of the sedimentary stack.
High sedimentation rates favor fluid overpressure build-up, and in the literature previous physical simulations illustrate how rapid deposition can create weakness zones in deeply buried parts of the sedimentary wedges [35,46]. However, weakness and rock failure are also observed in margins with slow sedimentation, and other published simulation studies showed that additional preconditioning factors than overpressure should be considered [33,49]. In this regard, the models presented in this work are a way to specifically evaluate the impact of the basal angle tilt during sedimentation. Table 1 presents the three sedimentation rate scenarios simulated in this study. In each case, a specific duration is assigned to all sedimentation events regardless of their thickness, which makes the sedimentation rate vary with time. However, the resulting values are consistent with the ones estimated on actual passive margins around the world [26,81,82]. Table 2 presents the three tilting scenarios simulated. In every model, the bottom horizon starts horizontal and then tilts at an increment angle at each transgressive event, until reaching the target value in the present day.  Table 2. Values of subsidence tilt angle in the models set.

Present-Day Grids
For all models, a grid was built from the present-day geometry by dividing the 40 layers by 70 vertical pillars. The total number of cells is then 2760, but 47% of them are flattened because of the discordant nature of the sedimentation.

Past Geometries
The restoration of past geometries for all models was performed with a backstripping method [83,84], implemented in the TemisFlow TM software (2016 version, Beicip-Franlab, Rueil-Malmaison, France) [85]. The uncompacting used a different Porosity vs. Depth table for each lithological type. Shalier facies have larger deposit porosities but compact more quickly with burial than the sandier ones (plots available in Figure A1 in the Appendix A).
The model bottom was used as a reference horizon. Its shape was arbitrarily input for each deposition event in order to model a progressive subsidence and tilting for the margin basement, following Table 2. The paleo-bathymetry profiles were then computed by the algorithm from the compaction tables and the basement shape. Lastly, absolute sea level variations were input by translating the full model upward or downward from the 0 mark on the vertical axis. These variations were set in accordance with the sedimentary track being deposited, and keeping the full model immersed during the whole simulation.

Simulation Parameters
Once built and restored, all models were simulated with the same parameters.

Mechanical Properties
The Schneider's law and its tensor extension were used with specific coefficients for each lithological type, in order to model different compaction behaviors (plots available in Figure A1 in the Appendix A). These coefficients were all calibrated to be perfectly consistent in hydrostatic conditions with the Porosity-Depth tables used in the backstripping restoration. In addition, the whole model was given a solid density of 2.7 g·cm −3 and a Poisson's ratio of 0.32, seen as average values for well-consolidated siliciclastic sediments [86,87].

Permeability
Different Permeability-Porosity tables were associated with each of the six facies (plots available in Figure A2 in the Appendix A). These tables were designed to be representative of the actual petro-physical properties of layered siliciclastic sediments at the basin scale. Insights from various references were combined, notably experimental data and empirical relationships for shale rocks [88,89] and sandstones [90,91], literature review [92], and basin-modelling guidelines [57].
For the two sandier facies, permeability decreases with compaction following a classical convex-shaped profile. Still, for the three shalier facies, a high permeability drop is modelled between 0.35 and 0.2 porosity to represent the mineralogical transformations observed during the burial of shaly sediments, with kaolinite-rich clays converting to less permeable illite-rich ones [93]. In addition, permeability anisotropy was added to the model to represent the stratification of the sedimentation. The ratios between horizontal and vertical permeability range from 1 for the sandiest facies to 18 for the shaliest one, in accordance with diverse data available in the literature [57,94,95].

Boundary Conditions
When running the iterative coupling, the geo-mechanical solver was strongly constrained by the backstripping restoration, and no external tectonic loading was added. Practically, the input boundary conditions were fixed nodal displacements for the bottom horizon and normal displacements blocked on lateral boundaries.
The basis was given a very stiff behavior with very little deposit porosity, similar for instance to an intact Barre granite [96]. Its permeability was defined accordingly, from several laboratory measurements on such granites [96][97][98].

Simulation Results
In this section, we analyze the simulation results in three steps. First, we observe on the base case model the classical basin modelling outputs which are porosity and pore pressure, to verify their consistency with our modelling objectives. Second, we focus on the additional output provided by the iterative hydro-mechanical coupling technology: a computation of the full stress tensor in every cell of the model at every geological age. We characterize the multidimensional nature of the simulated stress states and appraise their variability through space and time. Lastly, we use the stress tensors simulated to predict the development of weakness areas in the sedimentary stack, seen as favorable locations for rock failure. We benefit from the full set of basin models to compare the influence of two geological factors impacting the sedimentary wedge stability: deposition rate and basal angle tilting. Figure 3 depicts the simulation results on the base case model in terms of present-time porosity and fluid overpressure evolution through geological ages. We can assess whether compaction and overpressure levels are sufficient to be representative of actual passive margins, and whether their changes through space and time can be related to the lithology and the geometry of the modelled deposition.

Pore Pressure and Porosity
Despite moderate burial, strong overpressure is simulated, with a maximum of 21.8 MPa in the present day. The highest values are located in a proximal and deep area, where the overlying sediment stack is the thickest. The simulated overpressure rises in this area during the deposition of highstand system tracks, increasing its burial. During the deposition of lowstand system tracks, the simulated overpressure rather rises beneath the continental slope, in the previously deposited lowstand wedges. The lateral facies variations then reverberate on the pressure profile: the simulated overpressure build-up concentrates in the inferior part of the wedges, presenting a shaly lithology, while the sandier superior part remains well drained.
Meanwhile, the simulated sediment compaction is significant but not extreme, with the lowest present-day porosity value at 25%, which corresponds to approximately half of the value at deposition. For a given burial, porosity estimations reflect both the lithological heterogeneities and the pore pressure profile, with higher values in the areas of sandy lithology or strong overpressures. At first-order and considering the strong differences in sediment stack thickness, the simulation results in the present day appear consistent with those previously predicted in real passive margins from seismic data, in situ measurements or idealised basin models, notably in New England [18,34,99], the Gulf of Mexico [36,100] or the Storegga slide area offshore Norway [46].

Stress States
With the iterative hydro-mechanical coupling, a stress tensor is computed in each cell of the basin model at each geological event. It represents the total stresses endured by the passive margin sediments during this event. In this study, we note S as this stress tensor, X the distality axis oriented from onshore to offshore, and Z the depth axis oriented downward. Consequently, S XX corresponds to the horizontal stress, S ZZ to the vertical stress, and S XZ to the shear stress. We use the sign convention with positive compressive stresses. The stress tensor is then diagonalized to determine the principal stresses. We note S MAX the maximum principal stress and S MIN the minimum one. In addition, effective stresses are computed from total stresses and pore pressure values using Terzaghi's theory [101]. The resulting tensor is noted S .
We note p the confinement pressure, or average effective stress, computed from S : With the confinement pressure, it is possible to isolate the deviatoric part Q of the effective stress tensor: In rock mechanics, the multidimensional nature of the stress states is then commonly characterized by the equivalent deviatoric stress q, defined as follows: The value of this equivalent stress is not necessarily meaningful on its own, and should be compared to the confinement pressure, notably in rock failure predictions. However, from a physical standpoint, a rise in deviatoric stress can be related to two phenomena: (1) development of shear stresses, leading to the rotation of the principal stress system, and (2) higher anisotropy in the principal stresses, leading to a more strongly extensional tectonic regime [102]. Consequently, in the following, we evaluate the significance of these two mechanical effects in our models and consider their variability through space and time, in relation to the lithology and geometry of the sediments deposited.

Shear Stress and Rotation of Principal Stresses
The shear stress in the present day simulated on the base case model is presented in Figure 4. The highest absolute values are concentrated in a deep area, located beneath the center of the continental slope. These values are quite moderate, as they do not exceed 3.5 MPa (Figure 4a). However, they appear much more significant when expressed as a ratio of the vertical effective stress. On the base case model, the shear stress in the present day represents up to 40% of the vertical effective stress (Figure 4b). The greatest ratios are simulated in the area of strongest shear but also in the inferior part of the overlying lowstand wedges, where the vertical effective stress is lowered by high overpressure. Shear stresses reflect on the local rotation of the principal stress system, measurable by the angle from vertical for the orientation of the maximal principal stress. Figure 4c presents these rotations in the margin sediments through geological time, as simulated on the base case model. Rotations up to 32 • are observed. The highest rotation angles are found in the area of highest absolute shear, deep beneath the center of the continental slope. However, strong rotations are also simulated in much shallower areas, due to low vertical stresses and sloping geometries for the deposits. These shallow zones of principal stress rotation mostly appear downslope during the deposition of the lowstand system tracks, and upslope during the deposition of highstand system tracks.

Extension and Compression
Even if substantial rotations of the principal stresses are simulated, the angle values displayed in Figure 4c show the normal efforts remains higher than the tangential ones in every part of the model. This means the tectonic regime of all sediments in the model remains extensional through their history. However, the ratio between minimal and maximal principal stresses varies with space and time, delineating areas and ages where the regime nears compression. These variations simulated on the base case model are represented in Figure 5. A clear zoning is visible in the simulation results: the proximal part of the model hosts the lowest principal stress ratios and the most extensional regimes, while the distal part gets very close to compression with ratios above 95%. This is in line with the structural style commonly observed after the collapse of passive margin sediments. Extensional normal faults usually appear in a proximal section of the wedge and compressional fault-and-thrust belts in a distal one. However, if the principal stress ratios simulated in the distal area remain quite constant through time, those in the proximal area varies significantly with the stratigraphic units being deposited. They appear significantly lower during the lowstand periods, with large areas under 80%, than during the highstand periods, with most of the model above 85%. This can be related to a more proximal and deeply immersed sedimentation during the highstand ages.

Weakness Prediction
The simulated stress states can be used to estimate a timing and spatial distribution for the development of weakness areas in the sedimentary wedge, considered as privileged locations for local natural fracturing or wider failure initiation. In this study, we introduce a weakness criterion based on the work of [103]. As presented in [71], this criterion corresponds to the shear plasticity threshold of the rheology used in the iterative hydro-mechanical coupling and can be assimilated to a critical state line. It is defined from the confinement pressure and the equivalent deviatoric stress. Concretely, a sediment is considered weak when: A and B are usually deduced from the internal friction angle and cohesive strength of the rock material. In this work, we used A = 1 (friction angle equal to 25 • ) and B = 0.5 MPa.
At each geological age in the model, the simulation results are post-processed to highlight the areas above the weakness criterion. It is worth mentioning the permeability of the sediments in these areas remain unchanged, as the criterion is considered to only represent a favorable pre-condition to failure and not directly natural fracturing. However, the corresponding cells in the model are tagged for the rest of the simulation, thus recording the progressive expansion of the weakness in the wedge. Practically, this methodology is only applied on the sediments qualified as sufficiently consolidated, which is previously defined with a porosity cut-off of 35%.
The results on the base case model are presented in Figure 6. The weakness zone in the present day corresponds to 144 cells of the base case model, covering approximatively 10% of the sedimentary wedge. Weakness appears quite late in the sedimentation history, at the 24th event out of 40. It initiates deep beneath the center of the continental slope, in a location of strong shear (visible in Figure 4). Next, the weak area grows downward and upward, forming a stripe roughly parallel to the continental slope, which seems in line with the slip plane orientations commonly observed in passive margin sediments [5,35]. Afterwards, a second weakness stripe initiates at the 33th event, also beneath the center of the continental slope, but in a shallower location. Then, both weakness areas expand through time, until joining in a single continuous weak zone. This result suggests that several levels of different shallowness can simultaneously be preconditioned to failure. This seems also in line with field observations, where superficial slides can coincide with more deeply rooted instabilities [6,26].

Influence of Sedimentation Rate
The influence of sedimentation rate on failure preconditioning is addressed with the variations introduced between the different basin models of the set. The iterative hydro-mechanical coupling was run on the three models described in Table 1 and the development of weakness areas predicted from the simulated stress states. The results are presented in Figure 7. As expected, a sedimentation rate increase reflects on higher fluid overpressure. For instance, the strongest present-day overpressure simulated in the high rate scenario is 22.7 MPa when it only reaches 20.5 MPa in the low rate one. This leads to a weakness area of earlier initiation and wider expansion. Weakness appears at the 15th event and covers 186 cells in the high rate scenario, when it initiates 10 events later and only covers 88 cells in the low rate one. These differences between scenarios can be explained physically. With high sedimentation rates, the increase in pore pressure engenders a decrease in confinement pressure, while the shear stresses remain identical. Consequently, the ratio between confinement pressure and deviatoric stress rises compared with lower sedimentation rate scenarios, which results in the sediments getting closer to the weakness criterion.
Our physical interpretation can be confirmed by local analyses-for instance, tracking the ratio between confinement pressure and deviatoric stress through time for a given cell and comparing scenarios. Such analysis was carried out for a cell corresponding to deep shaly sediments directly overlying the deepest turbiditic deposits (Figure 8a). In a petroleum exploration context, these sediments could form a sealing cap rock above a sandy reservoir. Assessing their likelihood of failure through geological ages would then be critical to appraise the economic interest of the prospect. The output is displayed in Figure 8b. The decrease in compaction originated by the increase in pore pressure is then clearly visible, with higher sedimentation rates conducing to stress paths shifted leftward in the p -q system. This leads the corresponding sediments to reach the weakness criterion at an earlier age in the high rate scenario than in the base case one, while the criterion is never reached in the low rate hypothesis.

Influence of Subsidence Tilt
Thanks to the variations in subsidence tilt angle introduced between the different basin models of the set, the influence of the basal slope on failure preconditioning can be compared with the influence of sedimentation rate. This aims at providing insights on the instability likelihood in situations of slow sedimentation, and at appraising the relative weight of geometrical factors (such as basal angle, clinoform slopes, or progradation rates) and hydro-dynamical ones (such as sedimentation rate or seal permeability) in the geological control of weakness development through the margin history. The iterative hydro-mechanical coupling was run on the three models described in Table 2, and weaker zones were delineated at each age from the simulated stress states. The results are presented in Figure 9.  An increase in tilt angle, contrary to an increase in sedimentation rate, barely changes the simulated pore pressure. However, it engenders stronger shear stresses and rotations of the principal stresses. In the high tilt scenario, the absolute value of the shear stress S XZ reaches 3.9 MPa when it does not exceed 3.0 MPa in the low tilt scenario. This also leads to an earlier initiation and a wider expansion of the weakness area. In this regard, the differences between our tilt scenarios are of similar magnitude as the differences between our sedimentation rate scenarios (visible in Figure 7), and even slightly more. Weakness appears at the 15th event and covers 192 cells in the high tilt scenario, when it initiates 10 events later and only covers 73 cells in the low tilt one.
The physical interpretation for the impact of subsidence tilt on weakness development differs from the sedimentation rates one. With high tilt angles, the average effective stress is merely impacted, but the increase in shear stress results in an increase in deviatoric stress. Consequently, the ratio between confinement pressure and deviatoric stress rises compared with lower tilt angle scenarios, which results in the sediments getting closer to the weakness criterion. This different way of reaching failure-prone conditions is confirmed by analyzing the local stress paths in the basal tilt scenarios, with exactly the same methods as done for the sedimentation rate ones. The output of this analysis is displayed in Figure 10. The increase in deviatoric stress originated by the increase in shear stresses is then clearly visible, with higher subsidence tilts leading to stress paths shifted upward in the p -q system. This upward shift has similar consequences as the leftward shift caused by high sedimentation rates. The studied sediments reach the weakness criterion at an earlier age in the high tilt scenario than in the base case one, while the criterion is never reached in the low tilt hypothesis. The simulations presented herein confirm that geometrical factors like the evolution of the basal angle can be decisive in the destabilisation of passive margin sedimentary wedges. Although the influence of sedimentation rate is more often invoked and more thoroughly studied in the scientific literature, our study indicates a slight tilting of the sedimentary stack can be of equal significance in failure preconditioning. This advocates for not overlooking geological and geometrical uncertainties in petroleum exploration or geo-hazard risk management, as they can crucially impact predictions of local natural fracturing or wider collapse for the passive margin sediments.

Stress States in Passive Margin Sediments
Altogether, our simulation results reflect the complexity of stress states in siliciclastic passive margin sediments. Notably, the geo-mechanical effects simulated are far from purely vertical, with significant rotations of the principal stress system and stress regimes nearing compression in the distal part of the wedge. It is noticeable that no external tectonic solicitation is needed to trigger the simulated non-vertical effects. They are only driven by the geometric evolution of the sedimentary wedge, with a moderately sloping clinoform deposition and a slightly tilting basis as the two key features.
Another finding is the variability of the non-vertical mechanical effects through space and time. Their magnitude differs considerably depending on the zone considered in the sedimentary wedge. For instance, in our models, the proximal part of the margin hosts the most extensional regimes, whereas the tangential stresses are close to the amount of the normal ones in the abyssal plain. Meanwhile, the highest shear stresses are concentrated beneath the continental slope. These lateral variations are directed by the geometry of the sedimentary stack and reflect the lithological heterogeneities, inferring a strong geological control on the stress distribution. The magnitude of the non-vertical mechanical effects simulated also changes with the stratigraphic age considered. As an example, the highstand periods coincide with higher principal stress ratios in the proximal area, while the lowstand periods coincide with higher downslope rotation of principal stresses.
The significance and variability of the simulated non-vertical effects suggest their integration in physical simulations worth the computation cost, even when studying passive geological contexts. Our work implies that considering passive margin sedimentation as a one-dimensional mechanical system or assessing non-vertical stresses with too simplistic methods may prove misleading. For instance, simply comparing the simulated pore pressure to 85% of weight of the overlying layers, which is often the default option of industrial basin modelling software, would result in predicting rock failure on the most proximal and deep areas of our models, which host the strongest overpressure values, rather than under the shelf break and the continental slope. Likewise, as the magnitude of non-vertical mechanical effects in our simulations appeared strongly controlled by the changing geometry of the wedge and the represented lithological heterogeneities, modellers should beware of excessive idealization of these two elements considering passive margin sediments.

Geological Control on Rock Failure Preconditioning
The importance of considering non-vertical mechanical effects in passive margin sediments was further illustrated by comparing the relative impact of sedimentation rate and subsidence tilt on our simulations results in terms of weakness development and failure preconditioning. Our models suggest that the rise in deviatoric stress caused by the tilting of the deposits can affect the stability of the wedge as significantly as the pore pressure rise caused by an increase in sedimentation rate. This could explain how rock failure arises in margins of slow sedimentation with relatively well-drained accumulations. However, appraising the impact of this kind of geological and geometrical uncertainty on rock failure predictions cannot be achieved easily with traditional basin models relying on a one-dimensional mechanical scheme. It requires alternative technologies accounting for the 3D nature of the stress states, as the one used in this study.

Modelling Hypotheses and Way Forward
This study relies on several modelling choices, and we find it useful to remind the main ones here. A first series of choices relate to the synthetic basin models used for the simulation. Notably, the slope of the modelled clinoforms corresponds to the steepest siliciclastic wedges actually observed. A second series of choices relate to the physical simulation itself. Average values were picked for the solid density and the Poisson's ratio of the rock material. Thermal effects were neglected for the sake of simplicity and to focus on the hydro-mechanical phenomena. We also chose not to increase the permeability of the sediments reaching our weakness criterion, as it only represents a favorable precondition to failure. Lastly, the rheology law implemented in the coupling remains quite simple and was essentially implemented for its practicability in terms of comparison with classical basin simulations and numerical cost.
The following step in this research would be to apply a similar methodology to an actual passive margin where the results can then be related to abundant geological knowledge. Doing so, the level of geological realism of the models could be further increased by linking the iterative hydro-mechanical coupling with a forward stratigraphic simulator [104]. First attempts in a carbonate platform context are addressed in [105].

Integration of Non-Vertical Stresses in Basin Simulators
The results of this study highlight the value of integrating non-vertical mechanical effects in basin simulations. In our passive margin context, it notably appeared paramount to appraise the impact of geometrical and geological factors like basal angle tilting on the timing and location of failure preconditioning in the sedimentary wedge. In a petroleum exploration context, this would facilitate estimates and risk assessment for reservoirs' location and cap rocks' integrity. Logically, integrating non-vertical mechanical effects is even more crucial in more structurally complex areas, enduring significant lateral tectonic input. Nonetheless, before meeting industrial needs, 3D hydro-mechanical simulators must reach an appropriate compromise between detailed physical models, easy integration in existing basin modelling software, and applicability on operational basin models. This still remains a stimulating topic for research and development activities.
In this work, we appraised a method based on a sequential and iterative coupling between a classical basin simulator and a finite element mechanical solver, as described in [70] and [74]. It showed satisfactory applicability on our set of 2D synthetic basin models, designed to be representative of siliciclastic sedimentation in passive margins, including laterally gradating deposition, discordant stratification, and lateral variations of lithological and mechanical properties. It also proved to be moderately expensive in computation cost as the total simulation time for a given model did not exceed 15 minutes with a standard workstation. Consequently, we consider this method as a promising way to capture and incorporate 3D mechanical effects in basin simulations, notably in an industrial oil and gas exploration context. Further appraisals must include heavier and 3D basin models as well as simulation of hydrocarbon maturation and migration. Figure A1. Porosity for the backward restoration and the forward basin simulator. At significant burial depths (e.g., 1 to 8 km), the shalier the facies, the less porous the sediment. However, this ranking is inverted at deposition.