Overburden deformation induced by dyke-fed conical sandstone 1 intrusions: insights from numerical experiments

6 Conical sandstone intrusions, as a distinct type of hydrocarbon reservoirs and carbon sequestration 7 sites, remain poorly understood regarding their emplacement mechanics. Here, we report a 8 numerical modelling study of conical sandstone intrusions using the two-dimensional discrete 9 element method. We built simplified numerical models that contain bonded elastic particles with 10 predefined mechanical properties in an open box as the overburden, and a thin tube filled with 11 unbonded particles as the feeder dyke that is connected to the upper box. The dynamic behavior 12 of the assembly was enabled by displacing of the driving wall that defines the lower boundary of 13 the dyke. The results show that the model composed of soft materials produced a pair of conical, 14 opening-mode fractures in the host sediments as the result of tensile stress concentrations in the 15 fracture tip zones. The overburden deformation was largely localised within the sediments adjacent 16 to the sandbody, without formation of a forced fold and significant uplift of the surface. Differently, 17 the model composed of stiffer materials produced conical fractures that have closed lower 18 segments and opening upper segments with a reverse sense of shear. The intrusion also caused a 19 forced fold in the overburden, with a vertical opening-mode fracture generated in the fold hinge. 20 The modelling results demonstrate that dyke-fed sand intrusions can significantly distort the local stress field, and the overburden can be subjected to fracturing and/or folding due to differential compaction during gradual inflation of the intrusive sandbody. Moreover, deformation patterns of the overburden in response to sandstone intrusion largely depend on mechanical properties of the host sediments.

The discrete particles displace independently from one another, and interact only at contacts 97 between the particles. A soft contact approach is used for particle contact, and the particles are 98 allowed to overlap one another at the contacts. The magnitude of particle overlap is determined by 99 the contact force via the force-displacement law. The mechanical behavior of such particles can 100 be then characterised by the movement of each particle and inter-particle forces acting on the 101 particle contacts. The relationship between the particle motion and it driving forces is provided by 102 the Newton's laws of motion. In addition, the particle contacts can be bonded together such that, 103 the particles act as linear-elastic springs in compression, and cohesive bonding that act in both 104 shear and tension (Fig. 2). The contact bond allows tensile stresses to develop at a contact when 105 there is no overlap between neighboring particles. The contact bond is broken when the inter- Deformation of a bonded aggregate of particles results from the movement of elastic, frictional 112 walls as the confining boundary for the particles. During deformation, particle interactions are seen 113 as a dynamic process with states of equilibrium developed whenever the internal forces reach a 114 balance. The dynamic process is represented by an explicit timestepping algorithm, which consists 115 of repeated applications of the law of motion to each particle, a force-displacement law to each 116 particle contact, and a constant updating of wall positions. At all times, the forces acting on any 117 particle depend exclusively on its interaction with the contacting particles.

