Halokinetic modulation of sedimentary thickness and architecture: A numerical modelling approach

Subsurface salt flow can deform overlying strata and influence contemporaneous sedimentary systems. Studying salt‐sediment interactions is challenging in the subsurface due to poor imaging adjacent to salt, and in the field due to the dissolution of halite. Discrete Element Modelling provides an efficient and inexpensive tool to model stratigraphy and deformation around salt structures, which is advantageous over other modelling techniques as it realistically recreates brittle processes such as faulting. Six 2D experiments were run representing 4.6 Myr to determine the effect of salt growth on syn‐kinematic stratigraphy. Halokinetic deformation of stratigraphic architecture was assessed by varying sediment input rates through time. Results show the realistic formation and evolution of salt‐related faults which define a zone of halokinetic influence ca. 3 times the width of the initial diapir. Outside of this, early diapiric and syn‐kinematic stratigraphy are undeformed. Within this zone, syn‐kinematic strata are initially isolated into primary salt withdrawal basins, onlapping and thinning towards the salt‐cored high. In most models, syn‐kinematic strata eventually thin across and cover the diapir roof. Thinning rates are up to six times greater within 350 m of the diapir, compared to further afield, and typically decrease upwards (with time) and laterally (with distance) from the diapir. Outputs are compared to a subsurface example from the Pierce field, UK North Sea, which highlights the importance of considering local fluctuations in diapir rise rate. These can create stratigraphic architectures that may erroneously be interpreted to represent increases/decreases in sedimentation rate. Exposed examples, such as the Bakio diapir, northern Spain, can be used to make inferences of the expected depositional facies, below model resolution. Our models aid the prediction of sedimentary unit thickness and thinning rates and can be used to test interpretations arising from incomplete or low‐resolution subsurface and outcrop data when building geological models for subsurface energy.

