# Spontaneous Initiation and Evolution of Polygonal Fault Systems During the Early Burial of Sediments in 3D Numerical Models

Alexandre Chemenda, G Ballas, A Gay

► To cite this version:

Alexandre Chemenda, G Ballas, A Gay. Spontaneous Initiation and Evolution of Polygonal Fault Systems During the Early Burial of Sediments in 3D Numerical Models. Journal of Geophysical Research: Solid Earth, 2025, 130 (7), pp.e2025JB031506. 10.1029/2025jb031506. hal-05300425

HAL Id: hal-05300425

https://hal.science/hal-05300425v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2025JB031506

**Key Points:**

* Numerical models show how diagenetically driven volume reduction in buried sediments creates depth-dependent polygonal fault pattern

* Fault spacing and maximum throw correlate with sediment layer thickness; fault displacement increases with volume reduction

* Model results suggest <10% volume reduction can explain typical fault attributes, matching natural observations

**Correspondence to:**

A. I. Chemenda,
chemenda@geoazur.unice.fr

**Citation:**

Chemenda, A. I., Ballas, G., & Gay, A. (2025). Spontaneous initiation and evolution of polygonal fault systems during the early burial of sediments in 3D numerical models. *Journal of Geophysical Research: Solid Earth*, 130, e2025JB031506. https://doi.org/10.1029/2025JB031506

Received 5 MAR 2025

Accepted 25 JUN 2025

**Author Contributions:**

**Formal analysis:** A. I. Chemenda

**Investigation:** A. I. Chemenda

**Methodology:** A. I. Chemenda

**Software:** A. I. Chemenda

**Validation:** A. I. Chemenda

**Visualization:** A. I. Chemenda

**Writing – original draft:** A. I. Chemenda, G. Ballas, A. Gay

**Writing – review & editing:**

A. I. Chemenda, G. Ballas, A. Gay

<figure><img src="image_2.png" /></figure>

# Spontaneous Initiation and Evolution of Polygonal Fault Systems During the Early Burial of Sediments in 3D Numerical Models

A. I. Chemenda¹, G. Ballas², and A. Gay²

¹Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, IRD, Geoazur, Valbonne, France, ²Géosciences Montpellier, Université de Montpellier, CNRS, Montpellier, France

**Abstract** Polygonal Fault Systems (PFS) or networks were generated in 3D finite-difference models within the upper (active) layer of two-layer model. The driving force of this process results from a progressive diagenetically induced volumetric contraction (porosity reduction) of the active layer during its burial. This contraction causes extensional strains, reduces horizontal compressive stresses, and increases deviatoric stresses, leading to elastoplastic yielding and strain-softening of the material in this layer. At a certain point, the initially homogeneous material loses stability resulting in deformation bifurcation and localization within narrow deformation bands (incipient faults). The material undergoes progressive failure within the bands, accompanied by the accumulation of the normal-sense displacement along the faults and the evolution of the fault system, with some faults dying and others forming. The fault architecture characterized in the models by depth-dependent fault density and pattern, is remarkably similar to the natural PFS, both buried and exhumed. The fault spacing (S) in the models scales with the active layer thickness (T) and increases as shear coupling ($\tau_0$) between the active and substratum layers decreases. The maximum throw $\delta$ along the faults linearly scales with T and increases with the volume or porosity reduction $\Delta\phi$. Under the chosen model parameters, the typical $\delta/T$ ratio of 0.045 for natural PFS, was obtained in the models for $\Delta\phi \approx 7\%$. However, a discrepancy exists between the location (depth) of the $\delta$ maximum in the model and nature, which is discussed in the paper and indicates avenues for further modeling development.

**Plain Language Summary** When fine-grained sediments are buried beneath the seafloor, they undergo natural changes that cause them to shrink in volume. Our research shows how this shrinkage can create regular patterns of fractures or faults in the sediments, forming polygonal shapes when viewed from above that change with depth. Using 3D computer models, we simulated how these polygonal fault systems (PFS) develop naturally during burial, without any external (tectonic) forces. Our models show that as sediments contract, they develop a network of faults that looks remarkably similar to those observed in nature, both in modern seafloor and subsurface sediments as well as in ancient rocks exposed on land. We found that the spacing between faults is directly related to the thickness of the sediment layer, and the amount of movement along these faults depends on how much the sediments shrink. This understanding is important because these fault systems can affect the flow of fluids through sediments, which is crucial for oil and gas exploration, CO₂ storage, and understanding how water moves through ocean sediments. Our simplified approach successfully reproduces the complex patterns seen in nature, helping us better understand how these intriguing geological structures form.

## 1. Introduction

Polygonal fault systems (PFS) have been documented in many basins worldwide (Cartwright, 1994; Clausen et al., 1999; Gay et al., 2004; Ghalayini & Eid, 2020; He et al., 2010; Lonergan et al., 1998; Sun et al., 2010; Turrini et al., 2017; Zeng et al., 2022). PFS typically affect sedimentary units composed of fine-grained, predominantly clay-rich (Cartwright, 1997), but also carbonate-rich (Hansen et al., 2004) sediments. These systems are usually restricted to stratigraphic intervals (tiers) a few hundred meters thick (~100–1,000 m), controlled by lithological and diagenetic variations in the host sediments (Cartwright & Dewhurst, 1998; Xia et al., 2022), and they abruptly terminate at specific stratigraphic horizons (Gay et al., 2007).

The same fine-grained and low-permeability sedimentary sequences that are prone to the PFS formation constitute major areal seal or caprocks in petroleum systems (e.g., Cartwright et al., 2007; Stuevold et al., 2003; Watterson et al., 2000). The formation of PFS can compromise sealing integrity by enhancing permeability, as

<footnote>© 2025. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Polygonal Fault System (PFS) in the Norway Basin. (a) 3D visualization of a portion of the PFS. (b and c) Two horizontal variance (coherency⁻¹) slices at different depths, with their positions indicated in panel (a). The vertical scale in panel (a) is given in milliseconds of Two-Way Travel Time (TWT). Note that the linear scale changes with depth due to variations in seismic wave velocity. For reference, the thickness of the faulted layer (tier) in panel (a) is approximately 1,000 m. The images were generated using seismic data from Gay and Berndt (2007).</figcaption></figure>

suggested by numerous cases where fluid flow indicators (e.g., fluid pipes, mud volcanoes, shown in Figure 1a, pockmarks, or gas hydrate accumulation) are associated with PFS (Gay & Berndt, 2007; Ho et al., 2016, 2018; Hustoft et al., 2007; Zeng et al., 2022). These structures therefore are important for petroleum exploration (Hantschel & Kauerhauf, 2009; Jackson et al., 2014) as well as for CO₂ sequestration and nuclear waste disposal (e.g., Cartwright, 2011).

PFS are rather enigmatic structures. They mostly form on tectonically inactive passive continental margins where regional dips are very small. PFS therefore are not of tectonic origin and are generated in absence of regional horizontal displacements due to the internal processes within a sedimentary sequence developing during its “passive” burial. Local displacements do occur along the faults of PFS, which are all normal faults (Cartwright & Dewhurst, 1998) and hence result from the horizontal extension of the bearing layers. Given the variety of orientations of polygonal faults (PFs) forming PFS, the horizontal extension appears radially isotropic.

The fault throws profiles, $\delta_z$ versus depth $z$, along PFs are C-shaped, with a maximum $\delta_z = \delta_z^{\max}$ in the central part of the faulted tier. $\delta_z$ values approach zero near the top and bottom of the tier (Hu et al., 2024; King et al., 2022; Laurent et al., 2012; Shin et al., 2010), which is typical for simple blind faults of tectonic origin (Barnett et al., 1987; Baudon & Cartwright, 2008; Walsh & Watterson, 1988). $\delta_z^{\max}$ roughly scales with the tier thickness $T$ as $\delta_z^{\max} \sim 0.045T$ (Shin et al., 2010). The vanishing of $\delta_z$ towards the basal fault tip is understandable. As for the upper tip, there are a number of examples where active PFs are very close to the seafloor or even cross it, as in the Hatton Basin (Berndt et al., 2012), the Norway Basin (King & Cartwright, 2020; Laurent et al., 2012), or the Grenada Basin (Gay et al., 2021). This suggests that, at least in these cases, faulting becomes synsedimentary at a certain point, propagating upwards through newly sedimented layers (King & Cartwright, 2020; King et al., 2022; Lonergan et al., 1998; Morgan et al., 2015; Stuevold et al., 2003), with the upper fault tips breaching the contemporaneous depositional surface.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Geometry of PFS. (a) 3D visualization of a portion of the Vøring Basin PFS, Norwegian Sea (from King et al., 2022). (b) Three horizontal coherence slices whose position is shown on the vertical seismic profile in panel (c) (from Cartwright, 2011).</figcaption></figure>

The case where PFs reach the seabed is presented in Figure 2a, which also shows higher fault density at deeper horizons. The fault pattern also varies with depth (Figures 1b, 1c, 2b, and 2c). Moreover, the PFS architecture appears as an imbrication of fault networks of different scales and maturities (Gay et al., 2004; Wrona et al., 2017), likely reflecting the processes responsible for the formation of these fault systems.

Several mechanisms and hypotheses have been proposed to explain the formation of PFSs. Among them: (a) bulk volume loss of pore fluids during burial (Cartwright & Lonergan, 1996; Gay et al., 2004); (b) gravitational sliding influenced by far-field tectonic stresses at a regional scale (Clausen et al., 1999; Ghalayini et al., 2017); (c) density inversion and associated hydrofracturing (Henriet et al., 1991; Watterson et al., 2000). The role of low residual

friction on forming fault planes has also been highlighted (Goult, 2008). These hypotheses have been discounted either because they contradict certain observations or because they are not applicable to many known PFS occurrences (see Goult, 2008; Cartwright, 2011 for a review).

The mechanisms that appear to be the most generic and fundamental are those involving volumetric contraction of sediments, although none of them is completely clear, comprehensive, or fully elaborated yet. Dewhurst et al. (1999) suggested that the contraction or syneresis of colloidal (very fine-grained) clay-rich sediments with high smectite content occurs due to attraction between clay flocs, causing expulsion of the fluid phase. Volumetric contraction can also occur due to a diagenetically induced reduction in grain volume (Cartwright, 2011; Shin et al., 2008, 2010) and porosity reduction caused by the opal-A to opal-CT transformation (Davies & Ireland, 2011; Neagu et al., 2010; Volpi et al., 2003).

The diagenetically induced volumetric contraction of a layer (tier) at fixed boundaries results in a reduction of the horizontal compressive stresses, which at some point become sufficient to meet one or another failure criterion and cause the formation of normal faults. We are aware only of three modeling studies of this process. Tuckwell et al., 2003 conducted 2D simulations of the propagation of opening mode fractures in plane vue under the remote isotropic horizontal tension within the framework of linear elastic fractur mechanics. These models have limited application to PFS as they do not generate normal faults, and do not take it to account inelastic

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Setup of 2-layer 3D models. The upper layer, Layer 1 (active layer), exhibits elastoplastic properties and undergoes volumetric contraction and faulting. The lower layer, Layer 2 (passive substratum layer), is predominantly purely elastic in the presented models, except for some cases where it is elastoplastic. The bottom of this layer is fixed in all directions. The lateral boundaries of the entire model are constrained normal to the boundaries (roller conditions). The model is prestressed to the initial stresses $\sigma_{xx}^{ini} = \sigma_{yy}^{ini} = K_0\sigma_{zz}^{ini}$, and $\sigma_{xx}^{ini} = \sigma_v$, where $K_0 = 0.7$. Vertical (lithostatic) effective stress $\sigma_v$ is also applied to the model surface. In certain models, the layers are separated by frictional-cohesive interfaces (denoted as 1 in the figure).</figcaption></figure>
deformations. In addition, the polygonal falting is essentially 3D phenomenon. Shin et al. (2008) conducted 2D simulations (in cross-section this time) aimed at demonstrating that the volume reduction due to grain dissolution can lead to shear faulting. Roberts et al. (2014) performed 3D finite-element modeling, integrating isotropic chemical compaction resulting from the opal-A to opal-CT transition. They obtained realistic networks of polygonal faults, which lends credibility to their approach. However, these authors did not address the initiation and evolution of the 3D fault systems. They used a complex constitutive framework, including phenomenological three-invariant Cam Clay-type yield and plastic potential functions, as well as temperature- and time-dependent chemical porosity (volume) reduction. This model is defined by several tens of parameters, most of which are very poorly constrained, making the results obtained by adjusting these parameters lack physical clarity and transparency.