119
The DEM method has been used to solve a wide range of problems in granular mechanics. For 120 structural geology and tectonics, the DEM method has been effectively applied to simulate the 121 deformation of upper crustal rock materials, including the development of normal fault (Schöpfer 122 et al., 2006(Schöpfer 122 et al., , 2007a(Schöpfer 122 et al., , 2007b(Schöpfer 122 et al., , 2017Hardy, 2011Hardy, , 2013Smart et al., 2011;Smart and Ferrill, 2018), 123 thrust fault (Strayer and Suppe, 2002;Naylor and Sinclair, 2007;Dean et al., 2013;Morgan, 2015), 124 fracture (Spence and Finch, 2014;Virgo et al., 2014Virgo et al., , 2016, strike-slip fault (Imber et al., 2004;125 Liu and Konietzky, 2018), detachment fold (Hardy and Finch, 2005;Vidal-Royo et al., 2011;126 Meng and Hodgetts, 2019), and fault-related fold (Finch et al., 2003;Hardy and Finch, 2006;127 Benesh et al., 2007;Hardy, 2018). These successful previous applications highlight the 128 appropriateness of the DEM method to investigate the development of faults, fractures and folds 129 associated with sandstone intrusions for this study. and a floor that contained 20,335 particles as the overburden sediments. The particle radii range 135 from 15 to 25 m, following a uniform distribution, which can help avoid hexagonal close packing 136 of these particles. A vertical tube with a width of 100 m, defined by two side walls and a horizontal 137 wall at bottom, Wd, is located underneath the box and connected to the central part of the floor.

138
Such a design satisfies the generally acknowledged model, in which a relatively thin, vertical dyke 139 acts as the fluid migration pathway and also the feeder of sands for conical sandstone intrusions 140 (Huuse et al., 2004;Cartwright et al., 2008). This model has been adapted for several analogue 141 and numerical modelling studies of sandstone and igneous intrusions (Mathieu et al., 2008;142 Galland et al., 2009;Bureau et al., 2014;Gorczyk and Vogt, 2018).

144
The particles in the feeder dyke has a radius range from 5 to 10 m, which also follows a uniform 145 distribution. The smaller size for particles in the dyke allows more particle contacts and a more 146 efficient stress transmission within the system. The particles in the overburden were assigned with 147 a density of 2600 kg/m3. The density for particles in the feeder dyke was assigned a value of 1500 148 kg/m3 to simulate the mixture of water and sand grains. The whole model was gravitationally 149 loaded by 1 g.

151
Unlike continuum methods, the discrete element method prevents one from directly ascribing the 152 macroscopic mechanical properties and the desired aggregate characteristics, due to the particle-153 based nature of the models (Benesh et al., 2007). Instead, we first specified the microscopic 154 parameters to individual particles and their contacts, and then iteratively varied these parameters 155 until the desired macroscopic behavior and characteristics are achieved. The appropriate values of 156 particle stiffness, friction and bond strength were then attained by conducting the numerical equivalent of a biaxial rock mechanics test (Cundall and Strack, 1999), which helped derive the 158 macroscopic mechanical properties (Fig. 3). Trials runs of the rock tests were monitored to 159 evaluate whether the particle assembly showed any nonphysical behavior. In such tests, the 160 synthetic rock sample was loaded in a strain-controlled fashion by displacing the top and bottom 161 walls at a sufficiently slow rate, so as to attain a quasistatic solution. The stresses and strains 162 experienced by the rock sample were determined in a macro-fashion by summing the forces acting 163 upon walls and tracking the relative distance between the walls. The samples were loaded until the 164 axial stress falls below 70% of the peak stress. Following these tests, we selected a particle stiffness  (Table 1), respectively, to reasonably represent soft 169 and stronger clay sediments (e.g. with a higher content of silts and sands) that commonly appear 170 to act as host for conical sandstone intrusions (Huuse et al., 2004). We chose a particle stiffness of 171 1 x 10 7 N/m for the sand particles, both zero for the particle bond strength and friction coefficient 172 to simulate the non-cohesive, fluidized sand grains.

174
The particles were allowed to settle to the bottom of the model and compact under their own weight.

175
Firstly, the particles in the feeder dyke were packed under the force of gravity. Once the mean 176 unbalanced force in the whole particle assembly dropped to a negligible value that indicates the 177 achievement of static equilibrium, the extra particles above the top of the dyke were removed. The 178 trim of the assembly resulted in a small amount of vertical elastic rebound, which elevated the 179 upper surface of the assembly. When the new equilibrium had been achieved, the trimming process was repeated. This operation was iterated until the assembly reached to the dyke top and was 181 thereby considered to be settled. The particles in the overburden were packed in the upper box, 182 following the same routine. At a critical point when the assembly was trimmed to 3 km high and 183 no more than five particles could be removed at a new equilibrium, the packing process was 184 considered finished.

186
Colors were applied to the particles in the overburden to produce visible layering for later bedding 187 correlations, however, the assembly was mechanically homogeneous. The driving wall that define  The stress field witnessed a strong stress perturbation in the overburden around the intrusive body 201 from a previously isotropic stress field prior to sand intrusion, which is characterised by tensile 202 stress trajectories aligned in a half circular manner in the outer zones of the sand body (Fig. 5).
Then, a pair of opening-mode fractures simultaneously nucleated in the layer above the 204 propagation front of the sand body at T1. The right fracture propagated upward at a faster rate than 205 the left fracture, and firstly reached the magenta layer at T2. Both fractures propagated by the 206 coalescence of neighboring en echelon fractures that formed sequentially in the areas above the 207 fracture tips. This was achieved by tensile stress concentrations within fracture tip regions that led 208 to continuous fracture propagation. Later on, the sand body gradually splitted into two branches, 209 with upward-tapering tips pointing towards the fractures. This was followed by sand particles 210 entering opening fracture channels. Notably, the fractures occur in a hybrid mode, exhibiting a 211 reverse sense of shear. The orange layer right above the intrusive sands was uplifted to a higher 212 level, whilst the overlying layers were much less affected. Both the models presented produced conical fractures that exhibit a distinctive conical geometry 236 with a well-defined apex, and are compatible with those observed in seismic profiles (Fig. 1a, 1b).  (Fig. 6).

246
The fractures produced in both models exhibit a hybrid-mode with a reverse sense of shear, which   If the host sediments predominantly consist of soft clays, the dyke, through which fludised sands 306 are transported from the parent body, will tend to grow radially at a critical point when the fluid 307 pressure drops to be inefficient to drive hydro-fracturing (Fig. 8a). The accumulation of sands 308 leads to inflation of the intrusive sandbody and will result in the formation of conical opening 309 fractures at the propagation front of the sandbody due to differential compaction of the encasing 310 sediments. The conical fractures will propagate upward as sandbody inflation continues. The sand 311 grains will penetrate the clays between the fracture channel and the sandbody, and subsequently 312 fill the entire fracture. It is because that the conical fractures are tensile, the sands are readily stored 313 in the fractures when fluid pressure becomes reduced. Under this condition, surface deformation 314 and forced folding can be rather limited.
sandstone intrusion can more likely cause forced folding of the entire overburden (Fig. 8b).

318
Opening fractures can occur above the neutral surface and reach the surface, creating clear scarps.

319
Sub-vertical tensile fractures can be generated due to stretching of the flexed overburden. Notably, 320 when remoblised sands enter the fractures, they will be only preferentially stored in the upper 321 segments of the fractures after fluid pressure drops, because areas of the lower segments of the 322 conical fractures are dominated by compressive stresses and will be likely to be closed.

324
In both cases, fludised sand grains will exploit the conical fractures as a mechanically easier option.

335
This study utilized the discrete element method to simulate dyke-fed sandstone intrusions and their 336 associated overburden deformations. We conclude the following:

337
(1) Differential compaction in the encasing sediments of intrusive sandbody can lead to conical 338 fractures in the overburden. above the intrusive body, whilst stiff host sediments favours forced folding the surface, formation 341 of closed shear fractures below the neutrual surface, opening hybrid fractues above the neutral 342 surface and pure tensile fractures in the hinge zone of the forced fold.

343
(3) The opening segments of conical fractures that are dominated by tensile stresses, can serve as 344 preferential storage sites for the fludised sands.          (b) Schematic model showing the development of conical sandstone intrusions in stiff sediments. State 1, upward propagation of a vertical dyke as a mode 1 hydraulic fracture; Stage 2, shear failure induced by differential compaction due to inflation of the intrusive sandbody; Stage 3, nucleation of opening fractures above the neutral surface of the force fold; Stage 4, conical fractures reached the surface; Stage 5, sand transport along the fractures; Stage 6, sand storage in the upper segments of the conical fractures.