Fracture formation due to differential compaction under glacial load: a poro-elastoplastic simulation of the Hugin Fracture

The Hugin Fracture, discovered in 2011, is an approximately 3.5 km long seafloor fracture in the North Sea. This fracture was unexpected and, due to the geology in the North Sea no obvious explanation could be found. In our study, we adopt the hypothesis that the Hugin fracture was formed by differential compaction controlled by glacial load. We construct a simplified 2D geomechanical model partly covered by top load (ice sheet) and test this hypothesis. We employ transient poro-elastoplastic simulation with a finite element method. For the simulations, we had to make assumptions regarding the material properties, because the fracture is located in-between well locations. We used descriptions from drilling site survey reports and literature values and performed seismic matching form well paths to the Hugin Fracture. Nearby well data were only partly useful due to incomplete logging in the first 400 m below seafloor. To overcome this problem, we introduced a mixing k-value which allows us to easily change the material properties from pure clay to sand. Changing the mixing k-value for each simulation provided information about the limits and robustness of the simulation results. Simulation results show isotropic stress and strain distribution in the horizontally layered, isotropic part of the model that is totally covered by the ice. In the central, channelized part of the model a composite stress and strain pattern develops with sub-vertical focus areas tangential to channel edges. Low stress, strain and deformation values under total load increase drastically soon after the load starts to decrease, resulting in the development of fractures along the focussed zones. Surface deformation such as formation of compaction ridges above stiff clay-filled channels and depression associated with plastic deformation is observed. A fracture and associated surface deformation develop above the shallowest sand-filled channel, very much resembling the observed geometry at the Hugin Fracture. The simulation supports the formation hypothesis for the Hugin Fracture as a compaction fracture and suggests that thin ice sheets may induce differential compaction to a depth of several hundred meters.


Introduction
Recent studies on the Quaternary history of the North Sea (Ottesen et al. 2012;Reinardy et al. 2017) and the sediments in the North Sea Basin, reflect a sustained interest from both academia and industry in better understanding the implications of the glacial history for petroleum exploration and CO 2 -storage (Chadwick et al. 2004). The discovery of active fluid discharge from the seafloor at the east-to-westoriented 3.5 km long Hugin Fracture on the Utsira High in the northern North Sea (Pedersen et al. 2013), shows that unexpected, Pleistocene age features can still be found in apparently well-mapped, highly industrialized offshore areas such as the Norwegian Continental Shelf (Fig. 1).
The geographic area of the Viking Graben and Utsira High is still under investigation as to the ice conditions that prevailed during the Quaternary glaciations (e.g. Graham et al. 2011). Ice evolution, thickness, movement and baseconditions are not yet fully understood.
In general, sediment deposition from glaciers leave poorly sorted low-porosity tills. On the other hand, hydrocarbon charged Quaternary sand depositions with ice-scouring 1 Page 2 of 22 marks in the Southern North Sea demonstrate the possibility of high-porosity interglacial sediment deposition (Haavik and Landrø 2014). Glacial reworking and loading-unloading cycles due to glacial advance and retreat change sediment properties of the overburden of proposed large-scale CO 2 storage formations in the North Sea, Europe (e.g. GCCSI 2015, for locations). Different sediment packages and heterogeneities will undergo different degrees of reworking and compaction that will increase the brittle behaviour of the sediment and may lead to fracturing (Aplin et al. 1999;Bjørlykke and Hoeg 1997;Bjørlykke 2006). Boulton et al. (1999) discuss various subglacial failure geometries down to 150 m subsurface to a Scandinavian ice sheet and how they change with rock properties (sand vs. Fig. 1 a Location of the study area (black rectangle) in block 16 on the Norwegian Continental Shelf (map adapted from The North Sea CO2 storage atlas, Norwegian Petroleum Directorate). The boundary of the Southern Utsira Fm. that is suitable for CO2-storage according to NPDs CO2 storage atlas (NPD 2011) is shown by the red broken line polygon. b Seafloor mosaic from the Hugin Fracture, acquired with highresolution interferometric synthetic aperture sonar (HISAS) at a resolution of 10 cm × 10 cm. Light shades indicate high backscatter from harder sediments and dark shades indicate less backscatter energy from softer sediments granite). They also discuss the effect of a permafrost shield in front of the glacier on failure such as shear failure and hydraulic fracturing. Shear failure is suggested to occur in depths of at least 100-200 m. In unlithified rocks they may occur to depths in excess of 1000 m. Considering the generally poorly lithified rocks in the North Sea Basin (about 900 m thick, Ottesen et al. 2014), their work suggests that shear fracturing can take place in the entire interval. Boulton et al. (1999) find also that hydraulic fracturing is restricted to depths of up to 150 m with vertical fractures occurring at shallow depths (1-10 m), low angle up-glacier dipping subglacial fractures in intermediate depth, and horizontal fractures in the ice marginal zone at 150 m. The vertical hydraulic fractures can be filled with injected rock masses, sand or till in their examples, according to the potential gradient and layering (Boulton et al. 1999). Sejrup et al. (1994) found a more than 10 m thick overconsolidated till with low water content (around 30%) and a carbon-age of about 32.750 kyears in a depth below 4.5 m in the B2001 Sleipner core. Peacock (1995) proved the existence of drastic palaeo-environmental changes in the Viking Bank area spanning from glaciomarine to shallow marine and periodical dry land at the Viking Bank. Local stress variations due to (glaci)-fluvial erosion and sediment accumulations (fluvial fan) may contribute to local failure. Detailed understanding of this process and its effects on Pleistocene sediments in the North Sea Basin is so far inhibited due to insufficient core data. This will likely change in the near future as new data are acquired (e.g. Barrio et al. 2015). One way to evaluate glacier-related processes with a potential to weaken sealing properties of the overburden, is to employ numerical simulations.
Numerical simulation of rock physics and fluid-rock interaction has become easier to facilitate with the increasing computing capacity of modern computers. In addition, algorithms for elastoplastic and poroelastic models have been improved and, recently, poro-elastoplasticity has become more widely feasible and accessible (Nikolinakou et al. 2012;Lee et al. 2015). Hence, modern geomechanical models can more accurately represent the field situation depending on the amount and quality of rock or soil property measurements available to constrain the model. Where distinct measurements are absent, discovered structures can, nonetheless, be better understood through simulations by applying literature properties.
Our hypothesis is that the Hugin fracture formed as a result of differential compaction due to loading and unloading by glaciers. The current paper attempts to test this hypothesis employing deterministic, transient poro-elastoplastic simulation. Our approach is an attempt to falsify the hypothesis. The idea being that the fracture formation hypothesis must be omitted if the model fails to produce robust simulation results in support for the hypothesis. If the model does produce robust simulation results in favour of the hypothesis, this will indicate that the hypothesis is still valid, but it will not prove that it is correct.