We use here a much simpler constitutive framework in 3D finite-difference modeling, which allows for capturing and accurately modeling constitutive instabilities, leading to the theoretically predicted deformation bifurcation and formation of deformation localization bands (Rice, 1976; Rudnicki & Rice, 1975), as well as their regular networks (Chemenda, 2007). These networks correspond to PFS and are modeled here using Mohr-Coulomb-type constitutive framework with a simple volume reduction algorithm. The focus is on the initiation, evolution, and architecture of the fault networks and their dependence on the model parameters.

## 2. Modeling Setup

A two-layer model is fixed in all directions at the bottom and has roller conditions along the lateral boundaries (Figure 3). The model surface can be free or subjected to a prescribed lithostatic pressure $\sigma_v$. Layer 1 or active layer of thickness $T$ represents a tier susceptible to volumetric contraction (e.g., smectite-rich mudstone) and Layer 2 of thickness $T_s$ is composed of a generic mudstone representing substratum. The horizontal sizes $L$ and $D$ of most models are $L = D = 1,600$ m, $T_s = 100$ m, and $T$ varies from 100 to 400 m.

The whole model is prestressed in the vertical ($\sigma_v = \sigma_{zz}$) and horizontal ($\sigma_h = \sigma_{xx} = \sigma_{yy}$) directions with the imposed ratio $K_o = \sigma_h/\sigma_v$ of 0.7, typical for fine-grained sediments (Casey, 2014; Goulty & Swarbrick, 2005). The orientation of stresses $\sigma_{ij}$ is in accordance with axes in Figure 3 ($i, j = 1, 2, 3$ or $x, y, z$). All stresses are considered effective, corresponding to drained conditions with hydrostatic pore pressure. We neglect both viscous and thermal effects. Principal effective stresses $\sigma_i$ are defined following the rock mechanics convention with $\sigma_1 \ge \sigma_2 \ge \sigma_3$, where positive values represent compressive stress and negative values represent tension.

In most models, Layer 2 is purely elastic and can be separated from Layer 1 by interface 1 (Figure 3) with a limited shear strength to allow mutual sliding of the layers. The interface conditions and parameters are summarized in Appendix A2. In some models, Layer 2 is allowed to deform inelastically. Layer 1 is elastoplastic and undergoes a volumetric contraction as described in the next section. This contraction results in the horizontal incremental extension straining and faulting (shear banding) of Layer 1, which is the subject of this study.

## 3. Constitutive Model for Layer 1

This model is formulated within the elastoplastic framework and encompasses both tensile and shear failure mechanisms, coupled in a manner where both tensile ($\sigma_t$) and uniaxial or unconfined compression ($\sigma_c$) strengths decrease simultaneously with accumulated inelastic volume strain $\epsilon^p$. These strengths reach zero at $\epsilon^p = \epsilon^0$, where $\epsilon^0$ characterizes the material's brittleness/ductility.

To reduce the mesh dependence of the modeling results, the mesh size and viscous regularization are integrated into the constitutive formulation, as detailed by Chemenda et al. (2021, 2022) and summarized in Appendix A1.

The diagenetic volumetric contraction of the material is simulated by incremental reduction of volumetric strain $\Delta\tilde{\epsilon}$ in all numerical zones of Layer 1 between calculation steps $n$ and $n + 1$. The elastic trial stress $\sigma_{ij}^{tr}$ is computed as

$$ \sigma_{ij}^{tr} = \sigma_{ij}^{n} + K \Delta \tilde{\epsilon} \delta_{ij} \quad (1) $$
where *K* is the bulk modulus, and $\delta_{ij}$ is the Kronecker delta. If $F(\sigma_{ij}^r, \epsilon^{pn}) < 0$, $\sigma_{ij}^{n+1} = \sigma_{ij}^r$, otherwise the elastoplastic solution $(\sigma_{ij}^{n+1}, \epsilon^{pn+1})$ is computed following classical procedure by projecting the stress state onto the current yield surface (F is the yield function defined by Equation A1).

Accumulated during cycling volume contraction (reduction in volumetric strain) $\gamma = \int \Delta\tilde{e}$ attributed to diagenetic processes corresponds to the porosity reduction, $\Delta\phi$, $\Delta\phi = 100\gamma\%$. The volume contraction can be isotropic as in Equation 1 or not, with, for example, larger vertical contraction. Both elastic and strength parameters can easily be introduced as functions of $\gamma$. This, however, would strongly complicate and increase the uncertainty of the analysis, which may become necessary if the above simple model does not provide realistic results consistent with the available data on PFS.

Parameter $\gamma$ can be viewed as geological time, or more exactly, is related to this time. The rate of $\gamma$ increase under natural conditions can vary with time and depth depending on the temperature, pressure, and possibly other parameters/factors as well. Analysis of these relations is beyond the scope of this paper.

The presented constitutive model has been implemented into the finite-difference dynamic time-marching explicit code FLAC3D (Itasca, 2024), employed in a quasistatic mode (i.e., a steady-state solution of a fully dynamic problem is obtained by a dynamic relaxation method).

## 4. Choice of Model Parameters

The mechanical properties of clayey marine sediments vary considerably with the burial depth, depending on their mineralogical composition, sedimentation rate, thermal gradient, diverse early diagenetic processes, as well as the applied approaches for determining these properties. These approaches can be theoretical (e.g., Cosenza et al., 2023; Ortega et al., 2007) but mostly they are experimental with employment of different techniques (e.g., Carrier et al., 2016; Daigle et al., 2020; Hamilton, 1971, 1976; Ikari et al., 2015, 2018; Minaeian et al., 2017; Mondol et al., 2007; Rabinowitz et al., 2018; Sanada et al., 2009; Stipp et al., 2013; Wang et al., 2001).

The parameter ranges become much narrower when considering only brine-saturated, clay-rich sediments with a high smectite content at burial depths $z$ from the seabed to several hundred meters that are prone to polygonal faulting. In this depth range (from ~0 to ~500 m), the porosity $\phi$ decreases from ~60%–80% at $z \approx 0$ m to 30%–50% at $z \approx 500$ m (Mondol et al., 2007 and references therein), and the density $\rho$ increases from ~1,000 to 1,800 kg/m$^3$ (Hamilton, 1976; Tenzer & Gladkikh, 2014). The material stiffness and strength increase as well and are much smaller than the usual values for completely lithified rocks. The Young modulus $E$ typically increases from very small values to a few GPa and Poisson's ratio $\nu$ is more than ~0.47, approaching 0.5 for large $\phi$ (Mondol et al., 2007; Mondol, Bjørlykke, & Jahren, 2008; Mondol, Jahren, et al., 2008; Vanorio et al., 2003). According to Hamilton (1971) these values are typical for a wider class of clayey sediments. This author reports very small $E$ values, $E < 0.3$ GPa.

In most studies, including those cited above, the elastic moduli were calculated from *P* and *S* velocities ($V_p$ and $V_s$, correspondingly) of ultrasonic waves traveling through specimens in the experimental tests and therefore are dynamic. It is well known that $V_p$ increases and $V_s$ decreases with the increase in water saturation level (King, 1966; Mondol et al., 2007; White & Sengbush, 1953 and references herein), leading to the increase in the $V_p/V_s$ ratio, which defines the dynamic $\nu$ value bringing it closer to 0.5 (e.g., Hamilton, 1971), with increase of porosity, saturation level, clay and smectite content. Cyclic loading/unloading in such dynamic tests is too fast to allow pore fluid to flow, so their results correspond to undrained conditions regardless of the experimental setup. Under these conditions the compressibility of the sediments is largely defined by practically incompressible water, which corresponds to $\nu \approx 0.5$.

The deformation leading to PFS formation is very slow (static) and occurs rather under drained conditions assumed in previous (King et al., 2022; Roberts et al., 2014; Shin et al., 2008) and present numerical studies. Unfortunately, we are not aware of the experimental tests defining the elastic moduli of the sediments in question under static drained conditions. We can only suppose that static $\nu$ under drained conditions is smaller than its dynamic counterpart. Indeed, the axial static compression of a specimen under drained conditions will be accompanied by the fluid escape from the pores, which leads to a reduction of its lateral expansion. This corresponds to a larger compressibility and hence a lower $\nu$. It is interesting that the dynamic $\nu$ of saturated sediments, which is close to 0.5 as indicated above, can be as small as 0.1 for dry sediments (e.g., Mondol et al., 2007).
The actual $\nu$ value characterizing slow deformation of soft, high-porosity saturated sediments should lie between these extreme values. For our modeling, we chose the reference values $E = 2$ GPa and $\nu = 0.3$, although we vary these parameters to assess their impact on the results.

The cohesion $c$ (related to $\sigma_c$ by condition Equation A4) of shallow sediments is very low (a few kPa) and increases with depth (Daigle et al., 2020; Ikari et al., 2015; Stipp et al., 2013). Measurements in static tests using a direct-shear apparatus have shown that the gradient of $c$ increase is about 30 kPa per 1 MPa of effective overburden stress $\sigma_v$ (Ikari & Kopf, 2011). This yields around 0.1 MPa for $c$ at depths of several hundred meters. This value was assumed for modeling.

The internal friction angle $\varphi$ (related to the friction parameter $\alpha_c$ of the constitutive model by condition Equation A3) of marine clayey sediments does not show clear $z$-dependence in the subsurface conditions of interest, and ranges from 6° to 46° (Ikari et al., 2018 and references therein; Okuda et al., 2021). These values are consistent with the internal friction measured on clay-bearing fault gouges (Morrow et al., 2017 and reference therein). However, considering the effect of fluid-saturation and high smectite content in sediments of interest, their internal friction angle $\varphi$ probably does not exceed 20° (Ikari et al., 2015; Ikari & Kopf, 2011; Saffer et al., 2001). The estimates of residual internal friction angle $\varphi^{\text{res}}$ vary from 3.4° to 10.3° (Goulty, 2008; King et al., 2022; Mesri & Shahien, 2003 and references therein). We assumed $\varphi = 19.5^\circ$ ($\alpha_c = 2$) and $\varphi^{\text{res}} = 5.2^\circ$ ($\alpha_c^{\text{res}} = 1.2$).

The tensile strength $\sigma_t$ ranges from ~20 kPa near the seafloor to ~700 kPa at $z \approx 500$ m (Daigle et al., 2020), with the average value over the depth of around 350 kPa assumed in our models.

The dilatancy factor $\beta$ (related to the coefficient $\beta_c$ by Equation A7) is set to a small initial value of 0.05 ($\beta_c = 1.1$) to ensure the material softening. $\beta$ reduces as $\varepsilon^p$ increases and stabilizes at $\beta = 0$ at $\varepsilon^p = \varepsilon^0$, that is, at the end of softening as defined by Equation A6.

The constitutive model for Layer 1 includes three more regularization parameters, $C$, $\tilde{\varepsilon}^0$, and $\eta$ (Appendix A3). $C$ and $\tilde{\varepsilon}^0$ define energetic and $\eta$ viscous regularizations; they are set to $C = 1$ m, $\tilde{\varepsilon}^0 = 0.012$, $\eta = 0.1$ (Chemenda et al., 2021, 2022). $\tilde{\varepsilon}^0$, related to $\varepsilon^0$ by Equation A11, also defines the brittleness/ductility of the material and is varied in the models.

In our simplified setup in Figure 3, the layer undergoing faulting (Layer 1) is homogeneous with uniform properties representative of the entire faulted tier. The thickness $T$ of this layer is constant, whereas in nature, it may increase during syn-sedimentation faulting due to additional deposits or decrease due to compaction. The representative $T$ value could be of 200 m, but we test $T$ values in the range 100–400 m.