Detailed depositional facies models of halokinetically influenced depositional systems (or portions of such systems) benefit from calibration with outcrop analogues (e.g. Banham & Mountney, 2013a, 2013b, 2014Lerche & Petersen, 1995;Madof et al., 2009;Rodriguez et al., 2020), such that they can be more widely and confidently applied to salt-influenced basins globally. However, outcrop examples possess their own set of uncertainties; exposed examples are often limited in the rock record, largely due to the dissolution of halite forming the core of the salt bodies . Rare exposures provide sub-seismic scale facies information for shallow-(e.g. Giles & Lawton, 2002;, deep-(e.g. Poprawski et al., 2014Poprawski et al., , 2016 and non-marine stratigraphy (e.g. Banham & Mountney, 2013a, 2013b, 2014Ribes et al., 2015Ribes et al., , 2017. Many field examples are small in size compared to subsurface basins and therefore provide only small-scale details of, for example, sedimentary structures and stratal stacking patterns, rather than the larger, basin-scale tectonostratigraphic context of salt-sediment interactions provided by integrated subsurface datasets. While useful, each subsurface or field example represents a unique record of the ratio of relative salt rise and sedimentation rate. Also, subsurface and field examples provide only one snapshot in time ( Figure 1). Physical models have an advantage in recreating the evolution of specific subsurface analogues Dooley et al., 2013Dooley et al., , 2015Dooley et al., , 2020Ferrer et al., 2017;Roma et al., 2018) and studying sediment gravity flow distribution and evolution adjacent to salt topography (Gaullier & Vendeville, 2005;Sellier & Vendeville, 2009;Soutter et al., 2021).
A number of remaining questions can be addressed by taking a numerical modelling approach that allows us to modify and isolate the key controls on salt-sediment interactions. These include: (a) How does salt topography influence depositional systems, and thus depositional facies, and how does this vary laterally and temporally? (b) How do unconformities, onlap contacts, and faults and fractures vary in salt-influenced settings? (c) How does sedimentation rate influence the width of the roof and basin salt-related deformation zones? and (d) How do stratigraphic thinning rates Highlights • We use a Discrete Element Model (DEM) to determine the effect of salt growth on stratigraphic architecture for different sedimentation rates and patterns.
• Results show realistic formation and evolution of salt-related faults which define a zone of halokinetic influence, outside of which stratigraphy are undeformed.
• Stratigraphy initially thins and onlaps' the salt cored structure and is isolated to individual minibasins, before being laterally extensive across the model but continuing to thin over the residual diapir's crest.
• Halokinetic modulation is up to six times more intense within 350 m of the diapir, compared to further afield, halokinetic influence is shown to reduce laterally and temporarily.
• The DEM results are compared to analogues from the Eastern Central Graben and Basque Cantabrian Basin.

|
associated with salt growth vary with sedimentation rate and distance from the salt structure? Most recent advances in the numerical modelling of saltrelated deformation use finite element models (FEM), which are based on continuum methods. Such studies have focused on the physical conditions required for the initiation and development of diapirism (Chemia et al., 2008;Fernandez & Kaus, 2015;Fuchs et al., 2011;Gemmer et al., 2004Gemmer et al., , 2005Hamilton-Wright et al., 2019;Nikolinakou et al., 2017;Peel et al., 2020;Poliakov et al., 1993), the stratigraphic architecture of subsiding minibasins (Fernandez et al., 2020;Sylvester et al., 2015;Wang et al., 2017), reconstructing the evolutionary history of salt-affected stratigraphy (Ismail-Zadeh et al., 2001Pichel et al., 2017Pichel et al., , 2019, and salt-related stress (and strain) analysis Luo et al., 2012Luo et al., , 2017Nikolinakou et al., 2012Nikolinakou et al., , 2014aNikolinakou et al., , 2014bNikolinakou et al., , 2018. FEM treats the overburden as a continuous frictional-plastic or viscous-plastic material, which prevents the development of realistic brittle deformation, for example, fracturing and faulting, in overburden stratigraphy ( Figure 1). Most FEM packages currently have limited capacity to generate faults during simulations, and therefore either have no faults, or faults that are pre-defined at the start of the model simulation (e.g. Heidari et al., 2016;Nikolinakou et al., 2014aNikolinakou et al., , 2018. FEM is, however, advantageous for studying ductile deformation and salt flow dynamics (Albertz & Ings, 2012).
Discrete Element Models (DEM), which essentially treat the contact between each element as a slip surface, are able to replicate spontaneous, realistic, localised fault nucleation and growth Pichel et al., 2017Pichel et al., , 2019 and are therefore appropriate for studying the interactions between salt-related topography, sedimentation, and stratigraphic evolution. DEMs do have limitations, including the need for careful calibration of element parameters (Abe et al., 2011;Botter et al., 2014) and the limited number of elements and duration of simulations (Zhu et al., 2008).
Despite this, DEM provides a quick, efficient and inexpensive method to investigate system evolution through time (Allen & Tildesley, 1987;Donzé et al., 1994;Finch et al., 2003Finch et al., , 2004. It is, therefore, possible to test several scenarios and collect structural growth and syn-kinematic sediment thinning rate data during deformation. This helps to improve the recognition of the processes involved with salt-sediment interactions. DEM studies have recently been adapted to salt tectonics so that elements representing salt behave as viscous-plastic materials (Pichel et al., 2017(Pichel et al., , 2019. Here, we use a DEM approach to understand how sedimentation rate affects stratal geometries in salt basins experiencing diapirism (Figure 1). First, we generate a baseline model with no sediment fill, to determine the effect of salt growth on the early diapiric 'overburden' sequence. For consistency in terminology, in line with recent work (Cumberpatch, Finch, et al., 2021), we define the 'early diapiric sequence' as the layers deposited prior to our simulation; they are discordant to the diapir and thus represent the earlier syn-kinematic strata related to an initial phase of diapirism that is assumed to have emplaced the diapir into our model . We then vary sedimentation rates and patterns to study how these control the stratigraphic record of halokinesis in the later 'syn-kinematic' stratigraphy. The aims of this study are to: (a) investigate variable syn-kinematic sedimentation rates adjacent to a dynamic salt diapir using a DEM; (b) quantify near-diapir thinning rates and how this is controlled by varying sedimentation rates and patterns; and (c) compare the results to subsurface and field analogues to test the validity of our approach and our model predictions.
DEMs offer advantages over other numerical methods in that scale is not a restriction, complex re-meshing is not required, and results are easily reproducible. DEMs are fundamentally discontinuous, and therefore each element simulates the specific physical properties of a given rock. These conditions make it a fit-for-purpose method to quantitatively study syn-kinematic deformation.
DEM treats objects as assemblages of circular elements, connected by breakable elastic bonds through a 'repulsive-attractive' force. Elements remain bonded until their breaking separation (defined as the relative strength of the assemblage) is exceeded (Donzé et al., 1994;Finch et al., 2004). Once these bonds break, previously connected elements experience no further 'attractive' force, if these elements return to contact with each other a 'repulsive' force acts between them, preventing the healing of bonds (Finch et al., 2003(Finch et al., , 2004Hardy & Finch, 2006). Motion of elements is frictionless and cohesionless with elastoplastic behaviour (Finch et al., 2003;Hardy & Finch, 2006). Forces are resolved in the x and y directions and elements are subject to gravity (F g ) (Finch et al., 2003). The equations that define the inter-relationship of all forces acting on the DEM are: where F i,n corresponds to the total elastic force acting on an element, V represents the viscosity and Ẋ and Ẏ correspond to the velocity of the element in the x-and y-directions. A viscous term is added to counteract the elastic behaviour within a closed system, making it ideal for studying quasisteady state tectonic processes (Finch et al., 2004;Pichel et al., 2017Pichel et al., , 2019. For a comprehensive description of the equations governing DEM, see Mora and Place (1994), Finch et al. (2003, and Finch (2005, 2006). Pichel et al. (2017Pichel et al. ( , 2019 recently used a DEM to model salt tectonics for the first time, studying regional-scale compressional salt tectonics (Pichel et al., 2017), and the effect of base salt relief on salt flow and overburden deformation styles (Pichel et al., 2019). Cumberpatch, Finch, et al. (2021) adapted these to focus on the modulation of stratigraphy by salt diapir growth. In these models, the elements representing salt were adjusted so they behave as a viscous-plastic material in order to represent rock-salt (Pichel et al., 2017(Pichel et al., , 2019. This requires inter-element interactions to be adjusted so they behave macroscopically as viscous-plastic materials and deform microscopically by dislocation creep, which is expected for dry rock salt (Pichel et al., 2017(Pichel et al., , 2019Spiers et al., 1990). This does not naturally scale to completely reproduce salt, which typically deforms on a spectrum of mechanisms including diffusion and dislocation creep Pichel et al., 2017;Spiers et al., 1990), and usually contains traces of brines (Warren, 1999(Warren, , 2006, but is considered a satisfactory assumption for studying salt tectonics. Pichel et al. (2017) tested breaking separations using biaxial compression tests. Values representing tenths of a model unit (e.g. 0.05) develop defined fault segments, and produced responses typical of brittle materials; these values are, therefore, used to represent overburden sediment in this study. Values representing hundredths of a unit (e.g. 0.001), however, show a minor elastic component (e.g. F i , α ≈ 0), representing ductile viscous-plastic materials that accumulate strain without significant stress variations. Consequently, a breaking separation of 0.001 for salt elements is used in this study (see also Pichel et al., 2017), and other physical (Spiers et al., 1990) and numerical (Li & Urai, 2016) experiments of salt deformation. By this value, we ensure salt element motion is entirely controlled by the viscosity and gravity of the system (viscous-plastic behaviour). The scaled viscosity of the salt is 1.1 × 10 9 Pa s, which is lower than its real-world viscosity (10 17 -10 18 Pa s) (Hudec & Jackson, 2007;, but works as a reasonable approximation when compared with physical models (e.g. Dooley et al., 2009Dooley et al., , 2012Vendeville et al., 1995).
The modelled media in this study consists of a simulated 4.5 km by 1.5 km box, to ensure that the outer boundaries do not influence the structural evolution of the model centre. The media consist of an undeformable base and ca. 44,500 elements with varying radii (0.175-0.35 units, representing 5.25-9.75 m); these are randomly distributed to reduce failure in preferential orientations within the matrix. A 150 m-thick salt layer is overlain by nine coloured 150 m thick early diapiric overburden layers. A salt density of 2.2 g/cm 3 is used, to mimic a slightly impure halite composition, comparable to many global salt basins (e.g. Grant et al., 2019;Hudec & Jackson, 2007;Warren, 1999Warren, , 2006 and previous models (Cumberpatch, Finch, et al., 2021;Pichel et al., 2017Pichel et al., , 2019. We do not investigate the initial stages of diapir evolution, which have been well-studied elsewhere (Costa & Vendeville, 2002;Hudec & Jackson, 2007;Trusheim, 1960;Vendeville & Jackson, 1992), and instead, focus on late stage diapir growth. Therefore, we simplify a complicated three-dimensional process into a two-dimensional model, where we assume a linear salt wall or radially symmetric diapir was emplaced by an earlier phase of diapirism; such an assumption is used in other numerical (Pichel et al., 2017(Pichel et al., , 2019 and physical models (Davison et al., 1993;Dooley et al., 2009Dooley et al., , 2012. This allows us to focus on the coupled deformation-sedimentation characterising the late stage of diapir growth, when a tall pre-existing diapir may be rejuvenated by compression (not modelled here) or rise due to buoyancy. During the experiment, the diapir is assumed to grow by passive diapirism (initially by halokinetic active rise; , driven by the pressure of the overburden on the salt source layer and to a lesser extent by the density difference between the salt and the overlying stratigraphy . Such growth can happen in the absence of regional tectonics, although mild far-field compression or extension, which can enhance deformation rates, are likely in most natural settings . Active rise is retarded by roof thickness and strength, and salt viscosity. Therefore, diapir height must be >66%-75% of the surrounding overburden thickness for substantial halokinesis to occur and the roof thickness must be <750 m (Schultz-Ela et al., 1993). In adhering to this rule, we invoke an individual sinusoid, 750 m (base width) wide, 1,050 m (70% of the 1,500 m overburden) tall diapir and thus a 450 m roof ( Figure 2). The diapir geometries | 2577 EAGE CUMBERPATCH ET Al.
used in our models are comparable to those observed in natural examples . Overburden breaking separation (relative strength) increases with depth linearly from 0.023 to 0.027, representing increasing rock strength with depth. An overburden density of 2.4-2.6 g/cm 3 is used, similar to natural and modelled examples and increases with depth (e.g. Dooley et al., 2009Dooley et al., , 2012Fuchs et al., 2011).

F I G U R E 3
Model 3, the intermediate aggradation experiment with a sedimentation rate of 0.3 mm/year. Diapir growth rate is continuous and constant throughout (0.023 mm/year). Outputs begin at the start of the experiment during a stage of no sedimentation (T = 0 Myr, circle S0 (sedimentation stage 0)). At T = 2.2 Myr (b) the diapir has generated surface topography and sedimentation begins to be added (S0/S1). Subsequent outputs are generated at 0.4 Myr intervals across the three equal stages of sedimentation (open circles S1-S3) from 2.2 to 4.6 Myr. Active sediment input is constant in stages 1-3 and the boundaries are highlighted by the grey circle (e.g. T = 3 Myr represents the end of sedimentation 1 and the beginning of sedimentation event 2). The final output is at 4.6 Myr (h), following the final sedimentation from S3. Due to minimal deformation in the outer section of all models, for clarity, subsequent figures will focus on only the central 3,000 m around the diapir (shown by the grey outline, In nature, diapir growth periods are hugely variable in duration, ranging from 100,000s of years to 100s of millions of years . Our focus on the latestage of diapir growth and using seismic stratigraphic observations from natural examples (Grimstad, 2016;Oluboyo et al., 2014) and run times of previous numerical models (Pichel et al., 2017(Pichel et al., , 2019, supports our experimental run times of 46,000 timesteps with a timestep equivalent to 100 years (4.6 Myr in total). We impose an upwards motion of 0.023 mm/year, based on North Sea diapirs (Davison, Alsop, Birch et al., 2000) to all elements representative of salt to mimic diapir growth rate ( Figure 2); this aims to replicate the volumetric salt supply rate (Q) described by Peel et al. (2020). The diapir grows for 2.2 Myr (22,000 time steps) to allow the model to equilibrate creating seabed or surface topography, prior to the addition of sediment ( Figure 3). Sediment is added from 2.2 Myr with a constant density of 2.3 g/cm 3 , in agreement with natural examples of near seabed sediment (Rider & Kennedy, 2018;Tenzer & Gladkikh, 2014) and a breaking separation of 0.023. Sediment is added in three 0.8 Myr (8,000 time steps) stages (S1-S3, Table 1). Sedimentation rates in nature are extremely variable (Sadler, 1981). Here, the sedimentation rate was varied between 0.15 and 0.45 mm/year to match Cenozoic rates measured in the North Sea (de Haas et al., 1996) and the North Atlantic (Whitman & Davies, 1979).
We present results from six experiments: a baseline zero sedimentation model, and five models with variable sedimentation rates and patterns (slow, intermediate and fast constant sedimentation, increasing and decreasing). Increasing and decreasing sedimentation rates are used to replicate the local advance and retreat of depositional sedimentary systems (progradation and retrogradation). Model set-up and parameters have been rigorously tested (Finch et al., 2003(Finch et al., , 2004Pichel et al., 2017Pichel et al., , 2019 and are summarised in Table S1.

| Model limitations
In addition to considering how modelled rock properties (density, viscosity, breaking separation; Table S1) scale to nature it is important to note that these values are often an oversimplification as they assume homogeneity in a given rock property for an entire layer (e.g. 4,500 m laterally or the entire 1,050 m tall diapir). The natural heterogeneity within stratigraphy, driven by depositional facies variability, differential diagenetic processes and products, and proximal to distal trends, are not incorporated in our models, which aim to replicate non-unique stratigraphic variations rather than specific, subtle, unique changes in depositional character. The complex three-dimensional processes occurring in salt basins are simplified into a two-dimensional model for this study. Therefore, we assume our models represent a cross-section through a three-dimensional linear salt wall or radially symmetric diapir (e.g. Cumberpatch, Finch, et al., 2021;Pichel et al., 2017Pichel et al., , 2019; this is an oversimplification based on the complex, often asymmetric, geometries of salt structures (Hudec & Jackson, 2007), but allows for the assumption that salt is continuously flowing in two dimensions, essentially being fed from out of the twodimensional plane (Tvedt et al., 2013(Tvedt et al., , 2016Vendeville et al., 1995). Modelling in two dimensions also assumes that all processes (such as salt withdrawal and stratigraphic bed rotation) are equal in all directions, which is an oversimplification (e.g. Dutta et al., 2016;Ismail-Zadeh et al., 2004;Mattson et al., 2020;. However, this is a suitable assumption in simple models that focus primarily on the role of sedimentation rate variability on the halokinetic depositional record. In order to prevent circular reasoning, our model inputs do not attempt to recreate a specific realworld diapir, but rather a simplified universally applicable structure. The absence of more complicated salt geometries (e.g. salt overhangs or welds) prevents direct comparison of the models to specific settings with complicated threedimensional salt structures (e.g. the Gulf of Mexico). However, our approach allows us to generate more general, possibly portable insights that are applicable to global salt basins. Finally, sedimentation rates are extremely variable and non-linear (Sadler, 1981); thus, when comparing to certain analogues, 'slow' and 'fast' sedimentation (stated in Table 1) should be taken as relative rates rather than absolute values. Sediments can be deposited above the diapir relief, suggesting that they are applicable to subaqueous settings, and are assumed to aggrade evenly, preventing direct comparison to deltaic systems (clinoforms).

| RESULTS
The model with an intermediate sedimentation rate (M3) is first presented (Figure 3) to examine the relationship between the rates of salt diapir rise and sedimentation. Subsequent sections describe and compare diapir growth, deformation and stratigraphic architectures across all models (Figures 4-8;  Tables 2 and 3).  Figure S2). Cold colours represent elements that are in contact with their original neighbour and hot colours indicate high displacement. This is used to show discontinuities and is a proxy for fault location. (c) (a) overlain with (b), interpreted with locations of discontinuities, highlighting the methodology used to interpret structures for all models. Note how the majority of the internal salt has remained in connection with its original neighbouring element and the radial faults are associated with salt withdrawal. The maximum displacement shows the relative movement of the salt, the neighbouring overburden, and the high mobility of layers in close proximity to the diapir. Note the lack of deformation outside of the salt withdrawal basin. Base of syn-kinematic sequence | 2581 EAGE CUMBERPATCH ET Al.

| Summary of temporal relationship between halokinesis and sedimentation
The initially horizontal basal salt layer thins adjacent to the diapir, and the diapir geometry changes from initially triangular/sinusoidal to vertically elongated during evolution (Compare H with A in Figure 3); this is consistent across all six models (Figures 3-6). Up to 45% thinning of the salt layer, ca. 450 m on either side of the diapir base, occurs in all final model outputs (Figure 3h), and is accommodated by thinning in the first 2.2 Myr of evolution, as the initial model equilibrates for diapir growth. This thinning appears not to continue throughout the model as the diapir is assumed to be fed by salt flow from out of the plane (in three dimensions), allowing diapir rise to be maintained, and preventing welding. The growth of the diapir is accompanied by withdrawal effects in the adjacent stratigraphy, which is indicated by thinning of the source layer and faulting in the basal layers (layers 1-3, Figure 3; M3 in Figure 5b). Salt withdrawal and evidence for upwards salt growth are shown by the basal part of the diapir narrowing between T = 0 and T = 2.2 Myr (0 to 22,000 time steps; A-B in Figure 3). Sediment is first introduced to all models after 2.2 Myr (22,000 timesteps; Figure 3b). Up to this time, the diapir has risen such that 'early diapiric' layers 4 and 5 are folded during its rise and thin dramatically towards the salt due to structural attenuation, assumed to be accommodated by layer-parallel slip. The early diapiric sequence, originally overlying the diapir (layers 6-9, Figure 3), are passively folded as the diapir rises. This results in 'postdepositional' layer thickening on the flanks and thinning over the crest (Figure 3b). Above the diapir, a topographic high is generated with associated faulting in layers 6-9 observed in T = 2.2 Myr (Figure 3b). The topographic high influences subsequent sedimentation (layers A-L) and is the focus of our study of stratigraphic modulation (Figure 3c-h). Sediment is added in Stage 1 (S1, T = 2.2 Myr, 22,000 time steps; Figure 3c), filling two salt withdrawal basins on either side of the diapir. For simplicity, following the methodology of Pichel et al. (2017), the first syn-kinematic layer, layer A, fills to a flat, standardised base level in all models, to allow consistent comparison across all models and subsequent layers are input by assigned sedimentation rate (in the case of M3, 0.3 mm/yr). Later in S1 (T = 3 Myr, 30,000 time steps, Figure 3d), the rate of deposition outpaces diapir rise. Layers extend across the salt rather than onlapping it (Figure 3c,d). Throughout diapir rise, the early diapiric sequence is rotated away from the diapir crest, thickening in the adjacent depocentres, again assumed to be due to structural layer-parallel slip. The diapir stem narrows throughout evolution (compare H and A, Figure 3). The thinning of early diapiric layers above the rising diapir continues into S2, where early diapiric layer 5 approaches vertical at the salt-sediment interface (Figure 3e,f). Upward movement of the diapir is associated with, and driven by, salt withdrawal underneath the basins, and increased displacement of faults at the base of the model (layers 1-4; Figure 3; M3 in Figure 5b).
In the final stage (S3), the early diapiric and syn-kinematic stratigraphy thin above the diapir crest, and are further rotated and thicken into the salt withdrawal basins (Figure 3g,h). Faulting is present above the diapir tip propagating through the early diapiric sequence and into the syn-kinematic sediment (M3 in Figure 5b). Throughout diapir evolution, the dip of the early diapiric layers increases towards the structure, and thus the overburden anticline steepens and narrows as sedimentation progresses (compare h with b in Figure 3). The deformation in the early diapiric overburden sequence, described here, is similar across all models (M1-M6).
The fault furthest from the diapir is taken as the edge of the halokinetically influenced part of the succession, which is ca. 1,150 m wide on either flank (from diapir centre to fault edge), resulting in an ca. 2,300 m zone of diapiric influence, in all models. Outside this zone, early diapiric and syn-kinematic strata appear undeformed. Salt mobility has a limited influence on sediments at the extremities of M3, which is consistent across all models, so in subsequent figures, only the central 3,000 m is shown (grey box, T = 4.6 Myr, Figure 3).
To permit comparison of stratigraphic variability across all models, subsequent figures (Figures 4-6) present all models (M1-M6) at the end of the experiment (T = 4.6 Myr). In the following section, we describe and compare diapir growth and roof folding, and stratigraphic architecture. In each case, we first present and discuss M1, the case where there is no sediment input for comparison with models in which sedimentation rate is constant (M2-M4), increasing (M5) and decreasing (M6).

| Diapir growth and roof folding
3.2.1 | Diapir rise with no sedimentation (M1) In M1, the diapir crest rises a total of 425 m, and the final width of the early diapiric folded roof, taken from the greatest change in dip in layer 9 is 961 m (Figure 4; Figure S1). These values represent the base-case to compare the effect of different sedimentation rates on the final geometry of the salt structure and the early diapiric and syn-kinematic stratigraphic architectures (Tables 2 and 3; Figures 5 and 6).

| Constant sedimentation rates (M2-M4)
Under slow (M2), intermediate (M3) and fast (M4) constant sedimentation rates the diapir rises by 393, 363 and 297 m, respectively ( Figure 5). Diapir growth compared to M1 is reduced by 8%, 15% and 30% for the different aggradation cases, respectively ( Table 2). The width of the early diapiric folded roof decreases from 961 m in M1 to 770 m (M2), 760 m (M3) and 734 m (M4), which accounts for a 20%-24% reduction relative to M1. Where syn-kinematic stratigraphy is present across the model (i.e. not in M2), the syn-kinematic folded roof is measured from the point within layer F where there is the greatest change in dip ( Figure S1). Layer F is chosen as it is the first layer that is laterally extensive across all models (M3-M6). The syn-kinematic folded roof is 839 and 890 m wide in M3 and M4, respectively. The syn-kinematic roofs are therefore 110% and 120% greater than the width of the early diapiric folded roof in the same models ( Figure S1).

| Variable sedimentation rates: Increasing
(M5) and decreasing (M6) sedimentation rates In M5, under increasing sedimentation rate, the diapir rises by 368 m, a decrease of 14% compared to M1. Under decreasing sedimentation rate conditions (M6), the diapir rises by 346 m, a 19% reduction when compared to the base case M1 (Table 2; Figures 6 and 7). The early diapiric folded roof width is reduced to ca. 750 m in both models, a decrease of 22% compared to M1 (Table 2; Figure 7). The syn-kinematic folded roof is 723 and 872 m in M5 and M6, respectively, representing a 4% reduction and 16% increase compared to the early diapiric folded roof in the same model. Note: Measurements are taken from model simulations (Figures 4-6). Initial diapir height (1,050 m at T = 0 Myr) is subtracted from final diapir height to give diapir growth. M2-M6 are compared to M1 to show the percentage reduction that occurs when different amounts of sediment are added. Early diapiric anticline measurements are taken from the stratal termination of layer A, in M1 they are taken from the top of the overburden at the location of the greatest change in the dip. For M2-M6, the width of the overburden is compared to M1 to show the percentage reduction in overburden anticline width that occurs when different amounts of sediment are added. Synkinematic folded roof measurements are taken from the location of greatest dip change in layer F ( Figure S1). Value is only calculated where layers extend across the entire structure (M3-M6). Ratio of syn-kinematic to early diapiric folded roof thickness is calculated. EAGE CUMBERPATCH ET Al.

| Effect of sedimentation rate on stratigraphic architecture
Here, we discuss the lateral extent, thinning rates, and termination styles of stratigraphy onto the topographic high for M2-M6. In all models, layer A fills to a fixed base level and therefore does not represent the sedimentation rate of S1, so the first layer described is layer B. Across all models, layer A heals the topography. This step reduces irregularities with the upper surface of layer 9 that might influence syn-kinematic sedimentation. This partial healing of topography, which is comparable to hemipelagic deposition prior to the activation of a depositional system, allows for direct comparison across all models (Cumberpatch, Finch, et al., 2021; Model U-C (%)  Figure 5, and are the same for all models to ensure direct comparison.  (Table 3) are discussed with reference to three points that remain fixed throughout all models: the crest, flank and undeformed zone (see uninterpreted M2 in Figure 5).

| Constant sedimentation rates (M2-M4)
Syn-kinematic layers (B-L) are not laterally extensive in M2; this model is defined by slow sedimentation rates. In this case, deposition is restricted to primary salt withdrawal basins, with no sedimentation occurring over the crest of the early diapiric anticline. All layers terminate adjacent to the diapir, with the uppermost layers (I-L) onlapping the remnant topography created by the layers below ( Figure 5). The entire stratigraphic package thins by 34%, at a rate of 0.029%/m from the undeformed zone to salt flank, before pinching-out towards the crest.
In M3, the model defined by intermediate sedimentation rates, the earliest syn-kinematic strata (layers B and C) are preserved only in the salt withdrawal basin, offset some distance from the diapir (Table 3). However, in contrast to M2, layers D-L are laterally extensive across the model, extending across the diapir crest (Table 3). The overall stratigraphic thinning for intermediate aggradation is 55% from undeformed zone to crest, at a rate of 0.037%/m, with over one-third of this occurring between the undeformed zone and flank, accounting for 19% thinning at a rate of 0.017%/m. Thinning rates of 0.1%/m, totalling 35% stratigraphic thinning (almost two-thirds of the overall thinning observed across M3) are observed from flank to crest. The thinning rate between the salt flank and the crest is 6.1 times greater than that between the undeformed section and the flank (Table 3).
In the fast sedimentation model, M4, all layers are extensive across the model except for layer B ( Figure 5; Table 3). The overall stratigraphic thinning for fast aggradation is 33% from the undeformed section to crest at a rate of 0.022%/m. Over one-third of this total thinning is between the undeformed section and the salt flank, accounting for 12% stratigraphic thinning at a rate of 0.01%/m, and the other almost two-thirds occur between the salt flank and the crest with a thinning rate of 0.06% accounting for 21% stratigraphic thinning. The thinning rate between the salt flank and the crest is 5.9 times greater than that between the undeformed section and the flank (Table 3). 3.3.2 | Variable sedimentation rates: Increasing (M5) and decreasing (M6) sedimentation rates In M5 (increasing sedimentation rate), layers B-D represent slow sedimentation; these layers are restricted to the salt withdrawal basin on either side of the structure and thintowards the salt-influenced topography before pinching out. Layers E-H were deposited under intermediate sedimentation rates. Layer E does not extend across the model and thins towards the topography between the undeformed section and the salt flank, before pinching out towards the crest. Layers F-H are extensive across the model. Layers I-L were deposited under the fast sedimentation rate, and are laterally extensive across the entire model ( Figure 6; Table 3). Overall thinning in the increasing sedimentation model accounts for 58% of total stratigraphic thinning at an average rate of 0.039%/m; two-thirds of this thinning takes place between the salt flank and the crest, accounting for 39% stratigraphic thinning at a rate of 0.11%/m. The remaining third of the thinning occurs between the undeformed section and the salt flank, at a rate of 0.017%/m accounting for 19% stratigraphic thinning. The thinning rate between the salt flank and the crest is 6.6 times greater than that between the undeformed section and the salt flank (Table 3).
In M6 (decreasing sedimentation rate), layers B-D represent fast sedimentation. Layer B is isolated to either side of the diapir and thins between the undeformed section and the flank before terminating towards the diapir (Table 3). Layers C and D are laterally extensive. Layers E-H represent intermediate sedimentation and are deposited across the entire model. Layers I-L are deposited under slow sedimentation and are extensive across the whole model, but thin markedly over the crest (Figure 6; Table 3). The overall stratigraphic thinning for the decreasing sedimentation rate model is 52% at a rate of 0.035%/m. Just over one-third of this thinning takes place between the undeformed section and the diapir flank, accounting for 19% stratigraphic thinning at a rate of 0.016%/m, and the remaining two-thirds of thinning occurs between the flank and the crest, at a rate of 0.097%/m, accounting for 34% stratigraphic thinning. The thinning rate between the salt flank and the crest is 5.9 times greater than that between the undeformed section and the salt flank (Table 3).

| Diapir growth and roof folding
In our model, we invoke a constant upward movement, or growth, of the salt diapir. Therefore, our models show how different sedimentation rates can dampen the late-stage growth of diapirs. The greatest upward movement of the diapir is observed in the model with no sedimentation (M1) (Table 2; Figure 7), because in this case there is no roof to resist the upward flux of salt. M1 is taken as the base case. Upward movement (growth) of the diapir is reduced with the addition of sediment ( increases the roof thickness towards the limit where diapir growth can occur (i.e. diapir height approaches <66% of overburden: Schultz-Ela et al., 1993). The amount of growth decreases, compared to the base case (M1), with increasing sedimentation rate (being limited to 70% in M4), showing that sedimentation rate is a key control on the burial of salt topography. This observation is in agreement with existing models (Fuchs et al., 2011;Giles & Lawton, 2002;Hudec & Jackson, 2007;Peel et al., 2020).
We also note that the anticline defining the early diapiric folded roof is widest in M1 and decreases with an increasing sedimentation rate (Table 2; Figures 7 and 9). In models with added sedimentation (M2-M6), the early diapiric roof anticline is 20%-24% narrower than in M1 (no sedimentation). The lack of variability between different sedimentation conditions may imply that sedimentation rate has only a minor control on early diapiric anticline width, with other controls such as salt supply, salt viscosity and regional tectonics (not modelled) being more important (Fuchs et al., 2011;Koyi, 1998). However, the anticline within the syn-kinematic folded roof is widest in M4, appearing to widen with increasing stratigraphic thickness (i.e. sedimentation rate). Synkinematic folded roof thicknesses are more variable across the models, because they are controlled by the sedimentation rate. Increasing the syn-kinematic folded roof width with sedimentation rate supports fold wavelengths being larger for thicker overburdens (Bonini, 2003;Duffy et al., 2018;Hudec & Jackson, 2011).

| Fault distribution and deformation zone extent
Fault distribution and salt withdrawal basin extent are interpreted using the nearest neighbour outputs (Figure 4; Figure  S2). Nearest neighbour outputs highlight the amount of displacement that has occurred during 4.6 Myr, relative to an element's initial neighbours; this is used to highlight discontinuities as a proxy for faults. Fault distributions are broadly similar across all models (Figures 4-6) and are summarised  Table 2 here. Numerous predominantly extensional faults, with variable dip directions, are identified in the early diapiric sequence in all experiments. These faults have displacements of metres to 100s of metres, with the greatest throws being observed between layers 1-7 (Figures 4-6). More contractional structures are observed in M4 and M6 compared to the other models. This corresponds with diapirs experiencing the greatest reduction in growth (Figures 5 and 6; Table 2), and could suggest that localised compressional stresses increase as the diapir is restrained beneath early thick sedimentation due to overpressure . However, this could also be an artefact of subtle differences in early diapiric and syn-kinematic anticlines and variations in salt topography. Steep structures appear to develop over the crest of the growing diapir, but are difficult to decipher in terms of slip style (i.e. normal or reverse) due to the relatively small number of displaced neighbouring elements. These crestal structures extend into the syn-kinematic strata F I G U R E 8 Comparison of the flank section stratigraphy across all models with sediment (M2-M6). See Figure 5 for the location of the flank profile, which is the same across all models. Thicknesses and bedding orientation highlighted. The flank location shows the most deformation in all models, greatest modulation is observed in M2, and least modulation is observed in M4. Generally, modulation decreases upwards, however, in M6 (layers A-J) overlying the overburden anticline (M3-M6); these discontinuities are largest and extend furthest into the syn-kinematic overburden under greater sedimentation rates (compare M4 with M3; Figure 5). In all experiments, layers 8 and 9 are dominated by small-scale faults that are localised to those layers. The edge of the salt withdrawal basin is taken as the distance of the furthest faults from the diapir; outside this zone, the strata are undeformed (Figures 4-6). The salt withdrawal basin, and associated deformation zone, is ca. 2,300 m wide in all models, accounting for three times the initial maximum width of the diapir (Figures 4-6). The similar extent of the salt withdrawal basin across all experiments suggests that syn-kinematic sedimentation only has a minor control on the deformation of early diapiric layers and structural configuration.

| Stratigraphic architecture variability with sedimentation rate
Here, the variability in thinning rate is compared between different models with syn-kinematic sedimentation. Layer A is excluded from descriptions in all models, and thus from our comparisons, as it fills to a linear, instantaneous base level and therefore does not always represent the sedimentation rate of S1. In the absence of layer A, more layers would onlap and less would be laterally extensive across the models as more relief would need to be healed. However, the inclusion of a consistently thick layer A across all models enables easier comparison and could represents a pelagic drape that partially heals topography prior to the onset of a coarser-grained depositional system (e.g. a turbidity current-fed channel and/ or lobe). Therefore, 11 layers (B-L) are described and compared in Table 3.
Under slow sedimentation rate (M2) all 11 layers of onlap topography and are not laterally extensive. Under intermediate sedimentation rate (M3) two layers onlap and nine are laterally extensive. In fast sedimentation (M4) and decreasing sedimentation (M6) this rises to 10 laterally extensive layers and one onlapping layer. Under increasing sedimentation rate (M5), the initial four layers (including three which are deposited under slow sedimentation) onlap topography and the remaining 7 are laterally extensive ( Figure 5). As expected, layers are more laterally extensive under higher sedimentation, suggesting that the effects of halokinetic modulation decrease more rapidly upwards under higher sedimentation rates (Cumberpatch, Finch, et al., 2021;Peel, 2014;Sylvester et al., 2015). The final stratigraphic thickness is greatest under fast sedimentation rate (M4, 570 m), least under slow sedimentation rate (M2, 228 m), and at a similar intermediate level for intermediate (M3, 401 m), increasing (M5, 441 m) and decreasing (M6, 382 m) sedimentation rates, logically showing that net sediment volume is the most important control on sediment thickness (Figures 7-9).
In almost all of the models, the thinning rate is between 5.9 and 6.1 times greater between the salt flank and the crest than it is between the undeformed section and the salt flank (Table 3; Figure 9). This rate is higher (6.6 times) under increasing sedimentation (M5), which suggests that salt structures have greater influence in models where sedimentation rate is initially slow (e.g. S1 in M5; Fuchs et al., 2011;Giles & Lawton, 2002;. Whilst the primary mechanism modulating the stratigraphic architecture is stratigraphic thinning, stratigraphic thickening into basins driven by salt withdrawal at depth must not be disregarded. Localised thickening into salt withdrawal basins is observed in all models. Such thickening accounts for 13% of the stratigraphic thickening in M5 (Cumberpatch, Finch, et al., 2021). This highlights how the presence of a growing diapir can be associated with localised additional accommodation (as well as a reduction of accommodation) due to the at-depth evacuation of salt from the source layer to feed the growing diapir, as is evidenced by subtle thickening into the basin (Figures 5 and  6). Accommodation reduction over the crest of the diapir is driven by diapir growth, which is recorded by stratigraphic thinning. Accommodation increases in salt withdrawal basins are accounted for by stratigraphic thickening and salt migration at depth.
Through time, thickening and thinning are eventually reduced as the halokinetic modulation on stratigraphy is minimised with the burial of the salt-cored topographic high and its flanking depocentres. In all experiments with layers that extend across the entire model, the thinning rate and bedding orientation change up-section (Table 3, Figure 8). In M2-M5, a decrease in thinning rate up-section is observed. Bedding orientations are variable but generally decrease (flatten) upwards in M2-M5. However, under decreasing sedimentation (M6) an overall increase in thinning rate upsection is noted ( Figure 6; Table 3), in addition to a slight increase in bedding dip between layers J, K and L. This suggests that initially thicker layers are less deformed, but that as diapir growth continues and sedimentation decreases, thin layers are still influenced by topography associated with rising salt (Giles & Lawton, 2002; | 2589 EAGE CUMBERPATCH ET Al. Hudec & Jackson, 2007;Soutter et al., 2019;Sylvester et al., 2015). Overall, halokinetic modification reduces with increasing sedimentation rate as halokinetic bathymetry is buried. Typically, such alteration decreases up stratigraphy, and laterally outwards from the diapir in agreement with outcrop and subsurface analogues globally (Figures 9-11; e.g. Banham & Mountney, 2013a, 2013b, 2014Doughty-Jones et al., 2017;Giles & Lawton, 2002;Kernen et al., 2012Kernen et al., , 2020Mayall et al., 2010;Oluboyo et al., 2014;Poprawski et al., 2014Poprawski et al., , 2016Pratson & Ryan, 1994;Ribes et al., 2015;Rodriguez et al., 2020;Wu et al., 2020).

| Comparison to natural examples
A key challenge for numerical models is ground-truthing against natural prototypes (Burgess, 2012;Oreskes et al., 1994). Here, the key findings and predictions from the DEM (Figure 9) are compared to published subsurface (Figures 10 and 11) and field analogues ( Figure 12) to understand their applicability and limitations.  Birch & Haynes, 2003). Pierce's tectonostratigraphic history spans ca. 200 Myr and is summarised as Jurassic reactive-active diapirism, followed by Cretaceous-Cenozoic passive diapirism, and contraction-driven active diapirism during Alpine compression (Scott et al., 2010). Despite a longer-lived and more complex evolution, the Pierce diapirs show geometrical similarities with several of our DEMs (Figures 9 and 10). In the Pierce example, stratigraphy is near horizontal ca. 2 km away from both diapirs, but is upturned adjacent to them, comparable with model results (Figures 5,6,9,and 10). An average seismic velocity of 2,000 m/s was used for approximate depth conversion for the entire time-migrated seismic data to calculate the thinning rate for stratigraphy adjacent to the Pierce field. This is an oversimplification of seismic velocity, which varies with depth and lithology. However, it is suitable to give a broad comparison to our modelled values (Table S2). The generation of brittle deformation throughout the Cenozoic stratigraphy over the crest of both north and south Pierce (Figure 11; Carruthers et al., 2013) corresponds to high zones of relative displacement across the crest of the diapir, in M3-M6, which extend into the syn-kinematic stratigraphy (Figures 5 and 6, Figure  S1). Similar crestal deformation has been demonstrated in physical models (e.g. Davison et al., 1993). Different model outputs may be applicable to different parts of stratigraphy due to the changing ratio of diapir rise rate and sedimentation rate. For example, stratigraphic architectures comparable to M2, M3 and M5 are observed in different stratigraphic packages around the Pierce diapirs ( Figure 11).
The unique observation from M2 (slow sedimentation) is non-extensive layers that thin towards and eventually onlap and pinch out against the long-lived, salt-related topographic high (Figure 9a). M2 results are analogous to the Top Cretaceous (lime green) to Mid Eocene (red) stratigraphy around the Pierce diapirs (Figure 11c), which is not laterally extensive across the diapirs (Birch & Haynes, 2003;Carruthers et al., 2013;Davison, Alsop, Birch, et al., 2000;Scott et al., 2010). The Top Cretaceous -Mid Eocene package, equivalent to S1-S3 in the models, thins significantly towards both diapirs, before onlapping the flanks (Figures 9a and 10a,c). The amount and rate of thinning reduce through time, from 51% at 0.044%/m in the Paleocene (Top Cretaceous -Top Lista interpretations) to 21.3% at 0.019%/m between in the Eocene (Tay and Eocene interpretations (Table S2)) in agreement with model observations showing thinning rates reducing upwards. These values are similar to the overall thinning of the slow sedimentation model (34%, at 0.029%/m; Table 3). Despite this apparent similarity in thinning rate values, regional sediment volumes are high throughout the Palaeocene (10,000 km 3 /my; Liu & Galloway, 1997). Specific lithostratigraphic units such as the Forties Sandstone are associated with ca. 200 m of sandstone deposited in ca. 1 million years (Eldrett et al., 2015;Kilhams et al., 2014). Sedimentation rates for the UK Palaeocene stratigraphy are uncertain, and spatially and temporarily variable around the Pierce diapirs, with rates ranging from 0.085 to 0.4 mm/yr (Eldrett et al., 2015;Kilhams et al., 2014;Liu & Galloway, 1997). The lower values are in broad agreement with our slow sedimentation input values. The upper end, however, indicates intermediate -fast sedimentation rates in our models, which have probably been coeval with very high diapir rise rates to form an overall geometry typical of slow sedimentation rates (Carruthers et al., 2013). This is further evidenced by the steep upturning of stratigraphy adjacent to the diapirs (Figure 11c; Giles & Lawton, 2002;Hudec & Jackson, 2007). In our model, we isolate and vary sedimentation rate, whilst in reality the dynamic ratio of sedimentation rate and diapir rise rate control stratigraphic architecture. It is suggested that during the Paleocene -Eocene sedimentation rates were high, but diapir rise rates were higher, likely driven by sediment loading (Carruthers et al., 2013) creating a relative effect of a 'low sedimentation rate' akin to our M2. In M3 (intermediate sedimentation rates), layers are initially non-extensive (S1), displaying early onlap and thinning (26.4% at 0.023%/m; Table 3; Figures 9b and 10b) similar to Eocene (orange) to Oligocene (peach) stratigraphy adjacent to the north Pierce diapir (Figure 11), which thin by 16% at 0.014%/m. Subsequent modelled layers (S2) are extensive but thin towards the diapir high (60% at 0.04%/m; Table 3), analogous to the Oligocene (peach) to Mid-Miocene (light peach) across north Pierce (65% at 0.044%/m; Table S2; Figure 11c; Carruthers et al., 2013). In our simulations, 65% of the total thinning occurs between the flank and crestal locations, similar to subsurface observations of 85% of the total thinning around north Pierce occurring between the flank and the crest. After this initial modulation (S1, S2), subsequent M3 stratigraphy (S3) records a reduction in thinning rates and halokinetic influence on stratigraphy upwards (reduction from 60% to 23% total thinning: Table 3; Figures 9b and 10b). This is also observed, albeit at a less extreme rate, between the Middle Miocene and Pliocene interval around the Pierce diapirs (reduction from 65% to 49% total thinning; Table S2; Figure 11). Much of the excess thinning in the Pierce example, compared to the model, occurs between the undeformed and flank position (28% of thinning, compared to 6% in M3; Tables 3 and Table  S2). This is likely due to Cenozoic compressional forces driving diapir growth and upturn of stratigraphy (Birch & Haynes, 2003;Carruthers et al., 2013;Scott et al., 2010), resulting in a less gradual reduction in halokinetic alteration upwards with respect to M3.
Based on this comparison of the lateral extent of layers, it is possible to infer that these examples represent intermediate sedimentation rates, relative to diapir rise rate. Sedimentation volumes for the Eocene to Oligocene are ca. 4,000 km 3 /Myr, lower than for the Top Cretaceous-Eocene. This suggests that the overall diapir rise rate had been slower during the Eocene to Oligocene, with respect to the Top Cretaceous-Eocene time, giving a stratigraphic architecture typical of lower diapir rise rates relative to sedimentation rates (Carruthers et al., 2013).
In M5 (increasing sedimentation rate), S1 and early S2 are initially isolated on either side of the diapir in saltwithdrawal basins and onlap salt-cored topography (Table 3). Late S2 stratigraphy is extensive across the model, thinning towards the topographic high, with modulation decreasing up-stratigraphy (S3, Figure 9d). When combining the Cretaceous to Eocene (M2 analogue, Figure 11) and Eocene to Mid-Miocene (M3 analogue) stratigraphy adjacent to the Pierce diapirs, we observe an increasing-upwards sedimentation rate, relative to diapir rise rate (Den Birch & Haynes, 2003;Hartog Jager et al., 1993;Jennette et al., 2000;Kilhams et al., 2014). The Cretaceous to Eocene was deposited when sedimentation rate <diapir rise rate, leading to the isolation of salt withdrawal basins and onlapping of the diapir flanks (S1, M5, Figure 11d; Carruthers et al., 2013). The Eocene -Mid Miocene (S2 equivalents) were likely deposited when sedimentation rates were higher relative to diapir rise rates, or when salt supply was equal to sedimentation rate (S2, M5, Figure 11d; Carruthers et al., 2013). Subsequent layers (Mid-Miocene -Pleistocene) are broadly extensive across the diapirs, similar to S3 in M5. The sedimentation volume for this interval is fairly low (ca. 2,000 km 3 /Myr; Liu & Galloway, 1997), so this observed reduction in halokinetic influence upwards is likely driven by a slowing or cessation of diapir growth rather than being driven F I G U R E 1 0 Tectonic framework of the North Sea rift system, and structural map of the Central Graben showing the location of salt diapirs related to major basin faults and Jurassic salt withdrawal basins (Carruthers et al., 2013). Red box locates subsurface analogue used for comparison to models and green line locates Figure 11  purely by high sedimentation rates (Carruthers et al., 2013). Pliocene to Pleistocene stratigraphy is extensive across north Pierce (S2 and S3, M5), but only the Pleistocene is present extensively across south Pierce, due to differential diapir growth histories and cessations (Carruthers et al., 2013;Scott et al., 2010). The upwards reduction in sedimentary thinning observed in the increasing sedimentation model (61% in S2 to 22% in S3) is similar to that observed between the Eocene to Mid-Miocene (88%) and the Mid-Miocene to Pleistocene (18%) in the subsurface example, highlighting the potential applicability of the results from our models. Whereas these subsurface thinning rates are similar to those in our models, it is important to note that our models do not take into account erosion, or periods of non-deposition once sedimentation begins at 2.2 Myr. Therefore, in nature, 'apparent thinning' could be derived from periods of erosion, removing stratigraphy and generating halokinetic unconformities. This is particularly noticeable at the Mid-Miocene unconformity where 'thinning' is at least partly accommodated by Cenozoic compression rejuvenating diapir growth and increasing diapir rise rate, causing erosion of the diapirs overburden, without a variation in sedimentation rate (Carruthers et al., 2013). Our models generate realistic halokinetic unconformities, with variable bedding orientations between the early diapiric sequence and the syn-kinematic sequence due to a period of uplift and non-deposition prior to the commencement of sedimentation ( Figure 8).
In nature, halokinetic sequence architecture is controlled by the dynamic ratio between sedimentation rate and diapir rise rate (e.g. , F I G U R E 1 1 Subsurface example of model application comparing some of the modelled results to stratigraphy from the Pierce Field, Eastern Central Graben, UK North Sea. (a, b and d) Models 2, 3 and 5 are analogous to different parts of the stratigraphy around the Pierce diapirs. (c) Interpreted time-migrated seismic profile across the Pierce diapirs. S1, S2, S3 highlight the likely sequences for comparison to models. The colour of the text represents which model those packages could represent (e.g. red represents M5). The location of the undeformed, flank and crest stratigraphic locations used for thinning rate calculations are shown for north Pierce and are the same spacing as those used for model calculations ( Figure 5; Table 3). Thinning values are calculated in Table S2. Seismic data courtesy of PGS (MegaSurvey Plus 3D seismic data) from Cumberpatch, Finch, et al. (2021)  such that an 'apparent' increase in sedimentation (reducing halokinetic influence upwards) could represent a slowing of diapir growth due to regional tectonic quiescence or depletion of the salt source layer. Our observations and comparison to DEMs with variable sedimentation rates are consistent with diapir rise to sedimentation rate ratios derived from halokinetic sequence studies (Carruthers et al., 2013). The role halokinesis plays in shaping the stratigraphic architecture of the north and south Pierce diapirs is reduced from the Oligocene and Pliocene respectively, due to source layer depletion, resulting in halokinetic bathymetry being gradually buried (Figure 11; Birch & Haynes, 2003). Integration of, and comparison between, DEM and subsurface data, demonstrates the importance of understanding local (salt layer variations) and regional (tectonic and sedimentation rate) controls when disentangling salt basin evolution.

| Comparison to outcrop: Stratigraphy around the Bakio diapir, northern Spain
The model simulations presented here document stratigraphic architectures and structural deformation, but do not invoke a specific sedimentary environment. Different compositions (e.g. carbonate or siliciclastic) and sedimentary environments (e.g. fluvial or deep-water) will respond differently to salt influence (Adams & Kenter, 2012). In our models, slow sedimentation (M2) leads to the natural exposure of the early diapiric sequence above the diapir in the absence of a syn-kinematic cover or roof ( Figure 5). This is likely to be common in marginal marine or terrestrial environments, where field data indicate the crest is often exposed, eroded and reworked into the syn-kinematic succession (Banham & Mountney, 2013a, 2013b, 2014Counts & Amos, 2016;Counts et al., 2019;Ribes et al., 2015Ribes et al., , 2017. If accommodation is continually available above diapir, as is the case in deep-water environments, a pelagic drape will likely cap the diapir, whereas coarser-grained deposition will be restricted to the flanks Poprawski et al., 2014Poprawski et al., , 2016Rowan et al., 2003). Our models are applicable to all depositional environments, but fundamental controls on depositional architecture (e.g. accommodation, sedimentary processes) must not be overlooked.
Below, we compare the observations of M5 (increasing sedimentation) to a well-exposed halokinetically influenced deep-water succession and describe how integrating DEM results with outcropping examples can help reduce uncertainty in subsurface energy reservoir prediction.
Exposed Aptian-Cenomanian strata adjacent to the Bakio Diapir, northern Spain document an increase in sedimentation rate associated with an increase in erosion of the hinterland (Figure 12; García-Mondéjar, 1996;Martín-Chivelet et al., 2002;Puelles et al., 2014), similar to M5 increasing from 'slow' to 'fast' relative sedimentation rates. The resulting progradational pattern manifests in the deepwater succession as an upward change from thin-bedded low-density turbidites, deposited in the lobe fringe, (equivalent to S1) to thick-bedded high-density turbidites deposited in a lobe axis (equivalent to S3) . Early diapiric (the Urgonian Group: Figure 12) and syn-kinematic (the Black Flysch Group) depositional elements are deformed closest to the diapir, and deformation intensity decreases away from the diapir, being minimal outside the halokinetically influenced sequence (ca. 700 m wide either side at Bakio); this observation is consistent with the predictions of our model (Figures 6 and 12). Outside the zone of halokinetic deformation, strata look similar to those deposited in an unconfined, salt-free setting, in terms of their architecture and facies distributions (e.g. Prélat et al., 2009: Figure 12d).
The Aptian-Albian was initially isolated on either side of the diapir, representing relatively slow sediment accumulation rates with respect to diapir growth. Individual thin beds representing distal deposition thin (by up to 1%/m; Figure 12e) and pinch out towards the salt-cored topographic high, in accordance with layers A-D in the increasing sedimentation model (S1, Figure 12) or the entirety of the slow sedimentation model (M2) (Cumberpatch, Finch, et al., 2021;. At Bakio, different sediment routing systems develop when stratigraphy is restricted to either side of the diapir due to different controls on deposition . Under increasing sedimentation rates, the Albian-Cenomanian (Black Flysch Group: Figure 12) stratigraphy is laterally extensive, representing 'intermediate' sedimentation conditions in our model (S2 in M5 and the entirety of M3). However, rather than simply representing high sedimentation rates, as might be inferred from our numerical models, this stratal architecture likely reflects both an increased sediment supply from the Landes Massif (García-Mondéjar, 1996;Martín-Chivelet et al., 2002;Puelles et al., 2014) and a reduction in salt rise rate (Poprawski et al., 2014(Poprawski et al., , 2016Roca et al., 2021). In the field, this stratigraphy shows a reduction in the numbers of mass transport deposits (MTD) upwards, which cannot be resolved in the models, but can be compared to the decrease in stratigraphic dip upwards (from 12° to 2° between layer E and I, Figure 8), reducing instability. The lower Black Flysch stratigraphy shows a reduction in halokinetic deformation and thinning rate upwards (S2) that is comparable to modelled results, where thinning rates decrease from 0.053% to 0.006%/m through time. Eventual deposition of the upper Black Flysch Group, across the diapir, akin to relatively high sedimentation rates in our model (S3 in M5 or the entirety of 2594 | EAGE CUMBERPATCH ET Al. F I G U R E 1 2 Conceptual facies diagram for a deep-water succession based on the integration of field-based facies analysis around the Bakio diapir, Basque Cantabrian Basin, Northern Spain (see  for summary), overlain onto the result of M5. (a) Interpretation for deep-water facies (based on field facies analysis) overlain on the upper part of the increasing sedimentation model (located in Figure 6). Syn-kinematic layers from the model are grouped (S1-S3) and coloured, and depositional elements derived from field studies are overlain to show facies variability. H and N are theoretical stratigraphic profiles in the halokinetically influenced and non-halokinetically influenced zones, respectively. (b) Location map of the field analogue, for a full geological discussion, see .
M4) heals the remnant diapiric topography, and shows almost no halokinetic influence (in similarity to layers I-L in the increasing sedimentation model), except for its depositional location within remnant diapiric topography (Cumberpatch, Finch, et al., 2021). The relative increase in sedimentation rate at Bakio is driven primarily by an increase in sediment supply due to the uplifting source area, which is comparable to M5. However, a coeval reduction of salt supply due to welding between the Bakio and Guernica diapirs cannot be ruled out Poprawski et al., 2014Poprawski et al., , 2016Roca et al., 2021).
In the syn-kinematic stratigraphy, particularly those equivalent to S2 in M5, thick-bedded sandstones deposited in channels and lobes dominate topographic lows, where sedimentary flows were focussed around salt-cored topographic relief (Figure 12d; Doughty- Jones et al., 2017Jones et al., , 2019Mayall et al., 2010;Rodriguez et al., 2020;Sylvester et al., 2012). Towards the highs, the lower-density part of the F I G U R E 1 3 Halokinetic zonation scheme shown for M3. The model is divided into five zones, four of which experience some form of halokinetic influence. Minimal deformation zones 1 and all stratigraphy outside of it show virtually no modulation by salt diapirism. Halokinetic influence increases towards the diapir-cored high, and changes from minor thickness changes in the 'withdrawal basin' zone to onlap and abrupt pinch out in the 'onlap' and 'salt flank' zones. Thinned and reduced stratigraphy are observed over the diapirs crest. The table highlights that while stratigraphic trap quality may be greater closer to the diapir, reservoir thickness and quality are likely higher further from the diapir, showing a 'trade-off' exists for subsurface energy exploration and production. Similar zonation is possible for all models [Colour figure can be viewed at wileyonlinelibrary.com] 2596 | EAGE CUMBERPATCH ET Al. flows may run-up topography depositing thinly bedded muddier sandstones towards the pinch out (Figure 12e; Kneller & McCaffrey, 1999;Soutter et al., 2019). At Bakio, crestal deposition is assumed to be limited to a thin pelagic drape during S1 and much of S2 due to elevation, which is comparable to M5, often leading to remobilisation and the formation of MTDs. During S3 halokinetic bathymetry is healed, and crestal deposition is extensive, with minimal modulation (the upper Black Flysch Group).

| Implications for subsurface energy
Despite advances in extent and resolution of 3D seismic reflection imaging, the salt-sediment interface remains difficult to image, yet determining its position and precise geometry is crucial in helping to appraise stratigraphic-structural traps for hydrocarbon, carbon storage, and geothermal prospects globally (Jones & Davison, 2014;Warren, 2016). Utilising outcrop analogues ( Figure 12) can help provide sub-seismic scale depositional facies information, helping reduce uncertainty in reservoir quality and distribution. Numerical modelling results do not represent specific analogue conditions nor a 'snapshot' in time, and they can, therefore, help to quickly identify generic depositional architectures, deformation patterns and sediment thickness relationships as a function of several forcing parameters, such as variations in sedimentation rates.
Using stratigraphic architectures from our DEM and sedimentological data from field examples, we can improve predictions of the likely architecture of syn-kinematic stratigraphy and sedimentology around salt structures, which are poorly imaged in seismic reflection data. Understanding this requires a multi-scalar and multi-technique approach. For example, models provide details about gross thickness changes and geometry, whereas field analogues enable inferences about reservoir quality and net-to-gross.

| Depositional reservoir quality
Regardless of the amount of sedimentation, our models show that stratigraphy thins as it approaches the diapir, suggesting a reduction in the amount of total reservoir close to the structure (Figures 8 and 13; . Siliciclastic depositional environments show a thinning of sandstone towards the topographic high and an overall concentration of high reservoir quality units at the base of topography (Figure 12a,d), such that a salt-related combined structural-stratigraphic trapping mechanism becomes unlikely ( Figure 13; Kane et al., 2012;Stricker et al., 2018). Muddier (lower reservoir quality) and thinner (lower netto-gross) units are expected closer to the diapir (Figures 8,   10e and 11; Banham & Mountney, 2013a;. These units are more likely to be overpressured due to upward rotation, creating a large pressure head, with the top seal rocks unable to hold back a significant hydrocarbon column (Figure 13; Heidari et al., 2017Heidari et al., , 2019Nikolinakou et al., 2014aNikolinakou et al., , 2014bNikolinakou et al., , 2018. In carbonate environments, shallow platform or reef facies with excellent reservoir potential may be present over salt highs (e.g. Burgess et al., 2013;Riding, 2002;Teixell et al., 2017). As salt growth continues, fractures are generated in the overburden, which could form significant secondary porosity within the carbonate reservoirs, increasing quality and producibility (He et al., 2014;Howarth & Alves, 2016;Saura et al., 2016).
Supplementing subsurface data with modelled stratal architectures and depositional facies observations from exhumed halokinetically influenced settings globally (e.g. Banham & Mountney, 2013a, 2013b, 2014Counts & Amos, 2016;Counts et al., 2019;Poprawski et al., 2014;Ribes et al., 2015Ribes et al., , 2017 is recommended as a useful workflow for building reservoir models for salt basins with limited data. Observations from multiple scales and types of models can be combined to further reduce the uncertainty associated with reservoir quality prediction, for example, recent FEM has shown porosity is lower than expected near the vertical parts of salt structures and higher than expected at the base of diapirs, due to mean principal stress variations (Nikolinakou et al., 2014a).

| Halokinetic zonation
The model results show that a deformation zone exists on either side of the diapir in all experiments (Figures 4-6). Outside of this zone, the syn-kinematic and early diapiric stratigraphy are undeformed ( Figure 13). The extent of this salt withdrawal basin is 1,150 m on either side of the diapir (2,300 m in total). Therefore, the total zone of halokinetic influence in all models is approximately three times the original diapirs' maximum width (750 m), with a deformed zone of 1.5 diapir widths on either side of the structure. The width of the deformation zones is comparable across all models (Figure 7), and therefore it is shown that sedimentation rate is unlikely to have a significant control on the width of the zone of halokinetic influence Hearon et al., 2014). Other factors such as salt supply, salt viscosity, and style and magnitude of regional tectonics (which are not modelled) will in nature, influence the width of the halokinetically deformed zone (Fuchs et al., 2011;Koyi, 1998). The model can be further divided into zones based on onlap geometry and thinning rates, which highlight the 'trade off' between reservoir thickness and stratigraphic trap potential in subsurface plays ( Figure 13). In the flank locations, bedding dips and thinning rates are shown to be greater under | 2597 EAGE CUMBERPATCH ET Al. slower sedimentation rates compared to higher sedimentation rates (Figure 8; Table 3), which is important when predicting hydrocarbon column height. Significant overpressures on reservoirs below can be created by fast sedimentation rates (Figures 7 and 8; Peeters et al., 2018).

| Fault distribution
Radial faults associated with salt diapirs are shown to cause the compartmentalisation of reservoirs (Birch & Haynes, 2003;Charles & Ryzhikov, 2015;Coleman et al., 2018;Peeters et al., 2018;Scott et al., 2010). Our DEM replicates localised fault growth, evolution and propagation because the contacts between elements are treated as potential displacement surfaces. Our models all document a similar fault pattern, suggesting that faulting, and therefore fault compartmentalisation is not heavily influenced by sedimentation rate. As well as seismically resolvable faults, outcrop and borehole data indicate brittle deformation is significant in salt basins (e.g. Koestler & Ehrmann, 1991). Our DEM replicates this brittle deformation (Figures 4-6); 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 compartmentalisation and seal integrity in the subsurface. Faults, when sealing, could act as lateral permeability barriers, especially if the faults and surrounding reservoir rocks become cemented with salt and salt-related breccia (Li & Urai, 2016;Van Bergen & de Leeuw, 2001). DEM is, therefore, advantageous in predicting potential traps, conduits and baffles, due to its replication of diapir-related brittle deformation.

| Future work
The DEM presented here is useful for predicting regional trends and studying generic interactions of salt structures and stratigraphy but further work is required to recreate specific complicated salt geometries. Other suggestions for further development of this model are to incorporate reactivation and cessation of halokinesis, studying the impact on stratigraphic architectures and the development of halokinetic unconformities due to periods of diapiric rejuvenation, uplift, erosion and non-deposition. Incorporating erosion and remobilisation into the model would also more realistically represent natural settings where entrainment and failure influence thickness patterns, stratigraphic architectures and can remove the diapiric roof. Other investigations may include assessing the influence of evolving salt geometries relative to salt thickness. In nature, salt diapirs are rarely isolated structures, and therefore subsequent studies will focus on the interaction of multiple salt structures, with different growth histories on different sedimentary successions. Many of the limitations of the model described above, are due to its two-dimensional nature. The ultimate aim of future work is to develop a three-dimensional DEM to better understand the four-dimensional variability in halokinetically influenced stratigraphy and associated subsurface energy reservoirs. 6 | CONCLUSIONS 1. DEM can form an integral part of the workflow when studying salt-sediment interactions. Here, a DEM is used to study the variability in stratigraphic architecture and deformation patterns around a growing salt structure under different sedimentation rates. 2. The models generate realistic salt-related faults. In all models, structural deformation and extent of halokinetic influence are similar, and syn-kinematic strata, at least initially, are isolated to either side of the diapir, thinning and onlapping towards the high. 3. Under slow sedimentation rate (M2) deposition is restricted to salt withdrawal basins on either side of the diapir throughout evolution, while in M3-M6 sedimentation eventually occurs over the diapir crest, often associated with significant lateral thickness variability. Diapir growth is most inhibited under fast aggradation (M4) and the halokinetic influence on stratigraphy reduces quickly with time. 4. Thinning of syn-kinematic stratigraphy from the undeformed section to the diapir flanks is greatest under slow aggradation (M2, 34%), and least under fast aggradation (M4, 12%). In all models, thinning is about six times greater between the salt flank and crest, compared to the undeformed section and the salt flank, indicating more intense deformation close to the diapir. 5. Thinning rate decreases through time (up stratigraphy), showing a reduction of halokinetic modulation with increased sediment thickness, as halokinetic bathymetry is 'healed'. This is true for all models except for decreasing sedimentation (M6), which experiences a slight increase in stratigraphy. 6. Our simplified two-dimensional models provide useful analogues for salt-influenced basins with complicated four-dimensional evolutions. Natural examples record the interplay between sedimentation rate and diapir rise rate, whilst models isolate and vary sedimentation rate. Comparison to the Pierce field diapirs, North Sea, shows how different models can be applicable to different parts of stratigraphy and suggest how interpreters could infer likely sedimentation rates and conditions from subsurface stratigraphic geometries.

|
7. Facies, and thus reservoir distribution, around salt diapirs will differ from unconfined settings due to halokinetic modulation, both vertically and laterally. A deep-water analogue from stratigraphy adjacent to the Bakio diapir shows that halokinetically influenced facies (e.g. the salt flank) contain thin beds, sandstone pinch outs and increased mass transport deposits in comparison with the allogenic deposits (e.g. the undeformed zone), which are difficult to decipher from deep-water strata in non-salt influenced settings. 8. Integrating DEM with subsurface and outcrop data helps to reduce reservoir and trap uncertainty in subsurface energy exploration and development.