Survey area and geological background
The study area is located on the western flank of the Utsira High in the Norwegian Sector of the central North Sea basin (Fig. 1a) and is an oil and gas province. The North Sea basin was formed because of Jurassic crustal extension associated with the opening of the Viking Graben to the east of the Utsira High. During the Plio-Pleistocene the North Sea basin subsided and was filled in by up to 1000 m of sediments derived from the British Isles and the rising Norwegian mainland bordering the basin (Ottesen et al. 2014 and references therein).
The sediments in the top 1000 m of the North Sea basin are generally poorly lithified and consolidation depends mostly on the content of fine particles or clay. The topmost ~ 400 m in the central North Sea basin have undergone extensive reworking due to late Quaternary glaciation cycles creating varying conditions, from open landscape with rivers to lacustrine, shallow marine and ice-covered terrain, for different time spans (e.g. Sejrup et al. 2016). Geotechnical investigations to find optimum locations for offshore platform foundations have frequently found buried channels of Mid to Late Pleistocene age and extensive heterogeneity attributed to glacial reworking of the sediment (e.g. Williams and Aurora 1982). Several authors describe overconsolidated clayey sediments at 15 to 25 m depth below dense silty-fine sands in the central North Sea basin, and argue this to be caused by glacial load rather than extensive erosion of overlying sediments (Sejrup et al. 1987;Williams and Aurora 1982).
The top ~ 800 m are part of the Nordland Group (Deegan and Scull 1977;Worsley et al. 1988) and consist mainly of sediments from an open marine environment, apart from the glacial deposits of the Pleistocene age. The North Sea experienced at least seven glaciation events during the Pleistocene age, with varying ice-sheet extents and geometries (Stewart and Lonergan 2011).
In 2011, a 3.5 km long seafloor fracture on the Utsira High was discovered (Fig. 1b, Pedersen et al. 2013). This large-scale, unexpected feature was found in an ostensibly well-mapped, highly industrialized offshore area. The discovery of the Hugin Fracture (Fig. 1b), named after the Autonomous Underwater Vehicle (AUV) Hugin that facilitated the discovery, represents a distinct structure with nearly no available geotechnical measurements (Pedersen et al. 2019). Seismic interpretation shows that the Hugin Fracture is located above an extensive (glaci-)fluvial channel network in the uppermost 400-500 m of Pleistocene test the hypothesis employing deterministic, transient poro-elastoplastic simulation for a 2D model based on the Fig. 2 a Interpreted seismic section used to build the simulation model. The red horizontal line at 0 TWT indicates the location of the sub-bottom mosaic in Fig. 1b). b Interpreted Quaternary sequence on a 2D seismic line from drilling site survey report for well 16/4-7 (modified from Fig. 4.1 in the report, Fugro Survey 2012). The units correspond to the layers in the geomechanical model and the model thicknesses are taken from the designated drilling site (broken red line). For a description of the seismic units see Table 1 interpreted seismic section in Fig. 2a and seismic from a drilling site survey report for the nearby well 3).

Method and mathematical model
Our hypothesis assumes differential compaction due to glacial load. To test this, we need to build a mathematical model that can describe such behaviour. First, we included poroelastic behaviour in the elastic material domain. Poroelasticity couples the porous fabric of the sediment to Dacian fluid flow through the pores, which will have an impact on how the sediment can be compacted. In order to introduce plastic sediment compaction into our model we included elastoplastic behaviour and a hardening rule (Fig. 4). Combining porous and plastic behaviour provides us the needed poro-elastoplastic model to simulate elastic compaction and over-compaction in seafloor sediments reported from the study area. It also allows simulating increasingly brittle behaviour of clayey sediments that obtain higher Young's modulus under stress.
There are many choices to consider when designing geomechanical simulations, one of which concerns whether to consider the gravity force. For our model, we assume sediment parameters for normally consolidated sediments, i.e. gravity and sediment accumulation. By employing such values for e.g. porosity, permeability and density, we implicitly include the constant gravity force. We focus on the dynamic change of the ice load rather than the total load (weight of ice and gravity).
The simulation presented in this study has been designed and computed with COMSOL Multiphysics. Multiphysics is a software platform of advanced numerical methods designed to investigate a wide variety of physical problems and questions (COMSOL Inc. 2016;Li et al. 2009). We employed this platform in order to investigate the response of a mechanically layered, heterogeneous geologic model subjected to vertical load.
The mathematic models presented in appendix A correspond in essence to the work flow in Multiphysics (COM-SOL Inc. 2016). The governing principles incorporated in the simulation tool are based on the work of Terzaghi (1943), Biot (1941) and Tang et al. (2015) among others. An interface couples two physical models, one poroelastic and one elastoplastic, to jointly compute a poro-elastoplastic simulation. The poroelastic model couples the porous soil and Dacian fluid flow in the soil skeleton. We want to simulate the geomechanical deformation sub-seafloor during a time-varying vertical load (Fig. 5) and employ the direct  Fig. 2b) and Arabic numerals indicate channels as in Fig. 2a. Boundary types and a linear increase in background cohesion at 7 kPa/m are indicated. The northern boundary of the model is defined by an infinity layer and symmetry is assumed at the southern boundary. Vertical movement is inhibited at the lower roller boundary while lateral movement is allowed. Sediment descriptions are listed in Table 1 and the different parameters for each layer are found in appendix B. Vertical exaggeration is 10 times, the same as for the seismic section in Fig. 2 1 Page 6 of 22 time-dependent MUltifrontal Massively Parallel sparse direct Solver (MUMPS) (Amestoy et al. 2001); this ensures an efficient coupling between the mathematical model and the geometric mesh.