Since Layer 1 is intended to represent a thicker natural tier, it does not necessarily crop out at the seafloor. It can be imagined to be located somewhere in the middle part of the tier and hence should experience the overburden pressure $\sigma_v$ (Figure 3). We assume that the top of this layer is located at depth $z_0 = 200$ m and hence undergoes effective overburden pressure $\sigma_v = (\rho - \rho_w)gz_0$ of 1.6 MPa ($\rho_w = 1,000$ kg/m$^3$ is the water density, and $g$ is the acceleration due to gravity).

Thus, the assumed reference parameters of Layer 1 are: $E = 2$ GPa, $\nu = 0.3$, $c = 0.1$ MPa ($\sigma_c = 0.28$ MPa), $\varphi \approx 19.5^\circ$ ($\alpha_c = 2$), and $\varphi^{\text{res}} = 5.2^\circ$ ($\alpha_c^{\text{res}} = 1.2$), $\beta = 0.05$ ($\beta_c = 1.1$), $C = 1$ m, $\tilde{\varepsilon}^0 = 0.012$, $\eta = 0.1$, and $\rho = 1,800$ kg/m$^3$.

For the substratum Layer 2 (Figure 3) of most models, we assume pure elastic behavior with the reference elastic moduli and density of Layer 1 indicated above. Thickness $T_s$ of Layer 2 in all models is 100 m.

All presented models are discretized with structured brick-shaped (cubic) numerical zones with the zone size $\Delta z$ of 8 m.

## 5. Modeling Results

In this section we present five series (Series 1 to Series 5) of numerical models run with parameters given in Table 1. In Series 1 to 4, we vary one parameter characterizing Layer 1 or coupling between the layers, with Layer 2 being purely elastic. In Series 5, Layer 2 is allowed to deform inelastically and has different strengths in different models. One or several evolutionary stages of each model are presented for different accumulated volumetric contraction strain $\gamma$. Only active layers undergoing faulting (Layer 1, Figure 3) are shown. Faults
<table><caption>**Table 1**
Model Parameters</caption><thead><tr><th rowspan="2">Model</th><th colspan="3">Layer 1</th></tr><tr><th>T (m)</th><th>ε<sup>0</sup></th><th>E (GPa)</th></tr></thead><tbody><tr><td>1, Figure 3</td><td>200</td><td>0.012</td><td>2</td></tr><tr><td>2, Figure 4a</td><td>100</td><td>0.012</td><td>2</td></tr><tr><td>3, Figure 4c</td><td>300</td><td>0.012</td><td>2</td></tr><tr><td>4, Figure 4d</td><td>400</td><td>0.012</td><td>2</td></tr><tr><td>5, Figure 5b</td><td>100</td><td>0.005</td><td>2</td></tr><tr><td>6, Figure 5c</td><td>100</td><td>0.003</td><td>2</td></tr><tr><td>7, Figure 6a</td><td>100</td><td>0.003</td><td>1</td></tr><tr><td>8, Figure 6c</td><td>100</td><td>0.003</td><td>4</td></tr></tbody></table><table><thead><tr><th rowspan="2"></th><th colspan="2">Interface strength</th></tr><tr><th>τ<sub>0</sub><sup>int</sup> (MPa)</th><th></th></tr></thead><tbody><tr><td>9, Figure 7b</td><td>100</td><td>0.003</td><td>2</td><td>0.15</td></tr><tr><td>10, Figure 7c</td><td>100</td><td>0.003</td><td>2</td><td>0.07</td></tr></tbody></table><table><thead><tr><th rowspan="2"></th><th colspan="3">Plastic Layer 2</th></tr><tr><th>c (MPa)</th><th>φ°</th><th></th></tr></thead><tbody><tr><td>11, Figure 8a</td><td>100</td><td>0.003</td><td>2</td><td>0.28</td><td>5.2</td></tr><tr><td>12, Figure 8b</td><td>100</td><td>0.003</td><td>2</td><td>0.28</td><td>0</td></tr><tr><td>13, Figure 8c</td><td>100</td><td>0.003</td><td>2</td><td>0.14</td><td>0</td></tr></tbody></table>

(deformation bands) are depicted in the presented contour plots in terms of the uniaxial compression strength $\sigma_c$, decreasing during inelastic deformation.

Indicated in the table are only parameter values that vary in different models. Other parameters have reference values. Reference values for Layer 1 are: $E = 2$ GPa, $\nu = 0.3$, $c = 0.1$ MPa ($\sigma_c = 0.28$ MPa), $\varphi = 19.5°$ ($\alpha_c = 2$), and $\varphi^{\text{res}} = 5.2°$ ($\alpha_c^{\text{res}} = 1.2$), $\beta = 0.05$ ($\beta_c = 1.1$), $C = 1$ m, $\bar{\varepsilon}^0 = 0.012$, $\eta = 0.1$, and $\rho = 1,800$ kg/m$^3$. Reference parameter values for Layer 2 are: $E = 2$ GPa, $\nu = 0.3$, $\rho = 1,800$ kg/m$^3$, and $T_s = 100$ m. For the elasto-plastic Layer 2 in Models 11 to 13, the additional reference parameters are: $c = 0.1$ MPa, $\varphi = 19.5°$, $\varphi^{\text{res}} = 5.2°$, $\beta = 0$, $C = 1$ m, $\bar{\varepsilon}^0 = 1$, and $\eta = 0.1$. Models 9 and 10 were run with interfaces between the layers. The interface shear strength $\tau_0^{\text{int}}$ is given in the table. The horizontal sizes of the models are 1,600 m.

**Series 1:** Variation of the Layer 1 thickness $T$ (Figures 4 and 5). Figure 4 shows evolution of Model 1 with reference parameters and $T = 200$ m from the very early stage of instabilities (deformation localization) at $\gamma \approx 0.001$ (Figure 4a) to a mature fault network at $\gamma$ up to 0.05 (Figure 4i). During the initial stages of faulting and fault network structuring, one can observe the reduction of the fault (deformation band) density, with some initial faults ceasing activity (Figures 4a–4d). At $\gamma \sim 0.005$ (Figure 4e), progressive material damage results in its complete failure (reduction of $\sigma_c$ to zero) along certain segments of deformation bands, although displacement along them continues as long as $\gamma$ increases. At larger $\gamma$, the fault pattern stabilizes, although the initiation and evolution of new faults can still be observed (Figures 4f–4i).

In Figure 5, in addition to the Model 1 with $T = 200$ m (Figure 5b), we present three more models: Models 2 (a), 3 (c), and 4 (d) with $T = 100, 300$, and $400$ m, respectively, deformed to $\gamma = 0.01$. One can see that the $T$ increase leads not only to an increase in fault spacing $S$, but also in the brittleness of faulting (or strain localization) process. At $T = 100$ m, the minimum $\sigma_c$ value within the faults does not reduce below 0.1 MPa, while at $T = 300$ and $400$ m, $\sigma_c$ already reached zero along certain fault segments. For $T = 400$ m, length/number of these segments, as well as the displacement or inelastic shear deformation along them, are larger than for $T = 300$ m. The plots in Figures 5e and 5f show that not only $S$ scales with $T$, but the fault throw $\delta_z$ and its maximum $\delta_z^{\text{max}}$ as well. The latter increases with $\gamma$ almost linearly (Figure 5f) as suggested by correlation coefficients $R^2$ (all greater than ~0.98) corresponding to the linear approximating functions given in the figure caption.

**Series 2:** Variation of the brittleness/ductility, controlled by parameter $\bar{\varepsilon}^0$ (Figure 6). This parameter is different in the three models (Models 2, 5, and 6) of this series. The other parameters are the reference ones with $T = 100$ m. The results in Figure 6 show that the brittleness of the model response predictably increases as $\bar{\varepsilon}^0$ decreases. This is reflected notably by the acceleration of the material failure within the faults and the earlier structuring of the final fault pattern at smaller $\bar{\varepsilon}^0$.

**Series 3:** Variation of the Young's modulus $E$ (Figure 7). This parameter differs across the three models (Models 7, 6, and 8) in this series. The other parameters are set to their reference values except for $\bar{\varepsilon}^0$, which is now 0.003. This parameter has been reduced in this and the following models with $T = 100$ m to decrease the calculation time. Figure 7 shows that an increase in $E$ has a similar effect on the deformation style as an increase in $\bar{\varepsilon}^0$ in the previous Series, that is, increasing material stiffness results in a more ductile (less brittle) response.

Unlike $E$, the variation of the Poisson's ratio $\nu$ in the range 0.2–0.48 has little effect on the faulting.

**Series 4:** Models 6, 9, and 10 (Figure 8) with reference parameters and a frictional-cohesive interface between the active layer (Layer 1) and substratum layer (Layer 2). The shear strength $\tau_0^{\text{int}}$ of the interface differs among these models and is provided in the caption of Figure 8 and Table 1. One can see that the decrease in $\tau_0^{\text{int}}$ (or in the coupling between the layers) leads to an increase in the fault spacing. The fault throws increase as well.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Model 1: Contour plots of the unconfined compression strength $\sigma_c$ showing the initiation and evolution of the polygonal fault system in Layer 1 for different volumetric contraction $\gamma$ indicated directly on the images of successive deformation stages from panels (a) to (i). The model has reference parameters with $T = 200$ m.</figcaption></figure>

Series 5: Layer 2 in these models (Models 11 to 13) is elastoplastic with different strength. There is no dilatancy ($\beta = 0$), and therefore, no strain softening occurs in this layer. Layer 1 has thickness $T = 100$ m and reference parameters except for $\bar{\epsilon}^0$, which is 0.003. In the models where strength parameters of Layer 2 were the same as those of Layer 1 (i.e., the unconfined compression strength $\sigma_c = 0.28$ MPa, and internal friction angle $\varphi = 19.5^\circ$), Layer 2 practically did not experience inelastic deformation. Therefore, the result was the same as that of Model 6 in Figure 6c with purely elastic Layer 2.

In Model 11 (Figure 9a), we reduced the frictional strength of Layer 2 by decreasing its internal friction angle from $19.5^\circ$ to $5.2^\circ$. This caused some inelastic straining in the uppermost part of Layer 2 but did not significantly alter the faulting process or fault pattern in Layer 1 (compare Figure 9a and the last deformation stage of Model 6 in Figure 6c). In Model 12 (Figure 9b), $\varphi$ is reduced to zero, and in Model 13 (Figure 9c), we additionally reduced the cohesive strength of Layer 2 by half, decreasing $\sigma_c$ from 0.28 to 0.14 MPa. This caused a progressive increase in the inelastic straining of Layer 2 and fault spacing in Layer 1 (Figure 9). Fault throws also increased. The results of this series (Series 5) are similar to those of the previous series in Figure 8, suggesting that reducing shear tractions at the base of a faulted layer leads to an increase in the fault spacing S.

Finally, we present one more model (Model, 14) with reference parameters of Layer 1, purely elastic Layer 2, and large sizes, $T = 400$ m, $L = D = 3.2$ km. This model with $\sim 2.5 \times 10^6$ numerical zones is used to present details of the fault network architecture (Figure 10), showing notably the network complication and densification toward the bottom of the active layer.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) to (d) Four models with reference parameters and different thickness *T* of Layer 1 at *γ* = 0.01. (a) Model 2, *T* = 100 m; (b) Model 1 (the same as in Figure 4), *T* = 200 m; (c) Model 3, *T* = 300 m; and (d) Model 4, *T* = 400 m. (e) Variation of fault throw *δ*<sub>z</sub> with depth *z* for different *T* at *γ* = 0.01, 1, 2, 3, and 4 in panels (e and f) correspond to Models 1 to 4, respectively. (f) Maximum throw *δ*<sub>z</sub><sup>max</sup> versus volumetric contraction *γ* in the models. Solid lines are linear functions approximating points for different *T* with correlation coefficients *R*<sup>2</sup>: for *T* = 100 m, *δ*<sub>z</sub><sup>max</sup> = 30*γ* − 0.21 and *R*<sup>2</sup> = 0.979; for *T* = 200 m, *δ*<sub>z</sub><sup>max</sup> = 140*γ* − 0.776 and *R*<sup>2</sup> = 0.989; for *T* = 300 m, *δ*<sub>z</sub><sup>max</sup> = 181.3*γ* − 0.083 and *R*<sup>2</sup> = 0.997; for *T* = 400 m, *δ*<sub>z</sub><sup>max</sup> = 260.2*γ* − 0.257 and *R*<sup>2</sup> = 0.998. *δ*<sub>z</sub> at different *z* was evaluated from the maximum (*u*<sub>max</sub>) and minimum (*u*<sub>min</sub>) *z*-displacements in the horizontal sections of Layer 1 as *δ*<sub>z</sub> = (*u*<sub>max</sub> − *u*<sub>min</sub>)/2.</figcaption></figure>

**6. Discussion**

Numerous physical, chemical, and biogenic processes define the depth- and time-dependent evolution of prop-
erties and deformations of sedimentary material deposited at seafloor and buried to depth. We undertook
modeling characterized by extreme conceptual simplicity, considering two layers with constant and uniform
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Models 2 (a), 5 (b), and 6 (c) with *T* = 100 m and reference parameters except for ε⁰, which differs in these models. This parameter is 0.012, 0.005, and 0.003 in panels (a), (b), and (c), respectively. Three evolutionary stages of each model are shown in three columns for the γ values indicated in the figure on the left of each row.</figcaption></figure>

properties over time. The only process driving the material failure and faulting of the upper layer is its isotropic contraction or volume reduction, which is supposed to be caused by diagenesis in clay-rich, water-saturated sediments. Despite this simplicity, the 3D networks (systems) of polygonal faults generated in the models are very similar to natural PFS (compare, e.g., Figures 2, 10, and 11). This similarity is not limited to the first-order geometrical attributes of the inclined normal faults that form polygonal systems in horizontal cross-sections. The change of the natural fault density and pattern with depth (Figures 1 and 2) is captured by the models as well (Figure 10). Closer to the bottom of the active (faulted) layer, the faults are not localized but rather distributed in arrays having rounded shapes (Figure 10c). Such details are generally difficult to resolve by seismic data (limit of resolution), but they can be guessed in Figure 2b III and are clearly visible in exhumed and outcropping ancient PFS (Figure 11b). Thus, it appears that the models capture the principal features of the PFS architecture down to certain details. This may suggest that they also correctly reproduce the first-order process of the PFS formation.

## 6.1. Origin of PFS in the Models

What is the mechanism of PFS initiation and evolution in the models? Since the upper or active layer and the entire model are fixed along the lateral boundaries normal to them (roller conditions, Figure 3), any tectonic effects or influence are excluded. The horizontal contraction of this layer (which is a part of the applied isotropic contraction) generates horizontal incremental extension strains and reduces the horizontal compressive stresses $\sigma_h = \sigma_{xx} = \sigma_{yy} = \sigma_2 = \sigma_3$ within the layer. Since this layer is attached to the lower (substratum) passive layer, the latter experiences increasing horizontal compression during the contraction of the upper layer. The coupling
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Models 7 (a), 6 (b), and 8 (c) with $T = 100$ m and reference parameters except for $\bar{\epsilon}^0$ and $E$. Here, $\bar{\epsilon}^0 = 0.003$, and $E$ varies among the models: it is 1, 2, and 4 GPa in Models 7, 6, and 8, respectively. Three evolutionary stages of each model are shown for the $\gamma$ values indicated in the figure. Model 6 (Figure 6c) is included here for the comparison.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Models 6 (a), 9 (b), and 10 (c) with $T = 100$ m and reference parameters for Layer 1, except for $\bar{\epsilon}^0$, which is 0.003. In Models 9 (b) and 10 (c), Layers 1 and 2 (Figure 3) are separated by the interfaces with imposed shear strengths $\tau_{0i}^{min}$, which are 0.15 MPa in panels (b) and 0.07 MPa in panel (c). In model 6 (a) shown here for comparison, the layers are fully coupled, corresponding to infinite interface strength. One deformation stage is shown for each model at $\gamma = 5.5 \times 10^{-3}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Models 11 (a), 12 (b), and 13 (c) with $T = 100$ m and reference parameters for Layer 1 except for $\bar{\epsilon}$, which is 0.003. One deformation stage of each model is shown for $\gamma = 5.5 \times 10^{-3}$. Layer 2 in these models is elastoplastic with the reference parameters for Layer 1 except for internal friction angle $\phi^{\text{lay2}}$ and unconfined compression strength $\sigma_c^{\text{lay2}}$. In panel (a), $\phi^{\text{lay2}} = \phi_{\text{res}} = 5.2^{\circ}$, $\sigma_c^{\text{lay2}} = \sigma_c = 0.28$ MPa; in (b), $\phi^{\text{lay2}} = 0$, $\sigma_c^{\text{lay2}} = \sigma_c$; and in (c), $\phi^{\text{lay2}} = 0$, $\sigma_c^{\text{lay2}} = 0.5\sigma_c$.</figcaption></figure>

between the contracting and passive layers defines the vertical gradient of the horizontal stress-strain evolution and, as a consequence, the change in the fault pattern with depth seen in Figure 10.

The inelastic straining or damage of the active layer starts when the differential stress $\sigma_v - \sigma_h$ reaches a sufficiently large value to meet the yield condition $F = 0$, where $F$ is defined by Equation A1. In our models with reference parameters, this occurs at total volumetric strain reduction $\gamma$ of $3 \times 10^{-4}$, that is, very soon after the contraction onset. This value decreases as the elastic compressibility of the active layer increases. A reduction in the initial horizontal stress or in $K_0$ also leads to the reduction in this $\gamma$ value.

The damage results in the reduction of the material strength (strain softening), which can occur at a different rate defined by the dimensionless hardening modulus $h$

$$h = \frac{1}{G} \frac{\partial F(\sigma_{ij}, \epsilon^p)}{\partial \epsilon^p} \qquad (2)$$

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Contour plots of cumulated inelastic volume strain $\epsilon^p$ in active layer (Layer 1) of the model with reference parameters and sizes $T = 400$ m and $L = D = 3,200$ m at $\gamma = 0.0045$. (a) General 3D view. (b) Vertical $x$-parallel sections (V1 to V5), see Figure 3 for axes orientations. (c) Horizontal sections (H1 to H4). The positions of the sections are indicated in panel (a).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Fault patterns in nature and the model. (a) Horizontal section H4 from the model in Figure 10 at depth of 360 m (T = 400 m). (b) Traces of polygonal faults in chalk mapped using satellite imagery; Khoman Formation, Western Desert, Egypt (from Tewksbury et al., 2014).</figcaption></figure>

where G is the shear modulus, and $\epsilon^p$ is the accumulated inelastic volume deformation, the parameter tracking the history of inelastic deformation. Submitting here the expression for $F(\sigma_{ij}, \epsilon^p)$ from Equation A1, obtain

$$h = -2(1 - \nu) \frac{\sigma_c^0 + \sigma_3(\alpha_c^0 - \alpha_c^{\text{res}})}{E\epsilon^0} \quad (3)$$

where $\alpha_c^0$ is the initial (at $\epsilon^p = 0$) value of the internal friction parameter $\alpha_c$ defined in Equation A3, $\alpha_c^{\text{res}}$ is the residual (at $\epsilon^p \ge \epsilon^0$) value of this parameter, and $\sigma_c^0$ is the initial (at $\epsilon^p = 0$) value of uniaxial compression strength, related to the cohesion c by Equation A4. During the elastoplastic loading of sediments, the initially positive h decreases, and when it becomes as small as the critical modulus $h_{cr}$, the continuation of a uniform deformation becomes impossible and it localizes within narrow deformation localization bands as a result of deformation bifurcation (Rice, 1976; Rudnicki & Rice, 1975). Critical modulus $h_{cr}$ usually have small negative values and depends on the constitutive and stress state parameters. It has different forms for different constitutive models, for the Drucker-Prager model (Rudnicki & Rice, 1975), the Mohr-Coulomb model (Ottosen & Runden, 1991) and a three-invariant model (Chemenda & Mas, 2016). The deformation localization or bifurcation results from a constitutive or material instability whose physical nature lies in the nonlinearity of the material properties contained in the constitutive description. Therefore, the presence of any sort of heterogeneities or weak zones is not necessary for the deformation localization to occur. There are no heterogeneities or fault seeds in the models here. Such heterogeneities can exist in natural materials and, if strong enough, can affect the locations of deformation bands, whose initial thickness in granular materials is several grain sizes (Desrues, 2003; Fortin et al., 2009; Sulem & Ouffroukh, 2006). With further deformation and displacement along the bands, their thickness can either decrease, resulting in clear-cut fractures or faults, or increase, producing thick homogeneous or clustered bands (e.g., Mair et al., 2000; Tran et al., 2019). The former case corresponds to dilatant strain-softening, and the latter to more ductile deformation regimes.

The smaller (the more negative) the $h < h_{cr}$ value (the smaller the $\epsilon^0$ value in our models, Equation 3), the more rapid the strain localization and hence the higher brittleness of the material response. On the contrary, the closer the $h < h_{cr}$ value to $h_{cr}$, the more ductile the response (faulting) is (Chemenda, 2007). At $h > h_{cr}$, the deformation localization does not occur at all, and the deformation remains uniform.

Figure 6 shows how the process of deformation localization changes with the $\epsilon^0$ (hence h) variation. In a more brittle material (Figure 6c), localization occurs more rapidly with larger spacing S between the faults. Then, S
decreases during further deformation. In a more ductile material (Figure 6a), the initiating fault (deformation band) network is very dense, and the material failure within the bands (the reduction in $\sigma_c$) is much slower.

One observes a similar change in the faulting process when increasing Young's modulus $E$ at constant $\epsilon^0$ (Figure 7). This is rather counterintuitive because one would expect a more brittle response of a stiffer material. The explanation follows from Equation 3, showing that $h$, defining the brittleness/ductility of the material response, depends on $\epsilon^0$ and $E$ in the same way. Therefore, to have the same (similar) response of the material with larger $E$, one has to reduce $\epsilon^0$ proportionally. Obviously, there is no total similarity between these two cases, as the indicated parameters and notably $E$, affect not only $h$ but also the stress state of the model in general.

## 6.2. Factors Controlling Fault Spacing

Figure 5 shows that the first-order fault spacing $S$ scales with the layer thickness $T$, as is the case for opening-mode fractures in sedimentary layers. It remains, however, to be studied what the exact form of the $S(T)$ relationship is for PFS and whether it changes with $T$ increase from linear to strongly nonlinear as is the case for fractures (Chemenda, 2022 and references herein). Future investigations will show whether the physics of $S(T)$ relation in both cases is the same or not.

Models in Figure 8 reveal that the fault spacing $S$ depends also on the coupling between active and substratum layers or on the shear strength $\tau_0^{\text{int}}$ of the interface between these layers defining the possibility and magnitude of sliding of one layer past another. A reduction in $\tau_0^{\text{int}}$ results in an increase in $S$. When $\tau_0^{\text{int}}$ is sufficiently large and sliding is impossible (Figure 8a), the shear tractions at the layer interface can only be limited by the strengths of the layers, or only of the active layer, when the substratum or passive layer is purely elastic. If the substratum layer can behave elasto-plastically and is weak enough, the shear stress at the layer interfaces will be limited by the strength of this layer, which is the case in Models 11 to 13 (Figure 9). One can see in this figure that the weaker the substratum layer, the larger the fault spacing. The reduction of the strength of the substratum layer thus has similar effect on the fault spacing as the reduction of $\tau_0^{\text{int}}$ in Figure 8, although these two situations are not equivalent.

It would be challenging to establish a relationship between the first-order fault spacing $S$ and the layer interface and/or substratum layer strengths in natural PFS systems, simply because these strengths are not known. On the other hand, $S$ does seem to scale with $T$, although, to our knowledge, there are no studies focused on this relationship so far.

## 6.3. Fault Throw

Fault throw $\delta_z$ increases from the bottom to top of the active layer. This increase is fast in the vicinity of layer interfaces, representing a sort of boundary layer in fluid mechanics. The thickness $\Delta$ of this layer seems to be weakly dependent on $T$ and can be evaluated from Figure 5e to be several tens of meters. The $\delta_z$ increase slows down towards the top of the layer, where $\delta_z$ reaches its maximum (Figure 5e). The variation of $\delta_z$ with depth $z$ is thus nonlinear. The larger the $T$ value, the stronger this nonlinearity, which is explained by the reduction of the influence of the shear tractions generated at the layer interface with the increase in the distance from this interface. This distance is larger for larger $T$. In the absence of these tractions, $\delta_z$ is constant through the entire layer thickness.

