Subsurface salt movement in the absence of external tectonic forces can affect contemporaneous sediment deposition, mask allocyclic signals, and deform older strata. We used a discrete element model (DEM) to better understand salt-related modification of a sedimentary sequence with an increasing sedimentation rate. This permitted quantification of thinning rates and analysis of the lateral extent of synkinematic layers. Results show realistic evolution of salt-related faults, defining two salt-withdrawal basins, beyond which strata are undeformed. Thinning of stratigraphy is four times greater between the salt flank and crest than between the undeformed zone and flank, confirming an intense zone of halokinetic modulation adjacent to the diapir. Early, slowly aggrading layers are isolated within the salt-withdrawal basin and strongly influenced by salt growth, whereas later, quickly aggrading layers are more laterally extensive, matching inferences made from subsurface and outcrop data. Halokinetic modulation reduces up the stratigraphic section, mirroring observations around the Pierce diapirs, in the North Sea, offshore UK. Our DEM provides quantitative insights into the dynamic interplay between halokinetic and allocyclic controls on salt-stratigraphic relationships.
Seafloor salt topography affects the routing of sedimentary systems and their resultant stratigraphic architecture (e.g., Davison et al., 2000). Thickness variations and stratigraphic pinch-outs adjacent to salt are observed in the subsurface (e.g., Trusheim, 1960), at outcrop (e.g., Giles and Rowan, 2012), and with physical models (e.g., Dooley et al., 2007). However, such stratigraphic relationships remain hard to image and quantify in seismic data due to poor velocity control, steep dips, and deformation, and they are rare and incomplete in outcrop due to dissolution (Jones and Davison, 2014). Salt movement is often proposed to cause stratigraphic modification of allocyclic signatures (Mayall et al., 2010; Giles and Rowan, 2012). However, because natural analogs represent a “snapshot” in time of a complicated evolution controlled by the dynamic interplay of parameters (sedimentation rate, salt supply, etc.), it is often impossible to discern exact controls. Models help to isolate such controls.
Finite-element modeling (FEM; models based on continuum methods) has enhanced our understanding of salt flow (Gemmer et al., 2005; Hamilton-Wright et al., 2019), minibasin dynamics (Fuchs et al., 2011; Fernandez et al., 2020), and stress evolution (Nikolinakou et al., 2018). FEM treats the overburden as continuous frictional- or viscous-plastic material, which is useful for studying ductile deformation, but it prevents the development of brittle structures.
To address this, and complement current methods and models, we used a novel discrete element modeling (DEM) approach to isolate sedimentation rate (SR) from salt rise rate. DEM enables us to study dynamic system evolution and is an efficient way to test multiple scenarios without requiring computationally intense remeshing, which is typical in FEM (Abe et al., 2011).
We tested the hypothesis that halokinesis can modulate stratigraphy and mask allocyclic signals. Using a DEM, we modeled salt growth under an increasing SR. We aimed to (1) identify and quantify thickness variations and pinch-out geometries in a salt-influenced depositional system; (2) test diapir response to increasing SR; and (3) analyze lateral and temporal variations in stratigraphic architecture and compare results to natural analogs.
DEM has been used to study deformation processes (Hardy and Finch, 2006), and salt tectonics (Pichel et al., 2017, 2019). For a discussion of the equations governing DEM, and DEM limitations, see Finch et al. (2003, 2004) and Botter et al. (2014). Our modeled parameters build on previous work (summarized in Table S1 in the Supplemental Material1; Finch et al., 2003, 2004; Pichel et al., 2017, 2019).
Our model consists of a 4.5 × 1.5 km rectangle containing ∼44,500 randomly distributed elements of varying radii (5.25–9.75 m), to reduce failure in preferred orientations. A 150 m salt layer was overlain by nine 150-m-thick overburden layers (Fig. 1). This work was not intended to investigate initial diapir evolution (e.g., Vendeville and Jackson, 1992). To simplify a complex three-dimensional process into a two-dimensional model, we assumed a radially symmetric diapir, emplaced by an earlier phase of diapirism (Pichel et al., 2017, 2019). This diapir was initially assumed to actively rise, driven by loading of the salt and a salt-overburden density contrast, without regional tectonics (Jackson and Hudec, 2017). Active rise is restricted by roof thickness; diapir height must be more than 66%–75% of the surrounding overburden for growth to occur (Schultz-Ela et al., 1993). We modeled a sinusoidal 750-m-wide (maximum width), 1050-m-tall (70% of the overburden) diapir.
We imposed an upward motion of 0.023 mm yr–1 to all salt elements, based on North Sea diapir rise rates (Fig. 2; Davison et al., 2000), similar to a volumetric salt supply approach (Peel et al., 2020). The model was run for 4.6 m.y., where one time step represented 100 yr. The initial 2.2 m.y. interval was used for model equilibration, creating topography before addition of sediment in three 0.8 m.y. stages (SR1, SR2, SR3). SR increases through time (Fig. 1), from 0.15 mm yr–1 (SR1) to 0.45 mm yr–1 (SR3), were based on North Sea Cenozoic SR (De Haas et al., 1996). Increasing SR could reflect absolute increases (in this model) or relative increases with respect to salt rise (in the North Sea analog). Following addition of sediment (after 2.2 m.y. in our model), diapir rise is akin to passive diapirism and is modulated by synkinematic sedimentation (Rowan and Giles, 2020).
In our model, layer relative strength is defined by its breaking separation, where particles are bonded until this is exceeded (Donzé et al., 1994). Particle motion is frictionless and cohesionless (Finch et al., 2003). Overburden breaking separation and density increase linearly with depth from 0.023 to 0.027 and from 2.4 g cm–3 to 2.6 g cm–3, respectively, to mimic compaction. Salt and sediment have densities of 2.2 g cm–3 and 2.3 g cm–3, respectively. Ductile rock-salt is represented by a viscosity of 1.1 × 109 Pa·s, and viscous-plastic behavior is mimicked by weakening the breaking separation by an order of magnitude (0.001 for salt; Pichel et al., 2017). Over 100 models were run; here, we present a two-dimensional DEM to analyze diapiric modulation of stratigraphy.
Fault distribution and deformation extent were interpreted using a relative displacement algorithm (Fig. 2). Linear zones of high displacement represent faults, most of which are located in the early diapiric sequence and are extensional with a maximum throw of 135 m (Fig. 3). The faults furthest from the diapir mark the edge of two salt-withdrawal basins (∼1150 m wide; Fig. 3). These “deformation zone”–bounding faults, offsetting layers 1–7, are associated with smaller synthetic and antithetic structures; outside withdrawal basins, strata are undeformed. Above the crest, steeply dipping reverse faults related to inner-arc shortening of 24% are localized in layers 6–9 and extend into the synkinematic strata between the flank and crest (Figs. 2 and 3). Away from the crest, layers 8 and 9, and the lowermost synkinematic strata, are dominated by small-scale normal faults.
Layers A–L represent the synkinematic sedimentary sequence. Layer A fills to a base level and was excluded from analysis. Layers B–D represent low SR and thicken into the withdrawal basin by 3.5% (over 900 m from the undeformed zone), at 0.004%/m, before pinching-out toward the diapir. Layers E–H represent intermediate SR. Layer E is not extensive, thinning toward the flank (across 1250 m) by 0.053%/m. Layers F–H extend across the model but thin at 0.045%/m over 1500 m, from the undeformed zone to crest (Fig. 3). Layers I–L represent high SR and extend across the crest. Thinning rates halve, from 0.032%/m in I to 0.015%/m in L, showing a rapid decrease in salt influence with time. This reduced modulation of the late synkinematic sedimentary sequence is a product of increased SR.
The entire synkinematic package (A–L) thins by 59% across 1500 m at 0.04%/m. Thickening of 12% (over 900 m) at 0.014%/m is observed between the undeformed zone and the withdrawal basin. Thinning of 0.11%/m between the withdrawal basin and flank (350 m) and thinning of 0.16%/m between the flank and crest (250 m) are observed. This reflects thickening into the withdrawal basin and rapid thinning toward the crest.
DEM is not intended as a substitute for FEM, but to complement these approaches in salt tectonics. DEM can replicate localized fault growth, evolution, and propagation (Figs. 1–3), making it appropriate for studying the interactions of halokinesis, stratigraphy, and tectonics (Pichel et al., 2017, 2019). FEM packages have limited capacity to generate faults during simulations, unless they are predefined (e.g., Heidari et al., 2016; Nikolinakou et al., 2018). DEM treats the contacts between elements as potential displacement surfaces; this is realistic when studying areas with multiscalar salt-related faulting.
As well as seismically resolvable faults, outcrop and borehole data indicate brittle deformation is significant in salt basins (e.g., Koestler and Ehrmann, 1991). Our DEM replicates this brittle deformation (Fig. 2); as in nature, extreme thinning and termination of layers are in part accommodated by small-scale displacements. Understanding sub-seismic–scale fault distribution is important for predicting reservoir compartmentalization and seal integrity in the subsurface. DEM is therefore advantageous due to its replication of diapir-related brittle deformation, while FEM is useful for modeling salt flow (e.g., Albertz and Ings, 2012).
Halokinetic Modulation of Stratigraphy
In the absence of salt tectonics, loading and compaction are considered to be the drivers for subsidence, and structures are controlled by regional tectonics. The undeformed zone shows synkinematic strata thicken upwards during an allocyclic increase of SR (Fig. 3B). Here, salt topography results in stratigraphic onlap and thickness variations, which are absent in models without salt influence (Figs. 1 and 3), confirming assumptions from subsurface and outcrop examples (e.g., Davison et al., 2000; Giles and Rowan, 2012).
Lateral variations in synkinematic geometry are significant, defining the relationship between early diapiric and synkinematic strata. Synkinematic strata (441 m, 2.4 m.y.) conformably overlie early diapiric strata in the undeformed zone (Fig. 3C) with no evidence of salt influence. Subparallel strata (496 m, 2.4 m.y.) thicken into the salt-withdrawal basin to form a subtle, low-angle unconformity over the early diapiric sequence (Fig. 3C). Synkinematic strata (303 m, 1.6 m.y.) dip at a shallower angle than early diapiric strata, creating an angular unconformity at the flank, which is difficult to differentiate from a regional unconformity (Trudgill, 2011). Parallel deposition (183 m, 1.4 m.y.) is observed over the crest. Here, the early diapiric to synkinematic contact is conformable, but condensed (Jackson and Hudec, 2017).
When viewed in isolation, the crest suggests increasing SR occurs later and over a shorter time period than the rest of the model, indicating that salt growth can modify and mask allocyclic signatures. However, from the undeformed zone to the crestal regions (lateral distance of ∼1500 m), 258 m of strata representing 1.0 m.y. are missing, as sediment was either not deposited, or it was removed by erosion or failure during growth. This missing stratigraphy clearly demonstrates the effect of halokinesis on allocyclic signal preservation.
Thinning rates vary laterally, according to diapir proximity. For the withdrawal basin to flank and the flank to crest, thinning rates are 3 and 4 times greater, respectively, than from the undeformed zone to crest. Growing salt is associated with additional accommodation, as evidenced by thickening into salt-withdrawal basins (Fig. 3B). Up the stratigraphic section, late synkinematic thickness variability is reduced as halokinetic modulation is also reduced by the increased SR, healing topography (Fig. 3D). Such relationships are assumed in data where strata cannot be discerned close to the diapir. Our model confirms this and predicts lateral and temporal thinning rate evolution and stratigraphic geometries in the presence of growing salt. This is vital in crest and flank regions where imaging of overlying and onlapping synkinematic strata is difficult (e.g., Jones and Davison, 2014; Charles and Ryzhikov, 2015).
Subsurface Analog Comparison
The Pierce diapirs, in the North Sea offshore UK, show evidence for halokinetic modulation of Jurassic–Pleistocene strata (Fig. 4; Birch and Haynes, 2003), which, despite a longer, more complex evolution, show geometrical similarities to our DEM. The tectonostratigraphic history of the Pierce diapirs, spanning ∼200 m.y., includes: Jurassic reactive-active diapirism, Cretaceous–Cenozoic passive diapirism, and Alpine contractional growth (Carruthers et al., 2013). Strata are near horizontal ∼2 km away from the diapirs, but they are upturned adjacent to the diapirs, comparable with model results (Figs. 2 and 3). Like the DEM, layers are initially isolated on either side of the diapirs and then thin and become upturned adjacent to them (Jurassic–Eocene). Rapid thickness changes and pinch-outs suggest diapir growth was fastest during the Paleocene–Eocene (Carruthers et al., 2013). Core from flanks shows a reduction in reservoir quality toward the diapirs, similar to outcrop observations (Cumberpatch et al., 2021). Generation of faults over the crest (Fig. 4) and in Cenozoic strata (Carruthers et al., 2013) supports the results of physical models (e.g., Davison et al., 1993) and our DEM.
Strata eventually extend across the Pierce diapirs (Pliocene–Pleistocene), analogous to DEM thinning rates reducing with time. Modeled thinning rate reductions are due to increased SR, healing salt-cored topography. The Pierce diapirs show a similar evolutionary trend but reflect a relative increase in SR with respect to diapir rise rates (Carruthers et al., 2013). The halokinetic influence over the diapirs is reduced from the Oligocene to Pliocene; the controlling variable is difficult to isolate but could be source layer depletion, resulting in burial of salt bathymetry (Fig. 4; Birch and Haynes, 2003).
Hence, our model results can apply to basins with absolute or relative increases in SR and could represent a slowing of diapir growth due to regional tectonic quiescence or depletion of the salt source layer, as well as absolute increases in SR. Combined analysis demonstrates the importance of understanding, and where possible isolating, local (salt layer variations) and regional (tectonic and SR) controls when disentangling salt basin evolution.
Our two-dimensional model shows the stratigraphic architecture expected in a salt-influenced sedimentary system with increasing sedimentation rates. Layers initially onlap and form depocenters isolated by salt-cored topography. With increasing sedimentation, layers traverse the entire structure but continue to thicken into the salt-withdrawal basins and thin toward the diapir. Stratigraphically, an upward reduction in halokinetic influence is demonstrated, which is observed in natural analogs.
Laterally, we observe a zone of 1150 m to either side of the diapir that is influenced by halokinesis, and beyond this zone, strata are undeformed. Synkinematic strata on the flanks thin toward the diapir, but the thinning rate is four times greater between the flank and crest compared to that between the undeformed zone and crest. While the overall modulation style close to the diapir is stratigraphic thinning, stratigraphic thickening is observed in the salt-withdrawal basins.
Similar increased stratigraphic modulation is noted toward the Pierce diapirs, UK North Sea, despite a different tectonostratigraphic history to our DEM, confirming previous assumptions of halokinetic modulation of allocyclic processes. This DEM approach can be employed to test different scenarios for the development of halokinetic stratigraphy and fault patterns.
Cumberpatch's Ph.D. is part of the Natural Environment Research Council (NERC) Centre for Doctoral Training (CDT) in Oil and Gas (grant NEM00578X/1). We thank Petroleum Geo-Services (PGS, Oslo, Norway) for permission to publish an image from MegaSurvey Plus 3D seismic data; Mads Huuse and Ross Grant for discussions; Bruce Trudgill, Jürgen Adam, Frank Peel, Evey Gannaway Dalton, and Michael Hudec for insightful reviews; and editor James Schmitt for handling of this manuscript.