Model description
Interpreted horizons and units are available from the final well report for 16/4-2 (Norsk Hydro AS 1990). A comparison with interpreted horizons from the seismic drilling site survey for well 16/4-7 (Fugro Survey 2012), about 10 km to the southeast, showed minor differences (m-scale) in layer thickness. Therefore, interpreted horizons for well 16/4-7 ( Fig. 2b) along with the description of the sediment content ( Table 1 in Appendix B) were taken as a basis for model construction.
The 2D geometry and layer parametrisation (Fig. 3) is based on interpreted seismic sections from site surveys at wells 16/4-2 and 16/4-7 (Norsk Hydro AS 1990 and Fugro Survey 2012, for details see table B-4). The seismic section in Fig. 2a reveals mostly horizontal layers to a depth of ~ 400 m, with some small-scale (< 2 km) variations. Local variations in the interval thickness between the two well sites may occur, but the layer sequence is continuous throughout a much larger area (Ottesen et al. 2014).
Numerous channel structures have been identified and the numbered channel structures in Fig. 2 are included in the 2D model. For simplicity, all layers are assumed to be horizontal and the channels are modelled as symmetric halfellipses, with the exception of channel 1. In nature, some Pleistocene layers show locally pronounced deviation from the horizontal model due to erosion. Due to lack of local or regional Pleistocene horizon maps and to avoid unnecessary complexity in the model, the local variations are ignored. However, channel 1 is modelled as asymmetric for two reasons: (1) an asymmetric channel bottom could comply with an alternative seismic interpretation than that presented in Fig. 2, and (2) an asymmetric channel will help to illuminate the significance of axis-symmetry and steep channel walls when compared to the other axis-symmetrical channels.
Sediment properties used in the model are specified in Tables 4 through 6 in Appendix B. Unfortunately, there are no geotechnical measurements available through the Pleistocene unit in the Norwegian North Sea. The shallow overburden on the Norwegian Continental Shelf has not yet been investigated in large detail. Research sampling resulted mostly in cores of only a few meters below seafloor, and petroleum industry has traditionally focussed only on reservoir rocks and the caprock immediately above. Numerous wells penetrate the Quaternary, and well logs, typically resistivity and gamma-ray, are available, that can at least give some indication of the sand or clay content along the well path. The interpretation of a sandy alluvial fan was based on resistivity and gamma-ray logs from well 16/4-2 (see NPD database). Lacking precise values from measurements, the Plastic deformation and hardening rule according to equation A12. Compared to ideal plastic deformation (red curve), where strain accumulates at constant stress from the yield point, the simulations in this study assume hardening behaviour, i.e. higher stress is required for further strain accumulation (blue curve).
T iso governs the gradient of the stress/strain curve (see also Eq. A17)  (Hamilton 1971); Young's modulus has been derived from bulk and shear modulus in the same source, employing the following relationship (Guéguen and Palciauskas 1994): Ikari and Kopf (2011) studied the influence of cohesion on clay-rich sediments by means of laboratory experiments. They found cohesion values in the range 9-70 kPa for clayrich unsheared, i.e. nearly uncompacted, sediments and in the range 15-133 kPa for clay-rich sheared (compacted by 0.9-2 MPa shear force) clays, while Sejrup and Aarseth (1987) report a shear strength of 200 MPa for the shallow till unit in B2001, near the Hugin Fracture. We therefore assume a value of 200 MPa for pure stiff clay in our model and assign a small, apparent cohesion to pure sand to avoid numerical instability.
The angle of internal friction is taken from a website for geotechnical information (Koliji 2013). Permeability values are taken from Bolton et al. (2000) and values for dynamic viscosity and compressibility of water are taken from Stanley and Button (in Sharqawy et al. 2010;Pierrot and Millero 2000) respectively). Ice thickness is based on values from Boulton and Dobbie (1993) who inferred ice sheet thickness from over-consolidation measurements in tills and clays in the Netherlands and the adjacent North Sea. Values for the Biot-Willis coefficient have been assigned in accordance with non-linear relationships between the coefficient and porosity (Wang 2000). Tests with different values for the Biot-Willis coefficient showed negligible impact on simulation results. This is in agreement with a discussion in Nermoen et al. (2013). Table 4 lists the parameters used in the layered background model and Tables 5 and 6 list the applied parameters for channel 1 and channels 2-9 respectively.

Mesh and boundary conditions
In our simulation the maximum free triangular element size was set to 100 m and the minimum to 1 m with a growth rate of 1.15. The curvature factor was 0.25. This automated setup provides good element coverage in the area of interest, and aims to minimize the number of elements and, hence, computing time. This is achieved by a high density of small elements in areas of expected large changes, i.e. heterogeneities and edges, while fewer, larger elements sufficiently represent the homogeneous areas. The simulations were run on a 16-core cluster computer with 128 GByte RAM. An average number of 188,487 degrees of freedom (DOF) were solved and the computing time was 39 min on average (for mesh statistics see appendix B). (1) The boundary on the left side of the 2D model ( Fig. 3) circumscribes an infinite element domain that virtually stretches the layer to a very large distance from the region of interest to reduce boundary effects. This part of the model will not be affected by the glacial load on the channelized part but allows for stress relief as it naturally occurs over the distance away from glacier fronts. The advantage of using an infinity layer over a lateral extension of the model is a lower number of mesh cells and, hence, reduced computing time. The right-side boundary is set to symmetry, virtually mirroring the model conditions, in order to provide numerical stability and reduce element numbers. The bottom is set to a horizontal roller boundary to allow only for lateral movement with zero vertical displacement. This way, the sediment in the model can compress and relax elastically, deform plastically and move laterally under loading, reflecting to a good extent the natural processes taking place below and in front of an ice sheet (e.g. Boulton and Dobbie 1993).