The maximum value of $\delta_z$, $\delta_z^{\text{max}}$, increases practically linearly with $\gamma$ (see Figure 5f and its caption, which includes approximating linear functions and correlation coefficients $R^2$, all exceeding ~0.98). The rate of this increase (the slope $\alpha$ of the $\delta_z^{\text{max}}(\gamma)$ curves) scales with $T$, with $\alpha(T)$ relation being approximated by a linear function with $R^2 = 0.97$. One can hypothesize that $R^2$ will increase, approaching unity with increasing $T$, while $\Delta$ becomes negligible compared to $T$. It follows, that the $\Omega = \delta_z^{\text{max}}/T$ ratio in the models is a function only of $\gamma$, which seems to be the case for natural PFS as well. For these systems, the average value of $\Omega$ derived from seismic data is 0.045 (Shin et al., 2010). Using results from Figure 5 and functions approximating numerical points in Figure 5f (given in the figure caption) obtain, that $\Omega = 0.045$ is attained at $\gamma$ around 0.07.

In other words, to generate natural PFS, the porosity reduction $\Delta\phi$ caused by the diagenetically induced volumetric contraction of sediments should be around 7%. This value is comparable to the porosity reduction during the Opal A-CT transition, where $\Delta\phi$ reaches 15% (Ishii et al., 2011; Rice et al., 1995), and the smectite-illite transition, where $\Delta\phi$ is about 10% (Peltonen et al., 2009). Although the latter typically occurs at temperatures
greater than 100°C, the exact temperature threshold depends on the duration of heating, solution chemistry, liquid/solid ratio, pressure, and interlayer cation composition (e.g., Cuadros & Linares, 1996; Eberl et al., 1993).

Note that an increase in the residual internal friction angle φ^res and/or material ductility results in a higher γ required to achieve the same Ω value.

The $\delta_z^{\max}$ value depends not only on $\gamma$ but on a number of other factors that affect fault spacing $S$. For a given $\gamma$ value, $\delta_z^{\max}$ is larger for larger $S$, which increases with the increase in $T$ and the reduction in shear coupling or shear stresses $\tau_0^{\text{int}}$ between the active and substratum layers. $S$ also increases with the reduction of the negative hardening modulus $h$ (increase in $|h|$) that occurs with the increase in the brittleness (or $\varepsilon^0$ value) of the active layer (Figure 6), or reduction of the Young's modulus (Figure 7) as follows from Equation 3. Note that the increase in $\alpha_c^0 - \alpha_c^{\text{res}}$ (e.g., the reduction in $\alpha_c^{\text{res}}$ at constant $\alpha_c^0$), as well as the increase in $\sigma_c^0$, also lead to the increase in the magnitude of negative $h$, hence a more brittle response. This has been confirmed by the modeling.

The fact that $\delta_z^{\max}$ scales with $S$ at a given $\gamma$ has simple explanation. Imagine the active layer of the model in Figure 3 is released along the lateral borders and bottom after the uniform elastic volumetric contraction $\gamma$. This would result in the layer shrinkage and reduction of its horizontal size from $L$ to $L - \Delta L$, where $\Delta L = \gamma L/3$. For $\gamma = 0.07$ and $L = 1,600$ m assumed in the models, $\Delta L = 37$ m. Now, if instead of releasing lateral boundaries of the layer, we cut it in the middle along an inclined plane, then the two segments of the layer will slide past one another along this plane (fault) during the layer contraction. The horizontal component of the normal-sense displacement along this fault will be $\Delta L$. If the fault dip angle is 45°, the vertical component of the displacement along the fault (throw $\delta_z$) will also be $\Delta L$. The equivalent $S$ value in this case is $S = L/2$. If instead of one, $n$ faults are created in the x or y direction (or both), then $S = L/(n + 1)$ and the throw along each fault will be $\Delta L/n$.

Understanding physical mechanism defining relation between $S$, $T$, and $\tau_0^{\text{int}}$ requires a separate investigation, although it should not be very different from that for opening-mode fractures characterized by similar relations between the indicated parameters (Chemenda, 2022).

Another issue to be addressed in future studies is the form of the relation $\delta_z(z)$, which in natural PFS exhibits a C-shaped profile with a maximum at $z$, approximately 2/3 of the thickness of the faulted tier (Gay et al., 2021; King et al., 2022; Laurent et al., 2012; Turrini et al., 2017; Wrona et al., 2017). Near the bottom and top of the tier, $\delta_z \approx 0$. In the models, $\delta_z$ is also near zero at the active layer bottom but increases up to the surface where $\delta_z$ reaches its maximum (Figure 5e). The models thus need further developments by integrating new elements for more complete reproducing the natural process. One possibility is introducing a z-dependence of the volumetric contraction rate, which is high up to a certain depth and then reduces to zero toward the layer surface, as there are no reasons to assume that this rate is constant throughout the entire tier. However, the most important factor responsible for the observed shape of $\delta_z(z)$ relation appears to be a syn-faulting sedimentation, leading to a continuous involvement of new (fresh) material layers into contraction and faulting (e.g., Gay & Berndt, 2007; Neagu et al., 2010).

## 7. Conclusions

The 3D architecture (pattern) of polygonal faults generated here in the numerical models has remarkable similarity with the natural PFS down to certain details (compare Figures 2, 10, and 11), suggesting that the origin and formation mechanism of fault formation in both cases are similar as well.

The force driving faulting arises from the diagenetically induced volumetric contraction of the material within a layer attached (completely or partly) to the substratum, which results in the generation of incremental horizontal extension strains, the reduction of horizontal compressive stresses, and consequently, an increase in deviatoric stresses within this layer.

This leads to the initiation of elastoplastic yielding and strain-softening of the material, which is an unstable process, and at a certain stage, results in deformation bifurcation and strain localization within narrow bands that are the incipient faults. As was previously predicted theoretically, in the material (layer in our case) with uniform properties subjected to a uniform straining, this instability results in the formation of a regular network of deformation bands. The PFSs appear to represent a perfect 3D natural example confirming this prediction, as the
band (fault) network develops within a sedimentary layer (tier) that was deposited in a calm environment and therefore possesses uniform properties over vast areas. The extension straining, internally generated within the layer due to its volumetric contraction, is also uniform (is the same along a horizon) over large areas.

The fault spacing S in the models scales with the layer thickness T and increases with the reduction of shear tractions $\tau_0^{int}$ at the interface between active and substratum layers. The relation between these parameters is similar to that for opening-mode fractures, although its exact form (strongly nonlinear for opening-mode fractures) is yet to be investigated.

The maximum throw $\delta_z^{max}$ along the faults in the models scales with S regardless of the cause of S variation, and increases with the volume contraction $\gamma$. The average ratio $\Omega = \delta_z^{max}/T$ of ~0.045, documented in the natural PFS, is reached in the models with maximum coupling between the layers at $\gamma \approx 7\%$, consistent with the available data. At low coupling and/or weaker substratum layer, $\gamma$ will be smaller. On the other hand, an increase in the residual internal friction angle $\phi^{res}$ and/or material ductility will result in a higher $\gamma$ required to achieve the same $\Omega$ value.

The maximum fault throw in the models is observed at the active layer surface, while in nature it is located below the middle of this layer. One of the reasons for this discrepancy could be the fact that faulting in nature is synsedimentation, at least partly. The variation of the material contraction rate with depth could be another possibility. This indicates avenues for further modeling development.

# Appendix A

## A1. Elastoplastic Model

To model the inelastic straining and its localization in Layer 1 (Figure 3), we use the constitutive model from Chemenda et al. (2021) summarized here. The model incorporates the following composite yield function

$$F = \begin{cases} F_t = \sigma_t + \sigma_3 & \text{if } \sigma_1 \le \sigma_{\text{int}} \\ F_s = \alpha_c \sigma_3 + \sigma_c - \sigma_1 & \text{if } \sigma_1 > \sigma_{\text{int}} \end{cases} \qquad (A1)$$

where $F_t$ and $F_s$ are the yield functions for the tensile and shear failure mechanisms, respectively; $\sigma_t$ and $\sigma_c$ are the tensile and unconfined compression strengths, respectively; $\sigma_{\text{int}}$ is the abscissa of the intersection point of the envelopes $F_t = 0$ and $F_s = 0$ in $(\sigma_1, \sigma_3)$ space; $\alpha_c$ is the dimensionless material (friction) parameter.

Stresses are defined following the rock mechanics convention with $\sigma_1 \ge \sigma_2 \ge \sigma_3$, where positive values represent compressive stress and negative values represent tension. All stresses are considered effective. The parameters $\sigma_t$, $\sigma_c$, $\sigma_{\text{int}}$, and $\alpha_c$ are functions of the accumulated inelastic volume strain $\epsilon^p$

$$\begin{aligned} \sigma_t(\epsilon^p) &= \sigma_t^0 \kappa(\epsilon^p), \quad \sigma_c(\epsilon^p) = \sigma_c^0 \kappa(\epsilon^p), \quad \sigma_{\text{int}}(\epsilon^p) = \sigma_{\text{int}}^0 \kappa(\epsilon^p), \\ \alpha_c &= (\alpha_c^0 - \alpha_c^{\text{res}})\kappa(\epsilon^p) + \alpha_c^{\text{res}} \end{aligned} \qquad (A2)$$

where $\kappa(\epsilon^p) = \max[0, (1 - \epsilon^p/\epsilon^0)]$; $\epsilon^0$ is the $\epsilon^p$ value at which both tensile and compression strengths reach zero during failure; $\sigma_t^0, \sigma_c^0, \sigma_{\text{int}}^0$, and $\alpha_c^0$ are the initial values of $\sigma_t, \sigma_c, \sigma_{\text{int}}$, and $\alpha_c^0$; $\alpha_c^{\text{res}}$ is the residual (at $\epsilon^p \ge \epsilon^0$) value of $\alpha_c$ related to the Coulomb internal friction angle $\phi$ as

$$\alpha_c = \frac{1 + \sin(\phi)}{1 - \sin(\phi)}. \qquad (A3)$$

the Coulomb cohesion c is related to $\sigma_c$ as

$$c = \frac{\sigma_c}{2\sqrt{\alpha_c}} \qquad (A4)$$
The tensile and shear failure mechanisms are thus coupled in this model by synchronous evolution of the strengths $\sigma_t$, $\sigma_c$, and $\alpha_c$ during the inelastic straining; at $\epsilon^p = \epsilon^0 \sigma_t$ and $\sigma_c$ reach zero and $\alpha_c$, its residual value $\alpha_c^{res}$.

The plastic potential function $\Phi$ reads

$$ \Phi = \begin{cases} \Phi_t = \sigma_3 & \text{if } \sigma_1 \le \sigma_{\text{int}} \\ \Phi_s = \beta_c \sigma_3 - \sigma_1 & \text{if } \sigma_1 > \sigma_{\text{int}} \end{cases} \qquad (A5) $$

where $\Phi_t$ and $\Phi_s$ are the plastic potentials for the tensile and shear failure mechanisms, respectively; $\beta_c$ is the material coefficient depending on $\epsilon^p$ and linearly decreasing with $\epsilon^p$ from its initial value $\beta_c^0$ to 1

$$ \beta_c = \beta_c^0 + (1 - \beta_c^0)(1 - \kappa) \qquad (A6) $$

This coefficient defines the dilatancy factor $\beta = d\epsilon^p/d\bar{\gamma}^p$, where $d\bar{\gamma}^p$ is the incremental inelastic equivalent or maximum shear strain. For the latter case

$$ \beta = \frac{\beta_c - 1}{\beta_c + 1} \qquad (A7) $$

as can be obtained from (Equations A5, A6), and the flow rule

$$ d\epsilon_i^p = d\lambda \frac{\partial \Phi_s}{\partial \sigma_i}, \qquad (A8) $$

where $d\lambda$ is the nonzero scalar function ($d\lambda$ is different for the different failure mechanisms introduced), $i = 1, 2, 3$. It follows from Equation A7 that when $\beta_c = 1$, $\beta = 0$.

## A2. Interfaces

In Models 9 and 10 (Figure 8), we introduced interfaces between Layers 1 and 2 with the imposed shear strength $\tau_0^{int}$, which is defined by the interface friction coefficient $\mu^{int}$ and cohesion $c^{int}$ as

$$ \tau_0^{\text{int}} = \mu^{\text{int}}\sigma_n + c^{\text{int}} \qquad (A9) $$

where $\sigma_n$ is the effective normal-to-interface stress, which is equal to the normal vertical stress $\sigma_{zz}$. In the models presented here, $\mu^{int}$ is set to 0, so $\tau_0^{int} = c^{int}$.

The normal ($k_n$) and shear ($k_s$) interface stiffnesses must also be specified. Since the interface's elastic deformation should be minimal, these stiffness values should be very large. However, the increase in $k_n$ and $k_s$ slows down the solution due to FLAC 3D's mass scaling based on stiffnesses. Therefore, one should use the smallest stiffnesses consistent with the small interface deformation. Using stiffness values that are around 10 times the equivalent stiffness of the stiffest numerical grid zone along the interface appears to be a good rule of thumb (Itasca, 2024)

$$ k_n = k_s = 10 \max \left[ \frac{E(\nu - 1)}{(\nu + 1)(2\nu - 1)\Delta z} \right] \qquad (A10) $$

where $E$ and $\nu$ are Young's modulus and Poisson's ratio, respectively; $\Delta z$ is the numerical grid element (zone) size, and the “max” function signifies that the maximum value over all zones adjacent to the interface is to be used.

## A3. Regularization of the Models

To reduce the mesh dependence of strain localization modeling, we incorporate the mesh size $\Delta z$ into the constitutive formulation through the following simple relationship (Chemenda, 2019)
$$ \epsilon^0 = \bar{\epsilon}^0 \left( \frac{C}{\Delta z} \right) \qquad (A11) $$

where $\bar{\epsilon}^0$ is the reference $\epsilon^0$ value that can correspond to real material, and $C$ is a constant coefficient. Similar conditions have been previously introduced by Bažant & Oh, 1983, Crook et al., 2003, and others, and helps mitigate the influence of $\Delta z$ on the rate of inelastic energy dissipation within deformation bands. This condition however does not allow to eliminate the loss of hyperbolicity or ellipticity of the field equations in dynamic or static regimes, respectively, at the bifurcation point. This loss results in the ill-posedness of the initial value problem during constitutive instability leading to strain localization. Treatment of this issue requires enrichment of constitute models using nonlocal/gradient (e.g., Chen & Schreyer, 1987; Needleman, 1988; Vardoulakis, 1989) or viscous (e.g., Duvaut & Lions, 1972; Needleman, 1988; Simo et al., 1988; Wang et al., 1997) regularizations. We apply here a viscous regularization using Duvaut-Lions scheme (Duvaut & Lions, 1972) as detailed by Chemenda et al. (2022). Integration of this regularization into the constitutive formulation leads to the emergence of a dimensionless nonnegative parameter $\eta$, which influences the inelastic deformation in a manner opposite to the classical viscosity coefficient. The viscosity effect is maximum when $\eta \to 0$, and $\eta = 0$ corresponds to purely elastic behavior. The viscosity effect is inexistent when $\eta = \infty$, which corresponds to a purely elasto-plastic or inviscid model. The $\eta$ value to be applied is defined by preliminary numerical simulations, as described, for example, by Chemenda et al. (2022).

Nomenclature

<table><tr><td>E</td><td>Young's modulus</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>G</td><td>Shear modulus</td></tr><tr><td>σ<sub>i</sub></td><td>Principal effective stresses, compression positive, i = 1, 2, 3</td></tr><tr><td>x, y, z</td><td>Axes defined in Figure 3</td></tr><tr><td>σ<sub>xx</sub>, σ<sub>yy</sub>, σ<sub>zz</sub></td><td>Effective normal stresses</td></tr><tr><td>σ<sub>v</sub> and σ<sub>h</sub></td><td>Effective lithostatic and horizontal stresses</td></tr><tr><td>σ<sub>t</sub>(ε<sup>p</sup>)</td><td>Tensile strength which is a function of ε<sup>p</sup></td></tr><tr><td>σ<sub>t</sub><sup>0</sup></td><td>Initial (at ε<sup>p</sup> = 0) tensile strength</td></tr><tr><td>σ<sub>c</sub>(ε<sup>p</sup>)</td><td>Unconfined (uniaxial) compression strength</td></tr><tr><td>σ<sub>c</sub><sup>0</sup></td><td>Initial value of σ<sub>c</sub></td></tr><tr><td>c</td><td>Cohesion defined by Equation A4</td></tr><tr><td>α<sub>c</sub>, β<sub>c</sub></td><td>Material parameters in Equations A3 And A6</td></tr><tr><td>β<sub>c</sub><sup>0</sup></td><td>Initial value of β<sub>c</sub></td></tr><tr><td>α<sub>c</sub><sup>0</sup></td><td>Initial value of α<sub>c</sub></td></tr><tr><td>α<sub>c</sub><sup>res</sup></td><td>Residual value of α<sub>c</sub></td></tr><tr><td>φ</td><td>Internal friction angle related to α<sub>c</sub> by Equation A3</td></tr><tr><td>φ<sup>res</sup></td><td>Residual value of φ</td></tr><tr><td>β</td><td>Dilatancy factor related to β<sub>c</sub> and defined by Equation A7</td></tr><tr><td>C, ε</td><td>Defined in Equation A11</td></tr><tr><td>ε<sup>p</sup></td><td>Accumulated inelastic (plastic) volume strain</td></tr><tr><td>ε<sup>0</sup></td><td>ε<sup>p</sup> at which σ<sub>t</sub>, σ<sub>c</sub>, and β reach zero and α<sub>c</sub> reaches α<sub>c</sub><sup>res</sup></td></tr><tr><td>ρ</td><td>Rock density</td></tr><tr><td>γ</td><td>Diagenetically induced volume contraction strain</td></tr><tr><td>φ</td><td>Porosity</td></tr><tr><td>Δφ</td><td>Porosity reduction due to the diagenetically induced volume contraction</td></tr><tr><td>η</td><td>Parameter controlling viscous regularization</td></tr><tr><td>T</td><td>Thicknesses of Layer 1 or active layer (Figure 3)</td></tr><tr><td>T<sub>s</sub></td><td>Thicknesses of Layer 2 or substratum layer (Figure 3)</td></tr><tr><td>L, D</td><td>Length and width of the models</td></tr><tr><td>S</td><td>Fault spacing</td></tr><tr><td>Δz</td><td>Numerical grid element (zone) size</td></tr><tr><td>δ<sub>z</sub></td><td>Fault throw</td></tr></table>
$\delta_z^{\max}$ Maximum fault throw

$k_n$ and $k_s$ Normal and shear stiffnesses of the interfaces

$\tau_0^{\text{int}}$ Interface shear strength

$F$, $F_t$, $F_s$ Yield functions defined in Equation A1

$\Phi$, $\Phi_t$, $\Phi_s$ Plastic potential functions defined in Equation A5

## Data Availability Statement

The Commercial Numerical Geomechanical Code FLAC3D (version 7.00.160, released on 22 August 2023) used in this study is available from Itasca Consultants at https://www.itascacg.com/software/flac3d. The FLAC3D project files and supporting data sets associated with this research are available at (Chemenda, 2025). These files include all necessary input scripts to reproduce the numerical simulations, as detailed in the manuscript. No external data sets were used or created for this research, except for the seismic data used in Figure 1. These data are available upon request from the Norwegian Offshore Directorate FactPages, a repository for all offshore exploration data in Norway, at https://factmaps.sodir.no/map/.

## References

Barnett, J. A., Mortimer, J., Rippon, J. H., Walsh, J. J., & Watterson, J. (1987). Displacement geometry in the volume containing a single normal fault. *AAPG Bulletin*, 71(8), 925-937. https://doi.org/10.1306/948878ed-1704-11d7-8645000102c1865d

Baudon, C., & Cartwright, J. A. (2008). 3D seismic characterisation of an array of blind normal faults in the Levant Basin, Eastern Mediterranean. *Journal of Structural Geology*, 29(6), 746-760. https://doi.org/10.1016/j.jsg.2007.12.002

Bažant, Z. P., & Oh, B. H. (1983). Crack band theory for fracture of concrete. *Matériaux et Construction*, 16(3), 155-177. https://doi.org/10.1007/BF02486267

Berndt, C., Jacobs, C., Evans, A., Gay, A., Elliott, G., & Long, D. (2012). Kilometre-scale polygonal seabed depressions in the Hatton Basin, NE Atlantic Ocean: Constraints on the origin of polygonal faulting. *Marine Geology*, 332-334, 126-133. https://doi.org/10.1016/j.margeo.2012.05.009

Carrier, B., Vandamme, M., Pellenq, R. J. M., Bornert, M., Ferrage, E., Hubert, F., & Van Damme, H. (2016). Effect of water on elastic and creep properties of self-standing clay films. *Langmuir*, 32(5), 1370-1379. https://doi.org/10.1021/acs.langmuir.5b04135

Cartwright, J., & Lonergan, L. (1996). Volumetric contraction during the compaction of mudrocks: A mechanism for the development of regional-scale polygonal fault systems. *Basin Research*, 8(2), 183-193. https://doi.org/10.1111/j.1365-2117.1996.tb00119.x

Cartwright, J. A. (1994). Episodic basin-Wide hydrofracturing of overpressured early Cenozoic mudrock sequences in the North Sea Basin. *Marine and Petroleum Geology*, 11(5), 587-607. https://doi.org/10.1016/0264-8172(94)90003-5

Cartwright, J. A. (1997). Polygonal extensional fault systems: A new class of structure formed during the early compaction of shales. In *Fluid flow and transport in rocks: Mechanisms and effects* (pp. 35-56). Springer.

Cartwright, J. A. (2011). Diagenetically induced shear failure of fine-grained sediments and the development of polygonal fault systems. *Marine and Petroleum Geology*, 28(9), 1593-1610. https://doi.org/10.1016/j.marpetgeo.2011.06.004

Cartwright, J. A., & Dewhurst, D. N. (1998). Layer-bound compaction faults in fine-grained sediments. *Geological Society of America Bulletin*, 110(10), 1242-1257. https://doi.org/10.1130/0016-7606(1998)110<1242:LCBFIF>2.3.CO;2

Cartwright, J. A., Huuse, M., & Aplin, A. C. (2007). Seal bypass systems. *AAPG Bulletin*, 91(8), 1141-1166. https://doi.org/10.1306/04090705181

Casey, B. B. A. (2014). *The consolidation and strength behavior of mechanically compressed fine-grained sediments* (Doctoral dissertation). Massachusetts Institute of Technology.

Chemenda, A. I. (2007). The formation of shear-band/fracture networks from a constitutive instability: Theory and numerical experiment. *Journal of Geophysical Research*, 112(B11), B11404. https://doi.org/10.1029/2007JB005028

Chemenda, A. I. (2019). Origin of regular networks of joints: Experimental constraints, theoretical background, and numerical modeling. *Journal of Geophysical Research: Solid Earth*, 124(8), 9164-9181. https://doi.org/10.1029/2019jb017454

Chemenda, A. I. (2022). Bed thickness-dependent fracturing and inter-bed coupling define the nonlinear fracture spacing-bed thickness relationship in layered rocks: Numerical modeling. *Journal of Structural Geology*, 165, 104741. https://doi.org/10.1016/j.jsg.2022.104741

Chemenda, A. I. (2025). FLAC3D input files for numerical modeling of Polygonal Fault Systems in sediments [Software]. Zenodo. https://doi.org/10.5281/zenodo.14964763

Chemenda, A. I., Hassani, R., & Fan, J. (2022). Numerical modeling of the opening mode fracturing emanating from deformation localization in layered rocks. *Computers and Geotechnics*, 147, 104774. https://doi.org/10.1016/j.compgeo.2022.104774

Chemenda, A. I., Lamarche, J., Matonti, C., Bazalgette, L., & Richard, P. (2021). Origin of strong nonlinear dependence of fracture (joint) spacing on bed thickness in layered rocks: Mechanical analysis and modeling. *Journal of Geophysical Research: Solid Earth*, 126(3), e2020JB020656. https://doi.org/10.1029/2020JB020656

Chemenda, A. I., & Mas, D. (2016). Dependence of rock properties on the Lode angle: Experimental data, constitutive model, and bifurcation analysis. *Journal of the Mechanics and Physics of Solids*, 96, 477-496. https://doi.org/10.1016/j.jmps.2016.07.007

Chen, Z., & Schreyer, H. L. (1987). Simulation of soil-concrete interfaces with nonlocal constitutive models. *Journal of Engineering Mechanics*, 113(11), 1665-1677. https://doi.org/10.1061/(ASCE)0733-9399(1987)113:11(1665)

Clausen, J. A., Gabrielsen, R. H., Reksnes, P. A., & Nysaether, E. (1999). Development of intraformational (Oligocene-Miocene) faults in the northern North Sea: Influence of remote stresses and doming of Fennoscandia. *Journal of Structural Geology*, 21(10), 1457-1475. https://doi.org/10.1016/S0191-8141(99)00083-8

Cosenza, P., Giot, R., & Hedan, S. (2023). Elastic moduli of clay minerals and their aggregates: A review. *Applied Clay Science*, 236, 106878. https://doi.org/10.1016/j.clay.2023.106878

## Acknowledgments

This work has been supported by the French government through the France 2030 investment plan managed by the Agence Notionale de la Recherche, as part of the Université Côte d'Azur's Initiative of Excellence (reference ANR-15-IDEX-01). We thank J. Cartwright and the anonymous reviewer for their useful suggestions.
Crook, A. J. L., Willson, S., Yu, J. G., & Owen, R. (2003). Computational modelling of the localized deformation associated with borehole breakout in quasi-brittle materials. *Journal of Petroleum Science and Engineering*, **38**(3-4), 177-186. https://doi.org/10.1016/S0920-4105(03)00034-1

Cuadros, J., & Linares, J. (1996). Experimental kinetic study of the smectite-to-illite transformation. *Geochimica et Cosmochimica Acta*, **60**(3), 439-453. https://doi.org/10.1016/0016-7037(95)00407-6

Daigle, H., Cook, A., Fang, Y., Bihani, A., Song, W., & Flemings, P. B. (2020). Gas-driven tensile fracturing in shallow marine sediments. *Journal of Geophysical Research: Solid Earth*, **125**(12), e2020JB020835. https://doi.org/10.1029/2020JB020835

Davies, R. J., & Ireland, M. T. (2011). Initiation and propagation of polygonal fault arrays by thermally triggered volume reduction reactions in siliceous sediment. *Marine Geology*, **289**(1-4), 150-158. https://doi.org/10.1016/j.margeo.2011.10.003

Desrues, J. (2003). Experimental characterization of localized deformation in geomaterials. In *Advanced mathematical and computational geomechanics* (pp. 77-106). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-540-45079-5_4

Dewhurst, D. N., Cartwright, J., & Lonergan, L. (1999). The development of polygonal fault systems by syneresis of colloidal sediments. *Marine and Petroleum Geology*, **16**(8), 793-810. https://doi.org/10.1016/S0264-8172(99)00036-9

Duvaut, G., & Lions, J. L. (1972). *Les inéquations en mécanique et en physique*. Dunod.

Eberl, D. D., Velde, B., & McCormick, T. (1993). Synthesis of illite-smectite from smectite at earth surface temperatures and high pH. *Clay Minerals*, **28**(1), 49-60. https://doi.org/10.1180/claymin.1993.028.1.06

Fortin, J., Stanchits, S., Dresen, G., & Gueguen, Y. (2009). Acoustic emissions monitoring during inelastic deformation of porous sandstone: Comparison of three modes of deformation. *Pure and Applied Geophysics*, **166**(5), 823-841. https://doi.org/10.1007/s00024-009-0479-0

Gay, A., & Berndt, C. (2007). Cessation/reactivation of polygonal faulting and effects on fluid flow in the Vøring Basin, Norwegian Margin. *Journal of the Geological Society*, **164**(1), 129-141. https://doi.org/10.1144/0016-76492005-194

Gay, A., Lopez, M., Berndt, C., & Seranne, M. (2007). Geological controls on focused fluid flow associated with seafloor seeps in the Lower Congo Basin. *Marine Geology*, **244**(1-4), 68-92. https://doi.org/10.1016/j.margeo.2007.06.003

Gay, A., Lopez, M., Cochonat, P., & Sermondadaz, G. (2004). Polygonal faults-furrows system related to early stages of compaction—Upper Miocene to recent sediments of the Lower Congo Basin. *Basin Research*, **16**(1), 101-116. https://doi.org/10.1111/j.1365-2117.2004.00219.x

Gay, A., Padron, C., Meyer, S., Beaufort, D., Oliot, E., Lallemand, S., et al. (2021). Elongated giant seabed polygons and underlying polygonal faults as indicators of the creep deformation of Pliocene to recent sediments in the Grenada Basin, Caribbean Sea. *Geochemistry, Geophysics, Geosystems*, **22**(12), e2021GC009809. https://doi.org/10.1029/2021GC009809

Ghalayini, R., & Eid, C. (2020). Using polygonal layer-bound faults as tools to delimit elastic reservoirs in the Levant Basin offshore Lebanon. *AAPG Bulletin*, **104**(3), 629-656. https://doi.org/10.1306/07151918155

Ghalayini, R., Homberg, C., Daniel, J.-M., & Nader, F. H. (2017). Growth of layer-bound normal faults under a regional anisotropic stress field. *Geological Society, London, Special Publications*, **439**(1), 57-78. https://doi.org/10.1144/SP439.13

Goulty, N. R. (2008). Geomechanics of polygonal fault systems: A review. *Petroleum Geoscience*, **14**(4), 389-397. https://doi.org/10.1144/1354-079308-781

Goulty, N. R., & Swarbrick, R. E. (2005). Development of polygonal fault systems: A test of hypotheses. *Journal of the Geological Society*, **162**(4), 587-590. https://doi.org/10.1144/0016-764904-118

Hamilton, E. L. (1971). Elastic properties of marine sediments. *Journal of Geophysical Research*, **76**(2), 579-604. https://doi.org/10.1029/JB076i002p00579

Hamilton, E. L. (1976). Variations of density and porosity with depth in deep-sea sediments. *Journal of Sedimentary Research*, **46**(2), 280-300. https://doi.org/10.1306/212F6E38-2B24-11D7-8648000102C1865D

Hansen, D. M., Shimeld, J. W., Williamson, M. A., & Lykke-Andersen, H. (2004). Development of a major polygonal fault system in Upper Cretaceous chalk and Cenozoic mudrocks of the Sable Subbasin, Canadian Atlantic margin. *Marine and Petroleum Geology*, **21**(10), 1205-1219. https://doi.org/10.1016/j.marpetgeo.2004.07.003

Hantschel, T., & Kauerhauf, A. I. (2009). *Fundamentals of basin and petroleum systems modeling*. Springer.

He, C., Tang, C., Huang, D., & Shi, S. (2010). Polygonal faults in the Sanzhao sag of the Songliao Basin: Their significance in hydrocarbon accumulation. *Mining Science and Technology*, **20**(2), 300-305. https://doi.org/10.1016/S1674-5264(09)60196-2

Henriet, J. P., De Batist, M., & Verschuren, M. (1991). Early fracturing of Palaeogene clays, southernmost North Sea: Relevance to mechanisms of primary hydrocarbon migration. *Generation, accumulation and production of Europe's hydrocarbons*, 1, 217-227. https://doi.org/10.1007/978-94-011-3186-8_16

Ho, S., Carruthers, D., & Imbert, P. (2016). Insights into the permeability of polygonal faults from their intersection geometries with linear chimneys: A case study from the Lower Congo Basin. *Carnets de Géologie*, **16**(2), 17-31. https://doi.org/10.1016/j.margeo.2007.06.003

Ho, S., Hovland, M., Blouet, J. P., Wetzel, A., Imbert, P., & Carruthers, D. (2018). Formation of linear planform chimneys controlled by preferential hydrocarbon leakage and anisotropic stresses in faulted fine-grained sediments, offshore Angola. *Solid Earth*, **9**(6), 1437-1468. https://doi.org/10.5194/se-9-1437-2018

Hu, S., Rotevatn, A., Jackson, C., Li, W., & Wu, X. (2024). Geometry and evolution of polygonal fault systems under a regionally anisotropic stress field: Insights from 3D seismic analysis of the Qiongdongnan Basin, NW South China Sea. *Basin Research*, **36**(1), e12855. https://doi.org/10.1111/bre.12855

Hustoft, S., Mienert, J., Bünz, S., & Nouzé, H. (2007). High-resolution 3D-seismic data indicate focused fluid migration pathways above polygonal fault systems of the mid-Norwegian margin. *Marine Geology*, **245**(1-4), 89-106. https://doi.org/10.1016/j.margeo.2007.07.004

Ikari, M. J., Kameda, J., Saffer, D. M., & Kopf, A. J. (2015). Strength characteristics of Japan Trench borehole samples in the high-slip region of the 2011 Tohoku-Oki earthquake. *Earth and Planetary Science Letters*, **412**, 35-41. https://doi.org/10.1016/j.epsl.2014.12.030

Ikari, M. J., & Kopf, A. J. (2011). Cohesive strength of clay-rich sediment. *Geophysical Research Letters*, **38**(16), L16309. https://doi.org/10.1029/2011GL048318

Ikari, M. J., Kopf, A. J., Hüpers, A., & Vogt, C. (2018). Lithologic control of frictional strength variations in subduction zone sediment inputs. *Geosphere*, **14**(2), 604-625. https://doi.org/10.1130/GES01543.1

Ishii, E., Sanada, H., Iwatsuki, T., Sugita, Y., & Kurikami, H. (2011). Mechanical strength of the transition zone at the boundary between opal-A and opal-CT zones in siliceous rocks. *Engineering Geology*, **122**(3-4), 215-221. https://doi.org/10.1016/j.enggeo.2011.06.006

Itasca. (2024). *Fast Lagrangian analysis of continua in 3-dimensions, version 5.0, manual*. Itasca. Retrieved from https://docs.itascacg.com/flac3d700/flac3d/docproject/source/flac3dhome.html

Jackson, C., Carruthers, D., Mahlo, S., & Briggs, O. (2014). Can polygonal faults help locate deep-water reservoirs? *AAPG Bulletin*, **98**(9), 1717-1738. https://doi.org/10.1306/06031413128

King, J. J., & Cartwright, J. A. (2020). Ultra-slow rates of polygonal fault systems. *Geology*, **48**(5), 473-477. https://doi.org/10.1130/G47221.1
King, J. J., Roberts, D. T., Cartwright, J. A., & Levell, B. K. (2022). Numerical modelling of the growth of polygonal fault systems. *Journal of Structural Geology*, **162**, 104679. https://doi.org/10.1016/j.jsg.2022.104679

King, M. S. (1966). Wave velocities in rocks as a function of changes in overburden pressure and pore fluid saturants. *Geophysics*, **31**(1), 50–73.
https://doi.org/10.1190/1.1439761

Laurent, D., Gay, A., Baudon, C., Berndt, C., Soliva, R., Planke, S., et al. (2012). High-resolution architecture of a polygonal fault interval inferred from geomodel applied to 3D seismic data from the Gjallar Ridge, Vøring Basin, Offshore Norway. *Marine Geology*, **332**, 134–151. https://doi.org/10.1016/j.margeo.2012.03.003

Lonergan, L., Cartwright, J., & Jolly, R. (1998). The geometry of polygonal fault in Tertiary mudrocks of the North Sea. *Journal of Structural Geology*, **20**(5), 529–548. https://doi.org/10.1016/S0191-8141(97)00094-7

Mair, K., Main, I., & Elphick, S. (2000). Sequential growth of deformation bands in the laboratory. *Journal of Structural Geology*, **22**(1), 25–42.
https://doi.org/10.1016/S0191-8141(99)00124-8

Mesri, G., & Shahien, M. (2003). Residual shear strength mobilized in first-time slope failures. *Journal of Geotechnical and Geoenvironmental Engineering*, **129**(1), 12–31. https://doi.org/10.1061/(ASCE)1090-0241(2003)129:1(12)

Minaeian, V., Dewhurst, D. N., & Rasouli, V. (2017). Deformational behaviour of a clay-rich shale with variable water saturation under true triaxial stress conditions. *Geomechanics for Energy and the Environment*, **11**, 1–13. https://doi.org/10.1016/j.gete.2017.05.001

Mondol, N. H., Bjørlykke, K., & Jahren, J. (2008). Experimental compaction of clays: Relationship between permeability and petrophysical properties in mudstones. *Petroleum Geoscience*, **14**(4), 319–337. https://doi.org/10.1144/1354-079308-758

Mondol, N. H., Bjørlykke, K., Jahren, J., & Høeg, K. (2007). Experimental mechanical compaction of clay mineral aggregates—Changes in physical properties of mudstones during burial. *Marine and Petroleum Geology*, **24**(5), 289–311. https://doi.org/10.1016/j.marpetgeo.2007.03.006

Mondol, N. H., Jahren, J., Bjørlykke, K., & Brevik, I. (2008). Elastic properties of clay minerals. *The Leading Edge*, **27**(6), 758–770. https://doi.org/10.1190/1.2939658

Morgan, D. A., Cartwright, J. A., & Imbert, P. (2015). Perturbation of polygonal fault propagation by buried pockmarks and the implications for the development of polygonal fault systems. *Marine and Petroleum Geology*, **65**, 157–171. https://doi.org/10.1016/j.marpetgeo.2015.03.026

Morrow, C. A., Moore, D. E., & Lockner, D. A. (2017). Frictional strength of wet and dry montmorillonite. *Journal of Geophysical Research: Solid Earth*, **122**(5), 3392–3409. https://doi.org/10.1002/2017JB014081

Neagu, R. C., Cartwright, J., & Davies, R. J. (2010). Measurement of diagenetic compaction strain: Quantitative analysis of fault plane dip. *Journal of Structural Geology*, **32**(5), 641–655. https://doi.org/10.1016/j.jsg.2010.03.005

Needleman, A. (1988). Material rate dependence and mesh sensitivity in localization problems. *Computer Methods in Applied Mechanics and Engineering*, **67**(1), 69–85. https://doi.org/10.1016/0045-7825(88)90069-2

Okuda, H., Ikari, M. J., Roesner, A., Stanisłowski, K., Hüpers, A., Yamaguchi, A., & Kopf, A. J. (2021). Spatial patterns in frictional behavior of sediments along the Kumano transect in the Nankai Trough. *Journal of Geophysical Research: Solid Earth*, **126**(11), e2021JB022546. https://doi.org/10.1029/2021JB022546

Ortega, J. A., Ulm, F. J., & Abousleiman, Y. (2007). The effect of the nanogranular nature of shale on their poroelastic behavior. *Acta Geotechnica*, **2**(3), 155–182. https://doi.org/10.1007/s11440-007-0030-3

Ottozen, N. S., & Runesson, K. (1991). Properties of discontinuous bifurcation solutions in elasto-plasticity. *International Journal of Solids and Structures*, **27**(4), 401–421. https://doi.org/10.1016/0020-7683(91)90131-X

Peltonen, C., Marcussen, Ø., Bjørlykke, K., & Jahren, J. (2009). Clay mineral diagenesis and quartz cementation in mudstones: The effects of smectite to illite reaction on rock properties. *Marine and Petroleum Geology*, **26**(6), 887–898. https://doi.org/10.1016/j.marpetgeo.2008.08.001

Rabinowitz, H. S., Savage, H. M., Skarbek, R. M., Ikari, M. J., Carpenter, B. M., & Colletti, C. (2018). Frictional behavior of input sediments to the Hikurangi Trench, New Zealand. *Geochemistry, Geophysics, Geosystems*, **19**(9), 2973–2990. https://doi.org/10.1029/2018GC007653

Rice, J. R. (1976). The localization of plastic deformation. In W. T. Koiter (Ed.), *Theoretical and applied mechanics: Proceedings of the 14th international congress on theoretical and applied mechanics* (Vol. 1, pp. 207–220). North-Holland Publishing Co.

Rice, S. B., Freund, H., Huang, W. L., Clouse, J. A., & Isaacs, C. M. (1995). Application of Fourier transform infrared spectroscopy to silica diagenesis: The opal-A to opal-CT transformation. *Journal of Sedimentary Research*, **65**(4a), 639–647. https://doi.org/10.2110/jsr.65.639

Roberts, D. T., Crook, A. J. L., Cartwright, J. A., Profit, M. L., & Rance, J. M. (2014). The evolution of polygonal fault systems: Insights from geomechanical forward modeling. In *ARMA US rock mechanics/geomechanics symposium* (pp. ARMA-2014). ARMA.

Rudnicki, J. W., & Rice, J. R. (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids*, **23**(6), 371–394. https://doi.org/10.1016/0022-5096(75)90001-0

Saffer, D. M., Frye, K. M., Marone, C., & Mair, K. (2001). Laboratory results indicating complex and potentially unstable frictional behavior of smectite clay. *Geophysical Research Letters*, **28**(12), 2297–2300. https://doi.org/10.1029/2001GL012869

Sanada, H., Niunoya, S., Matsui, H., & Fujii, Y. (2009). Influences of sedimentary history on the mechanical properties and microscopic structure change of Horonobe siliceous rocks. *Journal of MMJ*, **125**(10–11), 521–529. https://doi.org/10.2473/journalofmmij.125.521

Shin, H., Santamarina, J. C., & Cartwright, J. A. (2008). Contraction-driven shear failure in compacting uncemented sediments. *Geology*, **36**(12), 931–934. https://doi.org/10.1130/G25168A.1

Shin, H., Santamarina, J. C., & Cartwright, J. A. (2010). Displacement field in contraction-driven faults. *Journal of Geophysical Research*, **115**(B7), B07102. https://doi.org/10.1029/2009JB006572

Simo, J. C., Kennedy, J. G., & Govindjee, S. (1988). Non-smooth multisurface plasticity and viscoplasticity: Loading/unloading conditions and numerical algorithms. *International Journal for Numerical Methods in Engineering*, **26**(10), 2161–2185. https://doi.org/10.1002/nme.1620261004

Stipp, M., Rolfs, M., Kitamura, Y., Behrmann, J. H., Schumann, K., Schulte-Kortnack, D., & Feeser, V. (2013). Strong sediments at the deformation front, and weak sediments at the rear of the Nankai accretionary prism, revealed by triaxial deformation experiments. *Geochemistry, Geophysics, Geosystems*, **14**(11), 4791–4810. https://doi.org/10.1002/ggge.20293

Stuevold, L. M., Faerseth, R. B., Arnesen, L., Cartwright, J., & Möller, N. (2003). Polygonal faults in the Ormen Lange field, Møre basin, offshore mid-Norway. *Geological Society, London, Special Publications*, **216**(1), 263–281. https://doi.org/10.1144/GSL.SP.2003.216.017

Sulem, J., & Ouffroukh, H. (2006). Shear banding in drained and undrained triaxial tests on a saturated sandstone: Porosity and permeability evolution. *International Journal of Rock Mechanics and Mining Sciences*, **43**(2), 292–310. https://doi.org/10.1016/j.ijrmms.2005.07.001

Sun, Q., Wu, S., Lü, F., & Yuan, S. (2010). Polygonal faults and their implications for hydrocarbon reservoirs in the southern Qiongdongnan Basin, South China Sea. *Journal of Asian Earth Sciences*, **39**(6), 470–479. https://doi.org/10.1016/j.jseas.2010.05.002

Tenzer, R., & Gladkikh, V. (2014). Assessment of density variations of marine sediments with ocean and sediment depths. *The Scientific World Journal*, **14**, 823296–823299. https://doi.org/10.1155/2014/823296
Tewksbury, B. J., Hogan, J. P., Kattenhorn, S. A., Mehrtens, C. J., & Tarabees, E. A. (2014). Polygonal faults in chalk: Insights from extensive exposures of the Khoman Formation, Western Desert, Egypt. *Geology*, *42*(6), 479–482. https://doi.org/10.1130/G35204.1

Tran, T. P. H., Bouissou, S., Chemenda, A., Ambre, J., Vacher, P., & Michel, P. (2019). Initiation and evolution of a network of deformation bands in a rock analogue material at brittle-ductile transition. *Rock Mechanics and Rock Engineering*, *52*(3), 737–752. https://doi.org/10.1007/s00603-018-1641-8

Tuckwell, G. W., Lonergan, L., & Jolly, R. J. H. (2003). The control of stress history and flaw distribution on the evolution of polygonal fracture networks. *Journal of Structural Geology*, *25*(8), 1241–1250. https://doi.org/10.1016/s0191-8141(02)00165-7

Turrini, L., Jackson, C. A. L., & Thompson, P. (2017). Seal rock deformation by polygonal faulting, offshore Uruguay. *Marine and Petroleum Geology*, *86*, 892–907. https://doi.org/10.1016/j.marpetgeo.2017.06.021

Vanorio, T., Prasad, M., & Nur, A. (2003). Elastic properties of dry clay mineral aggregates, suspensions, and sandstones. *Geophysical Journal International*, *155*(1), 319–326. https://doi.org/10.1046/j.1365-246X.2003.02052.x

Vardoulakis, I. (1989). Shear-banding and liquefaction in granular materials on the basis of a Cosserat continuum theory. *Ingenieur-Archiv*, *59*(2), 106–113. https://doi.org/10.1007/BF00534852

Volpi, V., Camerlenghi, A., Hillenbrand, C. D., Reesco, M., & Ivaldi, R. (2003). Effects of biogenic silica on sediment compaction and slope stability on the Pacific margin of the Antarctic Peninsula. *Basin Research*, *15*(3), 339–363. https://doi.org/10.1046/j.1365-2117.2003.00213.x

Walsh, J. J., & Watterson, J. (1988). Analysis of the relationship between displacements and dimensions of faults. *Journal of Structural Geology*, *10*(3), 239–247. https://doi.org/10.1016/0191-8141(88)90057-0

Wang, W. M., Sluys, L. J., & De Borst, R. (1997). Viscoplasticity for instabilities due to strain softening and strain-rate softening. *International Journal for Numerical Methods in Engineering*, *40*(20), 3839–3864. https://doi.org/10.1002/nme.1620402012

Wang, Z., Wang, H., & Cates, M. E. (2001). Effective elastic properties of solid clays. *Geophysics*, *66*(2), 428–440. https://doi.org/10.1190/1.1444936

Watterson, J., Walsh, J., Nicol, A., Neil, P. A. R., & Bretan, P. G. (2000). Geometry and origin of a polygonal fault system. *Journal of the Geological Society of London*, *157*(1), 151–162. https://doi.org/10.1144/jgs.157.1.151

White, J. E., & Sengbush, R. L. (1953). Velocity measurements in near-surface formations. *Geophysics*, *18*(1), 54–69. https://doi.org/10.1190/1.1437901

Wrona, T., Magee, C., Jackson, C., Huuse, M., & Taylor, K. (2017). Kinematics of polygonal fault systems: Observations from the Northern North Sea. *Frontiers in Earth Science*, *5*, 101. https://doi.org/10.3389/feart.2017.00101

Xia, Y., Yang, J., Chen, Y., Lu, S., Wang, M., Deng, S., et al. (2022). A review of the global polygonal faults: Are they playing a big role in fluid migration? *Frontiers in Earth Science*, *9*, 786915. https://doi.org/10.3389/feart.2022.00786915

Zeng, Z., Zhu, H., Yang, X., & Cao, X. (2022). Three-dimensional seismic analysis of a polygonal fault system (PFS) in the Northern Carnarvon Basin, Australia: Implications for fluid flow migration and gas hydrate system. *Journal of Petroleum Science and Engineering*, *215*, 110602. https://doi.org/10.1016/j.petrol.2022.110602