Ice load
Quaternary glaciations in Northern Europe and the North Sea have long been studied to establish time lines of advancing and retreating ice shields as well as total ice coverage and thickness (see e.g. Løseth et al. 2013;Sejrup et al. 2000. Recent work by Sejrup et al. (2016) shows that around the Hugin Fracture the ice moved northwards, similar to ice movement in the Norwegian Channel. Accordingly, the ice load in the simulation arrives from the south and propagates northwards over a time span of 10 kyears before it retreats again. Grollimund and Zoback (2000) employed step-load functions in their simulations after applying a pre-consolidation phase of about one million years. With a focus on largescale (300 km) lithospheric response to thick ice shields (> 200 m) this seems a reasonable solution. They also chose to ignore the poorly consolidated top 1000 m North Sea sediments. Our study concerns smaller scale deformation, namely the response of heterogeneous sediments; hence, we assume a pre-consolidated model and use a ramp load function to avoid numerical problems and account for a glacier slope (Fig. 5, and localization of the load in Fig. 3). In our simulation, we assume linear ice growth and decay at a rate of 0.8 MPa/10 kyears. A typical ice density of 0.931 kg/m 3 was used, corresponding to a growth rate of about 8 m ice per 1 kyears and a maximum ice thickness of 80 m.
The ramp load function grows upwards rather than laterally; this may be counter-intuitive to modern glacial movement in rapid ice flows. There are two reasons for this approach, the first being the idea of permafrost conditions prior to glacier growth supporting a gradual load increase over a larger area, rather than the lateral advance of a glacier (load) front. The second reason lies in the fact that strain is 1 Page 8 of 22 accumulated in areas of high shear stress, i.e. the glacier front, rather than under isotropic load such as below a thick icesheet.
Model symmetry is assumed for both the sub-bottom and the ice load, this equates to a modelled ice sheet with a total width of 27 km, twice the length as the model presented in Fig. 3.

Hardening rule and cohesion
The simulations in this study assume strain hardening behaviour of the modelled sediments, as shown by the blue line in Fig. 4. This means that the Young's modulus changes for increasing strains during the simulation. The material properties at the location of the Hugin Fracture are uncertain due to lack of core material. However, strain hardening is assumed to be the dominating stress/strain behaviour for the mostly clay-dominated sediments (Wood 1990;Sejrup et al. 1995). Strain-hardening behaviour is related to increasing stiffness/brittleness of the sediment as the stress/ strain condition passes beyond the yield point. The higher the isotropic tangent modulus (E Tiso ), the larger the stress required to further increase the strain (see equation A-17). Increased stiffness renders the material more prone to brittle failure, i.e. fracturing. During the simulation, the original geomechanical properties are changed locally, wherever the stress/strain conditions reach the yield point according to the hardening rule. New and higher values for e.g. Young's modulus are assigned to the affected mesh areas and the model behaviour changes to increasingly stress-resistant, and brittle, behaviour.
In addition to strain-hardening (Fig. 4), a depth-dependent cohesion is assumed (Fig. 3), reflecting an increase in cohesion with depth due to increasing compaction of 7 kPa/m.

Definition of mixing rule (k-value)
Tests and measurements on the actual sediments being absent, material properties had to be determined based on descriptions from drilling site survey reports and literature values. Sediment properties for pure sands and clays such as porosities, densities and Young's modulus are taken and derived from Hamilton (1971). In the model, the layers are assumed to be either sand-dominated or clay-dominated with some mixtures in-between as indicated by sediment descriptions from drilling site surveys. For an easier adaptation of the pure end member sediment properties to assumed mixtures in the distinct layers, we defined a mixing relation to compute the needed parameters: (2) value layer = value sand + value clay − value sand ⋅ k layer , where value may be any of the desired material properties, e.g. porosity, density or Young's modulus, and k is a layerspecific mixing value, ranging from 0 to 1. Equation 2 defines a linear mixing relation with the end member values being pure sand and clay, respectively. In the case of pure sand, a value of k = 0 yields the desired value for sand. For k = 1 , Eq. (2) yields the desired value for pure clay, and for other values of 0 < k < 1 it yields a value inbetween these endmembers. The sand/clay mixing parameters presented in Appendix B Tables 1, 2 and 3 are derived from the sediment descriptions in the drilling site survey report from well 16/4-7 (Fugro Survey 2012, see Table 1). The Aberdeen Ground Witch formation is described as soft clay with several sand layers and embedded boulders. Therefore, it was assigned a k-value of 0.2, which is closer to sand than to clay.
The deeper channels have clay properties similar to layer IV (Fisher Fm.), hereby neglecting the sand fill at the base of the channels that is reported in most studies on tunnel valley fills in Northern Europe (e.g. Kehew et al. 2012). The shallow channel 1 was assigned the properties of pure sand, which slightly exaggerates the fill interpretation based on gamma-ray and resistivity logs from wells 16/4-2 and 16/4-1 (NPD database). These simplifications should not change the qualitative results of the simulations and they are motivated by maintaining a manageable model complexity. A detailed list of material properties can be found in Appendix B.

Results
The simulations were stable throughout the 20 kyears load/ unload cycle with a peak in maximum stress after 10.4 kyears. The curves in Fig. 6a represent the stress, plastic strain and volumetric strain development at the Hugin Fracture, on the surface (layer I) above the edge of channel 1 (red diamond in Fig. 3). Figure 6b shows the same values for a point just below the bottom of channel 7 (inside layer V, green diamond in Fig. 3). Stress values are scaled to 1/100 of original value.
Overall volumetric strain values (green line) show compaction (negative values and a downward slope) during growth of the top load (ice sheet) in the first 10 kyears of simulation and dilation (increasing values on an upward slope) during decay of the top load in the following 15 kyears of simulation. The effective plastic strain (red line) increases in steps, and the shear stress (blue line) shows a triangle-shaped increase and decrease. The stress and plastic strain development over time is governed by the defined hardening rule (see section Hardening rule and cohesion) for our constitutive model. Each step in effective plastic strain illustrates a new state of brittleness (new Young's modulus).
Decreases in shear stress show that the stress is transformed into plastic strain over a wider area. After all the material has reached the according state of brittleness, the shear stress increases again until the new yield point is reached. Then the stress decreases again as all the material in the neighbouring area undergoes plastic deformation.
At 10.4 kyears the total glacial load has reduced for 0.4 kyears at a rate of 0.8 MPa/10 kyears, i.e. at that time the maximum load of 8 MPa has been reduced by 32 kPa. However, this slight decrease in glacial load results in a dramatic increase in volumetric strain and (likely) fracturing towards the next time step (Fig. 6a). However, the model was not defined with any kind of fracture behaviour, so that the part of the simulation after fracturing occurred (after 10.4 kyears) is disregarded in the remainder of this paper.
A similar pattern of stress, plastic strain and volumetric strain development, without a peak, is shown for a point just below channel 7 (Fig. 6b, green diamond in Fig. 3). Figure 7a shows the final shear stress distribution for the whole model domain. As expected, at the infinity boundary (N) no stresses accumulate because the load influence is absent, and stresses can dissipate. On the southern side of the model, stresses accumulate uniformly inside the distinct layers according to the respective geomechanical parameters.

Load-induced shear stress
On the southern isotropic part (right side in Fig. 7a) no differential stress is induced during loading and unloading. On the northern end of the model (left side in Fig. 7a) the influence of the ice sheet is absent in the infinity layer. In the transition area, where the influence of the glacial load starts and grows southwards to its maximum thickness, strong changes in stress accumulation associated with a prominent surface flexure are observed (Fig. 7a). The top Pliocene unit (layer VI) shows the response of a homogeneous half-space to the ramp-form-defined ice load with high shear stress below the ice sheet and no shear stress in front of the ice sheet. The transition from ice sheet to free surface is marked with a gradual decrease in shear stress in layer VI (Top Pliocene). The largest stresses are found in units II and III, two over-consolidated clay layers (Swatchway and Coal Pit Fm., Fisher Fm.), reflecting their increased ability to withstand shear stresses. The sandier layers, i.e. units I, IV and V can partly adjust to the stress (compaction due to grain rearrangement and pore deformation) resulting in lower shearstress values. In the heterogeneous channelized central part of the model, stress accumulation is uniformly distributed in each layer and channel. All channels show similar large shear stress, although channel 1 contains loose sand and the others stiff clay. In addition, focussed sub-vertical stress accumulations are evident at most channel edges. This pattern of differential stress zones in units I through V, seemingly connected to the heterogeneities represented by the channels, developed under a uniform ice load (see Fig. 7b for more detail). Very high stresses in close vicinity to the lowest stresses are present in focussed zones tangential to the channel edges, indicating differential compaction and rebound around these heterogeneities. Layer VI seems unaffected by these differential stress lenses. Subsequent figures focus on the heterogeneous central part of the model.
A closer look at the shear stress distribution in Fig. 7b reveals that the focussed stress accumulations consistently coincide with the channel edges, except for the northern edge of channel 1. Strong changes in stress amplitude over short distances (meter scale, see Fig. 7b) from around 150 kPa to around 400 kPa and less, appear inside the focus areas. Stress shadows with lower stress values are evident below most channels, especially below channels 7 and 8 (Fig. 7a). Note that the channels are displaced downward and northward from their original location as indicated by the black lines and numbers. Some surface deformation such as gentle elevations and small depressions are barely recognisable in Fig. 7a; they are easier to see in Fig. 7b. Figure 7b also shows the stress distribution around channel 1 in more detail. The stress magnitude is governed by material properties as can be seen at the transition from the sandier unit IV (yellow-orange) to the more clayey (red) layer III and at the transition from the clayey layer II (red) to the sandy (orange) layer I near the top of the model (80 m depth). Note the maximum stress values in the topmost subvertical part of the stress zone connected to the southern edge of channel 1, interpreted to represent a failure area like the Hugin Fracture. Inside channel 1, the shear-stress distribution shows stress variations nucleating at the underlying  Fig. 7a of the shallowest channel 1 with sand fill and asymmetric shape. Focussed tangential stress zones are clearly visible at all channel edges, except for the northern part of channel 1 channel 8 that continue through channel 1 (Fig. 7b). Similar subtle stress variation from channel 4 is barely recognized to continue inside channel 1 but can be identified by slightly less stress near the surface above the southern edge of the displaced channel 4. The northern edge of channel 2 is also the origin of a tangential stress pattern that seems to cross inside channel 1. However, the stress variation inside channel 1 seems to disappear towards the surface.
In addition, Fig. 7a and b show stress shadows, i.e. less stress compared to adjacent areas, below channels 3, 7, 8 and 9. Stress values presented in Fig. 6b show that the stress below channel 7 is about 300 kPa compared to the 325 kPa in the more isotropic southern part of the model (Fig. 7). These stress shadows result from the heterogeneities in the layers, i.e. channel bodies, that accumulate more stress than if they were filled with the same material (homogeneous case). This leads to the observed stress relief in the surrounding layer(s).

Plastic deformation and volumetric strain
We follow Lewis et al. (2009) and use effective plastic strain and volumetric strain to predict the opening of fractures. Von Mises effective plastic strain distribution in the model at 10.4 kyears is shown in Fig. 8 and the total volumetric deformation is shown in Fig. 9. In general, plastic and volumetric deformation is largest for unconsolidated sand layers I and IV, reflecting volumetric adjustments to stress accumulation. The stiff clay or till-filled channels show lower plastic deformation, but higher stress accumulation. Plastic deformation accumulates in focussed areas tangential to the channel edges, except for the northern edge of channel 1. As plastic deformation increases the brittleness of the sediment (see hardening rule and appendix A), fracturing is most likely to nucleate in the focussed areas along the edges of the channels. Plastic and volumetric deformation occurs in the same focussed sub-vertical areas tangential to the channel edges as seen in the stress distribution (Fig. 7). Different layers exhibit different amounts of plastic deformation, the softer layers (I and IV comprising loose sand while V comprises soft clay) experiencing more deformation than layers II and III comprising stiff clay (see Fig. 3 for layer parametrization). Accordingly, the only sand-filled channel (channel 1) displays larger plastic deformation values than the other channels that are filled with stiff clay (Fig. 8). Because all the other channels are assigned properties of stiff, and hence, already strained clay, little (additional) plastic deformation is expected.
A closer look at the plastic deformation around channel 1 (Fig. 8a) reveals large values in focussed areas tangential to channels 4, 8 and 6 inside the sandy layer IV (Ling Bank Fm.) and the soft clay layer V (Aberdeen Ground Fm.). The amount of plastic deformation inside the focussed areas changes abruptly at the boundary of the shallower layer III (Fisher Fm) comprising stiff clay. At channels 2 and 3 the plastic deformation inside layer IV (Ling Bank Fm.) is 0.002-0.0025 dropping to only around 0.0012 inside layer IV (Fisher Fm.). The same sudden drop is observed at the southern edge of channel 1 (Fig. 8a, b); strong plastic deformation of about 0.003 in layer I (Witch Ground Fm.) decreases to below 0.002 in layer II (Swatchway and Coal Pit Fm.). The described changes in plastic strain occur over layer boundaries between sediments of different stiffness, suggesting a dependence on the necessary yield stress (see appendix A for explanation on the hardening rule in the constitutive model).
The volumetric strain displays, equivalent to the effective plastic strain displays, are shown in Fig. 9. The overall volumetric strain is negative at 10.4 kyears, reflecting compaction by the ice load that started to decrease only 0.4 kyears ago. Figure 9a shows volumetric strain accumulation zones tangential to channel edges consistent with stress and plastic strain accumulations (Figs. 7 and 8). Dark blue coloured layers in Fig. 9 exhibit largest compaction due to elastic rebound after partial unloading. The sand-dominated layers show less compaction (volumetric strain around −0.7 × 10 −4 ) while channel 1 (sand-filled) shows about the same volumetric dilation as the host layers consisting of stiff clay (about −1.1 × 10 −4 ). Channels 2 through 8 experience strongly negative volumetric strain, channel 1 and the surrounding layers experience medium negative volumetric strain and the sandier layers I and IV (Witch Ground Fm. and Ling Bank Fm.) as well as the soft clay layer V (Aberdeen Ground Fm.) experience less strain. The strain pattern reflects more elastic rebound in the stiffer clay layers and channel fills, inflicting additional stress on the surrounding material. The areas of focussed stress and strain along the channel edges manifest in focussed strain areas, revealing in detail the form of volumetric strain they experienced. Note that there are small areas of positive volumetric strain (dilation) inside predominantly negative strain (compaction) areas, best observed below channels 8 and 9 (southern edge) but also present below the southern edges of channels 6 and 7. Figure 9a shows the same sudden drops in volumetric strain over layer boundaries as for the plastic strain. Note that the change appears also at the boundary between layers IV and V (see below channels 2 and 3 in Fig. 9a) which is less evident in the plastic deformation (Fig. 8a). Apart from the focussed strain area at channel 1, additional focussed strain areas are visible further to the south. The first one, located above the gap between channels 2 and 3, appears to be a superposition of stresses and associated strains from the same channels that are best distinguished in the shear stress display in Fig. 7. Figure 8 shows plastic strain tangents from the channel edges throughout layer III while the volumetric

Strain and deformation patterns and surface deformation
The upper boundary of layer I displays subtle topography due to differential compaction and remaining ice load (Fig. 9b). The topography becomes increasingly pronounced during unloading (ice retreat). In the following, the strain and deformation patterns are described in more detail.
In Fig. 9b, a slight surface deformation related to the strong plastic deformation above channel 1 can be seen (red rectangle). The central depression is about 0.2 cm deep (100 Fig. 9 a Volumetric strain display with hot colours for compaction and cool colours for dilation. Original model geometry indicated by black lines; Arabic numerals indicate displaced channel locations and Roman numerals indicate displaced layers. b Enlarged area around the channel 1 fracture (Hugin Fracture). Both dilatational and com-pressional volumetric strain is present in vertical elongated zones, suggesting the opening of a fracture.Vertical displacement is 100 times exaggerated. The inset shows actual sub-bottom profiler data over the Hugin Fracture, with a similar seafloor topography as the topography predicted by the model 1 Page 14 of 22 times exaggerated vertical displacement) and the plastically deformed zone is located some ten meters south of the channel edge. The surface deformation increases during continued unloading. The Hugin Fracture shows similar surface deformation (inset in Fig. 9b) with a central depression on the order of 10 cm, and the Hugin Fracture is located up to 260 m south of the edge of an alluvial fan.
A detailed look at the southern edge of channel 1 (Fig. 9b) shows strong negative volumetric strain amplitude, representing compaction, in the top layer (Witch Ground Fm.) accompanied by slightly positive amplitude, representing dilation. In the deeper layer II (Swatchway and Coal Pit Fm.) the amplitude decreases but the pattern of negative and positive values continues. The depression in surface deformation above the edge of channel 1 coincides with the strongly negative volumetric strain amplitude ( −2.6 × 10 −4 ) and the neighbouring elevations with the positive amplitude values ( 0.3 × 10 −4 ). Applied to the Hugin Fracture, the model results suggest that the elevated seafloor is associated with dilation of the underlying sediments while the depression is associated with compaction and, hence, fracture opening. For comparison, the seafloor topography over the Hugin Fracture as shown by sub-bottom profiler data is shown as an inset in Fig. 9b. The actual seafloor topography is strikingly similar to the topography predicted by our model.
A complex strain pattern is observed near the surface above channel 3 that cannot be easily linked to a single channel edge. In shear stress and effective plastic strain distribution the same complex pattern appears near the surface (see Figs. 7b,8b). It seems to represent a superposition of stresses and strains originating at channels 3, 6 and possibly 2 that is not easily recognised in either display.
The volumetric strain accumulation near the surface above channels 2 and 3 shows another pattern; it resembles an hourglass shape with negative values (volume decrease/ compaction) in the upper part, and positive values (relaxation/volume expansion) in the lower part. Some very subtle surface deformation seems to be associated with the hourglass shaped and the southernmost strain accumulation. The expected deformation should be on the order of some ten centimetres. We do not expect fracturing associated with the hourglass shaped strain pattern, but there were no subbottom profiler data available from this area to support this. The sub-bottom profiler data (inset in Fig. 9b) acquired for studying the Hugin Fracture covers only a very local area around the seafloor fracture (Fig. 1b) and the analysis was focussed on fracture indications. More widespread highresolution backscatter data (HISAS) are available and might reveal additional subtle seafloor topography.

Discussion
The motivation for the presented model and simulation results was to test and, possibly, falsify the hypothesis of fracture formation for the Hugin Fracture by differential compaction due to glacial loading/unloading. However, the presented simulation results are robust and stable in a ± 15% range around the end member and layer values presented in the tables in appendix B. As such, the model results support the hypothesis of the Hugin Fracture being a compaction fracture, but they do not prove this hypothesis to be true. Geomechanical data from the Hugin Fracture are needed to confirm and constrain model parameters. The influence of channel geometry on the stress pattern should also be investigated and could shed light on the lateral changes in the surface appearance of the Hugin Fracture.
Based on the modelling and the assessed parameter variability, it is unlikely that the Hugin Fracture is a unique feature. The simulation results show that a complex pattern of focused strain accumulation is produced by this simplified model. The simulation results are robust over a range of realistic parameter combinations and for a realistic ice load. In other words, the simulation predicts that fractures like the Hugin Fracture should be rather widespread despite having previously not been discovered. This is a main result of our study.
The high-resolution backscatter and sub-bottom profiler data acquired at the Hugin Fracture cover only a very local area at and around the Hugin Fracture. More specifically, the northern edge of the interpreted alluvial fan (solid blue line in Fig. 2a) is located beyond the data coverage of the sub-bottom profiler data and the HISAS seafloor mosaic (Fig. 1b, 10 cm × 10 cm resolution of the HISAS mosaic). In our study, the alluvial fan is interpreted as slightly asymmetrical with a steep wall at the Hugin Fracture location and gentler slope towards the northern border (see solid blue line in Fig. 2a). Therefore, less stress and hence deformation is expected to have taken place over the northern part of the alluvial fan, most likely not enough to produce any critical strain. The present paper applies a slightly different geometry, where the northern wall is even gentler than the result from seismic interpretation would indicate. This choice is motivated by the idea of testing the effect of the wall slope, and the presented results indicate that a steep channel wall has a major influence. As argued above, a symmetric sand body results in symmetric strain accumulation and, in the case of the Hugin Fracture, an additional fracture above and to the north of the northern edge of the alluvial fan.
The symmetric channels in the model produce symmetric zones of deformation and strain, tangential to the channel edges. At channel 1, however, only the steep southern wall develops a deformation zone to the surface; the gentler northern wall shows no sign of stress accumulation and, correspondingly, no strain or deformation (Figs. 7,8,9). Channel 1 appears to experience higher shear-stress levels than layers I, IV and V, although channel 1 has material properties corresponding to pure sand whereas the layers represent different kinds of sand-clay mixtures. This is a direct result of the model geometry and the properties of layers II and III that surround channel 1. Layers II and III represent stiff clay properties and, hence, they can withstand higher shear stress than the sandier layers I, IV and V that respond with porosity reduction. In the poro-elastoplastic model, the pore pressure build-up in channel 1 is influenced by the compaction of the surrounding low-porosity layers II and III. In addition, the layers are hydraulically open to the north and south, while the channel is closed. This directly affects the embedded channel structure, which is forced to undergo similar shear stress as the surrounding stiff clay layers. The pore pressure cannot dissipate through the low porous layers II and III and, hence, the resulting, effective shear stress in channel 1 is lower than the surrounding layers but higher than for comparable sandy layers (Fig. 7).
Differences in host-layer properties only impact the amplitude of the accumulated stress, not the geometry (compare channels 2 and 3 to channels below in Fig. 7). The same is true for the channel thickness, as stresses at channels 2 and 3 inside layer IV (Ling Bank Fm.) are below 400 MPa. Properties of channel fill do not change the geometry for stress accumulation zones as seen from the southern edge of channel 1 (sand-filled) which is similar to the other stiff clayfilled channels. However, the gentle wedge-shaped northern edge of channel 1 seems to prevent stress accumulation. From this we conclude that steepness of the channel walls as well as symmetry play an important role in the formation and location of stress accumulation zones.
The maximum ice load in the simulation corresponds to a little over 80 m ice, assuming a density of 0.931 kg/m 3 . Considering the absence of horizontal stress in our model this should be understood as a minimum ice thickness. In reality, and under normal stress conditions, the horizontal stress increases with depth, counteracting additional vertical load (ice), and allowing for larger load before plastic deformation occurs. In addition, our assumption of present-day compaction, i.e. stiff clays as a result of earlier glacial override or deposition, might lead to an underestimation of the actual ice load or load duration needed to compact soft clay to the stiff and brittle configuration that eventually facilitates fractures. We did also not account for glacial erosion and sedimentation that might alter the local stress field. Boulton et al. (1999) discuss various subglacial failure geometries based on their simulations of sediment deformation due to glacier-induced water flow for essentially drained large-scale models. According to their discussion, the Hugin Fracture might also represent a hydraulic fracture, in a till layer, where sandier (more permeable) sediment from the top layer has been injected (see also Figs. 8,9). However, our simulation includes plastic and elastic deformation of porous unconsolidated rocks (soil) with highly heterogeneous geometry and is restricted to undrained conditions for the chosen physics model in the simulation software. This fits best our understanding of the geological history of the study area, as does the hypothesis of differential compaction as cause for the Hugin Fracture. As there are no samples from the sediments inside the Hugin Fracture we cannot completely discount the possibility of a hydraulic fracture, although this seems rather unlikely.
The parameters employed in our model have varying degrees of uncertainty, because there were no local geomechanical measurements available for our study. However, the presented results are robust in a range of ± 15% of the layers' k-value, i.e. the model did not have to be fine-tuned to produce the results. A variation of a layer´s k-value is equivalent to a joint variation of the respective layers' geomechanical parameters. We argue that this is a better approach in our case of uncertain and interdependent geomechanical properties than the usual method of changing a single parameter at a time (Ferretti et al. 2016). The results obtained with layer parameters in a ± 15% range of the figures presented in Table 4 were all similar: In all cases, similar stress patterns with focus areas tangential to channel edges formed, and a volumetric strain peak with associated stress drop occurred at some point during the simulation run (see Fig. 6a). Figure 6a shows a clear volumetric strain peak at the surface of layer I and above the edge of channel 1 where we would expect the Hugin Fracture to be formed (measuring point indicated in Fig. 3). A similar strain peak is not seen in the area of the stress shadow below channel 7 (Fig. 6b, measuring point indicated in Fig. 3). This volumetric strain peak with associated stress drop is interpreted as the opening of a fracture at the surface. At the same time (10.4 kyears) the effective plastic strain remains constant and a drop in shear stress occurs. We want to stress the fact that the model geometry did not include any predefined fracture plains but was solely designed based on layers and channel-shaped heterogeneities. However, the stress and deformation pattern in the simulation area did vary slightly in magnitude and angle, as well as in the time when the strain peak and stress drop would occur.
We did, however, notice that small changes in the end member values for both the cohesion and the angle of internal friction did affect the simulation results and the strength of the pattern. There are relatively few published studies of these closely interrelated parameters and they include large variations (e.g. Mondol et al. 2007;Marcussen et al. 2009;Hampton 2002). Changing end member values for the Young's moduli E also has an effect, but this parameter shows less variation in the literature and seems to be better constrained. The Poisson ratio shows only minor effects for realistic values close to 0.5 for the presented model of fully saturated sediments. The remaining parameters showed negligible effects for small changes in end member values.
The simulation results show surface deformation similar to and of the same magnitude as the observed centimetrescale seafloor elevations at the Hugin Fracture location (inset in Fig. 9b). High-resolution sonar data from an AUV show a similar pattern of micro-bathymetry of the same magnitude at the Hugin Fracture as predicted by the model (see Fig. 1b and inset in Fig. 9b). In addition, the simulation results predict gentle surface deformation above channel 8, which could be confirmed by high-resolution backscatter sonar data. Unfortunately, the acquired backscatter data do not cover the area above channel 8. Nevertheless, the good agreement between the modelled elevation caused by differential compaction and the actual elevation in the microbathymetry strengthens the model results.
We conclude that the Hugin Fracture represents the first example of a soft-sediment seafloor fracture on a glacially reworked sedimentary basin. Only the extremely high resolution of the AUV-mounted sonar revealed the Hugin Fracture at the seafloor (Fig. 1b). The fracture signal is confirmed on sub-bottom profiler sections (inset in Fig. 9b) but could have been missed due to different survey design and larger line spacing. We expect that similar fractures will be found in the future, when ultrahigh-resolution seafloor mapping becomes more widespread.

Conclusion
The poro-elastoplastic simulation supports the formation hypothesis for the Hugin Fracture as a compaction fracture and suggests that thin ice sheets may induce differential compaction to a depth of several hundred meters. However, this hypothesis was meant to investigate the possibility that the Hugin Fracture was the result of glacial load and not to identify sediment properties such as chemical compaction or clay mineralogy. The results of this simulation show, however, that variation in material properties have limited effect on the simulation outcome. As long as (glacial) erosion carves paths into older sediments, and these paths subsequently fill with younger sediment, differential compaction can take place.
Shallow horizontally layered, isotropic sediments subjected to ice sheet load develop isotropic stress and strain distributions according to the respective layer properties. Heterogeneities such as channels/tunnel valleys introduce disturbances in the developing stress field. Channelized Pleistocene sediments, which are common in the Central North Sea, experience differential compaction when subjected to ice load from glaciers according to the presented poro-elastoplastic simulations. The induced shear stress shows focussed zones of stress accumulation and equivalent strain build-up tangential to channel edges. Symmetric zones of high strain develop at the steep walls of symmetric channels irrespective of sediment fill properties, if they differ from the host layer properties. Heterogeneities to a depth of about 200 m are shown to develop strong deformation accumulations along tangential channel edge areas, suggesting that even thin ice sheets (< 100 m), given favourable geometry, may induce fractures in surprisingly large depths.
The vast bulk of stress and strain accumulates during glacier retreat and beginning sediment rebound, resulting in fractures along some of the focus zones, e.g. above the shallowest sand-filled channel. The simulation was able to reproduce a strain pattern interpreted as resembling a fracture opening connected to the edge of this channel with associated surface deformation similar to the Hugin Fracture. Hence, the formation hypothesis for the Hugin Fracture, being a compaction fracture, is supported by the simulation results.
Based on the application of standard literature values for our model, the modelling results indicate that similar fractures should be expected at other places with favourable geometry. The influence of the geometry on the simulation results, as well as the geomechanical parameters, should be further investigated.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix A
Fluid flow in a porous medium due to the hydraulic potential field is described by Darcy's law, with neglected gravity commonly written as: where κ is the permeability, μ the fluid dynamic viscosity, p the fluid pressure, q the Darcy velocity and ∇ the nabla operator. Formulated with the poroelastic interface equation (A1) becomes: Q m is the mass source term, φ the porosity and ρ the fluid density. In order to combine the solid mechanics with the Darcy flow, the storage coefficient S t needs to be introduced. The storage coefficient can be defined as: χ is the fluid compressibility. The poroelastic equation (A2) can be substituted by the equations in (A3) to: The linear solid mechanics theory describes that the deformations are proportional to stress and reversible. This assumption is well known as generalised Hooke's law: In Einstein notation and with σ = stress tensor, K the bulk modulus, τ the shear stress and G the shear modulus. A poroelasticity interface was used to couple the linear solid mechanics with Darcy's law and to account for poroelastic deformation. Based on Tang et al. (2015), the interface can be mathematically summarized as: and with the Biot-Willis coefficient α B , reference pressure p ref , the identity matrix I, p f and ρ f representing the fluid pressure and ε vol the fluid density strain. The Biot-Willis coefficient describes the interaction between confining stress and pore pressure.

Porous matrix deformation
For the poroelastic stress tensor we consider an isotropic porous material under plain strain conditions. Equation A8 describes the norm for a 2D poroelastic (Wang 2000): with E = Young's modulus, = Poisson's ratio of the porous material for the drained case, B = Biot-Willis coefficient and p = fluid pressure. The term B ⋅ p is often described as fluid-structure interaction. Since we added poroelasticity to our model, equation A8 contributes as a poroelasticity node in the simulation set up.

Elastoplastic model
The elastoplastic model is described through linear isotropic Young's modulus E, Poisson ratio ν and density ρ. Plasticity is assumed with small strain approximation yielding following equation for the elastoplastic model (e.g. Tang et al. 2015): with incremental strain dε, effective plastic incremental strain dε ep , μ = E/2(1 + ν) and λ = E ν/((1 + ν)(1 − 2 ν)). The relationship between the strain tensor ε and displacement u is: The soft soil in our simulation has a Poisson ratio of almost 0.5, causing instability in volumetric change calculations. This can be illustrated by the bulk modulus that measures volumetric change and tends to infinity as the Poisson ratio approaches 0.5: To avoid an ill-posed numerical problem, we employed a mixed formulation where we added a dependent pressure variable to the deviatoric stress tensor S and now need to solve a nearly incompressible problem: with C as fourth order constitutive Cauchy stress tensor; the " ∶"-symbol indicates a contraction over two indices ( = C ∶ = Cijkl ⋅ kl). Poroelastic fluid flow under static conditions and without body force can be described as: with σ e the effective stress tensor and p the pore pressure. Displacement in the poro-elastoplastic model can be described by inserting equations A3 and A4 in equation A1, yielding:

Yield criterion and hardening rule
Linear elastic materials deform under load and relax to its original form when the load is relieved. If a certain stress level (yield stress) is exceeded in an elastoplastic material, like soft sediment, irreversible plastic deformation takes place and increases during unloading.
In our simulation we consider the original sediment ductile rather than brittle because the porous sediment skeleton will deform under tensile stress without fracturing (until it reaches the maximum yield stress). Using this assumption, the von Mises yield function is a useful approximation. Von Mises yield criterion is defined as: where σ i is the main stress direction with i = 1,2,3. The resulting yield stress function F is then: with the yield stress σ y . In our simulation we also include hardening during plastic deformation and we assume linear isotropic hardening with the isotropic tangent modulus E Tiso .
The yield stress σ y then is the summation of the initial yield stress σ 0 and the hardening stress E hard ε ep which can be expressed as: