Models of overpressure build-up in shallow sediments by glacial deposition and glacial loading with respect to chimney formation

Chimneys and pipe structures have been observed in the caprock above the Utsira Aquifer in the North Sea. The caprock is of Pleistocene age and the chimneys appear to have been formed by natural hydraulic fracturing towards the end of the last glaciation. We study six different models for the pressure build-up in the Utsira Aquifer with respect to chimney formation. The first two models produce overpressure by a rapid deposition of glacial sediments. Using these two models, we show that the caprock permeability must be as low as 100 nD for sufficiently strong overpressure to develop. This value seems to be one order of magnitude lower than the measured permeabilities of the caprock. The four remaining models produce overpressure by a glacial loading of the caprock and the aquifer. This study shows that a 1-D model of a caprock with soil properties cannot produce conditions for chimney formation unless the least horizontal compressive stress is much less than the overburden. Furthermore, a 1-D poroelastic model of glacial loading of an aquifer and a caprock cannot produce conditions for chimney formation based on available geomechanical data. However, we demonstrate that a 2-D poroelastic model can produce conditions for chimney formation with glacial loads that partially cover the surface.


Introduction
The Utsira Formation is a sand aquifer in the Southern Viking Graben in the North Sea (see Fig. 1); it has as a caprock the Nordland Shale, which is characterised by seismic chimneys (Kartens and Berndt 2015). Seismic chimneys and pipes are vertical seismic anomalies interpreted as leakage structures, which allow fluids to escape from otherwise sealed aquifers (Berndt 2005;Kartens and Berndt 2015). Chimney structures are common features in sedimentary basins worldwide (Berndt 2005;Cartwright et al. 2007;Løseth et al. 2009). Several features of the chimneys are poorly understood. One important issue is their formation process, although chimneys and pipes can be explained by hydraulic fracturing of an impermeable caprock (Berndt 2005;Cartwright et al. 2007;Løseth et al. 2009;Kartens and Berndt 2015;Zhang et al. 2018). Assuming that natural hydraulic fracturing is the process that generated the chimneys rooted in the Utsira Aquifer, then the aquifer must once have had a high overpressure. Overpressure is defined as the difference between fluid pressure and hydrostatic pressure. The Nordland Shale caprock is of Quaternary age and the young formation age of the chimneys (Pleistocene to Holocene) hints at a connection between their formation and loading by ice (Kartens and Berndt 2015).
There have been observations of blowout structures on the prairies of North and South Dakota, USA (Christiansen et al. 1982;Zhang et al. 2018). These structures, despite their different geological setting, could be the result of the same hydraulic fracturing process that occurred in the caprock of the Utsira Formation. For the Dakotas, Christiansen et al. (1982) reported a hydrodynamic blowout structure with a diameter of nearly 300 m. It appears as a circular crater lake with a greater depth than other lakes in the same area. Observations indicate that the crater formed as a hydrodynamic blowout from an overpressured aquifer 400 m beneath the surface. The blowout has been dated to between 12,500 years and 12,000 years BP, coinciding with the end of the last glaciation in North America. Several similar blowout structures have been mapped in the same area (Zhang et al. 2018); therefore, the blowout structures 1 3 appear to result from overpressure build-up produced by glacial loading (Neuzil 2012;Zhang et al. 2018).
A model for Pleistocene pressure build-up in the Utsira Formation is needed to explain the chimney formation in terms of a hydraulic fracturing process. Chimneys are assumed to form by hydraulic fracturing when the reservoir pressure exceeds the least compressive horizontal stress (Løseth et al. 2009;Cartwright and Santamarina 2015). Furthermore, it is assumed that the vertical stress is the largest principal stress. Therefore, the aquifer fluid pressure must reach a substantial fraction of the overburden to produce chimneys. It should be mentioned that predicting pore fluid pressure and stress is a challenging task (Almalikee and Al-Najim 2018; Abdideh et al. 2020).
We studied six models for pressure build-up in the Nordland shaley caprock and the Utsira sand/sandstone aquifer. These models are labelled A to F. The first two models (A and B) produce overpressure by rapid deposition of Quaternary glacial sediments. The rapid deposition of low-permeability sediments is a process known for generating overpressure (Audet and Fowler 1992;Wangen 2010). These two models have different sediment rheologies. The next four models (C to F) produce overpressure build-up in the Quaternary sediments by glacial loading, another process responsible for overpressure (Neuzil 2012;Person et al. 2012;Zhang et al. 2018). One of these models (C) has soil properties, which is the same sediment rheology as in model B. The last three models are poroelastic. Poroelasticity allows for the computation of both the fluid pressure increase and the lateral stress increase. Models C and D are 1-D with uniform glacial loads, while models E and F are 2-D and have glacial loads that partially cover the surface.
The effects of glacial loading and unloading on hydrodynamics in sedimentary basins pose several challenges, as shown by numerical studies of Pleistocene glaciations in North America. These challenges include the thickness and timing of the ice sheets, pressure boundary conditions beneath the ice, groundwater recharge in aquifers, the stress and pore fluid pressure produced by compression and tension from the ice loads on a lithospheric scale, blowout structures generated by overpressure in shallow aquifers, brine mixing in aquifers by fluid flow, and present-day observations of underpressure (Neuzil 2012).
The aquifer overpressure may have leaked out through the chimney structures (Løseth et al. 2009(Løseth et al. , 2011, and these structures could be leakage pathways even today. The degree to which chimneys could be leaking at present is important to know for reservoir units storing CO 2 , such as Utsira. Therefore, it is necessary to understand the hydrodynamic processes that led to the formation of the chimneys. Although the Utsira Formation is only one specific aquifer in the North Sea, it could be relevant for other places in the North Sea area and in the Norwegian and Barents Seas. These places share a common history of shaping by glacial processes during the Quaternary.
This article is organised as follows. The six models for overpressure build-up are presented in turn, with the modelling approach and results presented for each model. The presentation of the six models is followed by a discussion of their ability to produce chimneys.

Method
Rapid deposition of sediments is a process that generates overpressure in sedimentary basins. The equation for the generation of overpressure is based on mass conservation of pore fluid and solid matrix in a compacting porous medium. Appendix A provides a derivation of the pressure equation. The Gibson model (Gibson 1958) is an exact solution of the equation for overpressure when the deposition is at a constant rate. The solution has been adapted for an arbitrary amount of sediment compaction (Wangen 1993(Wangen , 2010, and it is a simple means to predict overpressure from a deposition process. The Gibson model has  Fig. 1 The extent of the Utsira Aquifer in the North Sea. The axes have UTM-coordinates in units of 100 km the void ratio e = ∕(1 − ) , where is the porosity, as a linear function of the effective pressure p s = p l − p f , where p l is the lithostatic pressure (weight of the overburden) and p f is the fluid pressure (see Appendix A). The void ratio function of the Gibson model is in terms of the effective pressure p s , where the surface void ratio is e 0 (Gibson 1958). The permeability is assumed to be linear in the void ratio where k 0 is the surface permeability (Wangen 1993). Gibson's analytical expression for the overpressure, as a function of time and depth, is given by Eq. (26) in Appendix A. The thickness of the caprock of the Utsira Formation is in the range of 500 m to 1500 m, and the current water depth is from 90 to 120 m (Eidvin and Rundberg 2001;Bøe et al. 2002;Eidvin and Rundberg 2007;Halland et al. 2012). The caprock consists of clay-type glacial deposits and has a Pleistocene age (Schepper and Mangerud 2017). These low-permeability sediments were deposited during a time interval of less than 3 Myr. The models used to simulate pressure build-up during the deposition of the caprock sediments do not include the deposition of the aquifer. The fluid pressure in the aquifer is assumed to be equal to that at the base of the caprock.
Overpressure build-up is studied using the exact solution in five cases with different surface permeabilities: k 0 = 4.63 ⋅ 10 −21 m 2 , k 0 = 4.63 ⋅ 10 −20 m 2 , k 0 = 4.63 ⋅ 10 −19 m 2 , k 0 = 4.63 ⋅ 10 −18 m 2 a n d k 0 = 4.63 ⋅ 10 −17 m 2 , which are labelled 1, 2, 3, 4 and 5, respectively. The remaining parameters in the study are listed in Table 1. These five permeability coefficients were selected because they give the gravity numbers N g = 0.01 , 0.1, 1, 10 and 100, respectively. The gravity number is defined as where f is the pore-fluid density, is the deposition rate, g is the gravitational acceleration and is the pore-fluid viscosity. Sediment deposition takes place at a constant rate = max ∕t max , where max = 700 m is the present-day total amount of net (porosity-free) sediment in the caprock and t max = 3 Ma is the duration. The gravity number is a useful parameter because it distinguishes between three regimes of overpressure build-up. First, there is a near-hydrostatic regime with small overpressure, characterised by gravity numbers N g ≫ 1 . Next, there is an overpressure regime, where the fluid (1) e(p s ) = e 0 − C g p s , pressure is close to the lithostatic pressure, which is characterised by gravity numbers N g ≪ 1 . Finally, gravity numbers N g of the order 1 give moderate overpressure build-up.

Results
The fluid pressure at the base of the caprock must be a substantial fraction of the overburden to produce chimneys. Figure 2 shows the overpressure, porosity and permeability for the five different surface permeability cases. A permeability as low as k 0 = 4.63 ⋅ 10 −20 m 2 is needed for the deposition to produce high overpressure. This permeability gives a gravity number of N g = 0.1 . There are few available measurements of the permeability of the caprock of the Utsira Field. However, one of these indicates that the permeability is on the order of 1 ⋅ 10 −19 m 2 (Harrington et al. 2009). This measurement was made towards the base of the caprock, and it is uncertain to what degree it is representative of other parts of the caprock. A permeability of 1 ⋅ 10 −19 m 2 is sufficient to produce moderate overpressure, N g ∼ 1 , but it might be one order of magnitude too large for strong overpressure to develop. Another set of measurements, made by Springer and Lindgren (2006) on samples from the Sleipner Field, indicates that the caprock permeability is as high as 1 ⋅ 10 −18 m 2 . These permeabilities are too high for the deposition process to generate the necessary overpressure to trigger chimneys.

Method
Overpressure build-up by the rapid deposition of glacial sediments with soil-like properties was studied. Laboratory experiments involving clay compaction show that the void ratio in soils decreases nearly linearly with the logarithm of the effective pressure (Atkinson and Bransby 1978;Whitlow 2001;Harrington et al. 2009) where e v is the void ratio at the effective pressure p v , where the consolidation index C L controls the rate at which the void ratio decreases with increasing effective stress. The permeability function for the clay model is which is an extension of the permeability function (2) with a nonlinear term, which is similar to the Kozeny-Carman equation when the exponent is n 0 = 3 . The pressure equation was solved numerically as explained in Gutierrez and Wangen (2005). The clay model was applied to the same five cases of surface permeability as the Gibson model: k 0 = 4.63 ⋅ 10 −17 m 2 , k 0 = 4.63 ⋅ 10 −18 m 2 , k 0 = 4.63 ⋅ 10 −19 m 2 , k 0 = 4.63 ⋅ 10 −20 m 2 and k 0 = 4.63 ⋅ 10 −21 m 2 . Table 2 shows the other parameters used in the simulations. Therefore, these simulations have the same gravity numbers as for the Gibson solution in model A. Figure 3 illustrates the resulting overpressure, porosity and permeability. The compaction curve (5) in Fig. 3b shows that the porosity under nearly hydrostatic conditions decreases from = 0.5 at the surface to ≈ 0.15 at the base of the caprock when the compaction index is taken to be C L = 0.1 . This is a reasonable porosity development for shallow hydrostatic clays/shales (Lee et al. 2020). It should be noted that Harrington et al. (2009) reported a consolidation index for the Utsira caprock in the range of 0.004-0.013, with e v ≈ 0.6 for p v ≈ 4 MPa , but these data provide almost no compaction

Results
at the base of a 1000 m-thick caprock under hydrostatic conditions. The overpressure curve (4) in Fig. 3 shows weak overpressure build-up and corresponds to curve (3) in Fig. 2. The reason for this is that the permeability drops by an order of magnitude from the surface to the base of the caprock. The permeability coefficient must be as low as k 0 = 1 ⋅ 10 −19 m 2 for the overpressure to be nearly as high as the lithostatic pressure. The permeability at the base of the caprock is then 1 ⋅ 10 −20 m 2 . Therefore, the restriction on the permeability is the same as for the Gibson model, where the permeability is one order of magnitude too low when compared to the measurements of Harrington et al. (2009).

Method
The glacial loading of the surface compacts the sediments and creates overpressure. It is not evident what the boundary conditions should be at the seafloor in the case of glaciation. Three different scenarios are shown in Figs. 4. Figure 4a shows an ice sheet floating above the seafloor; in this case, the ice exerts no pressure on it. From the moment the ice touches the seafloor, it becomes a load, as shown in Fig. 4b. This load may be negligible, however, if most of the ice is supported by the water. Such a glacier is termed 'wet-based' if the temperature at its base is close to the pressure melting point (Hambrey and Glasser 2012). The situation is different for a dry-based glacier that is frozen to the seafloor. In this case, there is permafrost beneath the glacier, and there is no water providing buoyancy to the ice sheet, as indicated in Fig. 4c. The mapping of areas that were covered by cold-based ice during the Last Glacial Maximum is an unresolved problem, and the thickness of the Fennoscandian/Scandinavian ice sheet is poorly constrained (Olsen et al. 2013). In the following, it is assumed that the glacier is frozen to the seafloor. Therefore, there is no water above the seafloor adding to the pore-fluid pressure in the sediments, and the entire weight of the glacier contributes to the effective stress on the surface of the basin. For the sake of simplicity, the overpressure was assumed to be zero at the frozen seafloor. Only a fraction of the porosity lost during loading and compaction is recovered during unloading. The porosity function of effective stress that accounts for porosity rebound is a generalisation of the porosity function (4), and it is where e min is the minimum void ratio at the maximum effective pressure p s,max and C U is the recompression index (Fowler and Yang 2002). The recompression index is less than the compression index, C U < C L , since decompaction is less than compaction. Reloading is assumed to follow the same e − p s -path as unloading until p s once more reaches p s,max and normal compaction continues (Fowler and Yang (3)  (6) is shown in Fig. 5a for the parameters from Table 2. The aquifer pressure is assumed to be the same as the fluid pressure at the base of the caprock.

Results
The overpressure development was studied for the glacial loading history shown in Fig. 5b. The load increased from zero to its maximum value of 10 MPa over 45 kyr, followed by 5 kyr of unloading back to zero, and finally a 10 kyr pause, up to the present day. The load as a function of time is a simple representation of a complex loading history, involving several hotter and colder intervals.
The maximum load of 10 MPa corresponds to an ice sheet nearly 1200 m thick. The resulting overpressure history is shown in Figure 6 for a 1000-m-thick caprock. The three plots in Figure 6 differ in the caprock permeability coefficients k 0 , which are k 0 = 1 ⋅ 10 −19 m 2 , 1 ⋅ 10 −18 m 2 and 1 ⋅ 10 −17 m 2 in (a), (b) and (c), respectively. The overpressure in the caprock follows the load closely as it increases linearly with time. This is as expected, because the Darcy flow in the caprock is almost zero. A zero-flow term (undrained conditions) implies that the pressure Eq. (22) can be approximated as where L(t) is the load as a function of time. The compressibility of the rock dominates the fluid compressibility, |d ∕dp s | ≫ c f , which implies that p∕ t ≈ dL(t)∕dt , and p(t) ≈ L(t) , assuming that the overpressure was initially 0. Therefore, the maximum increase in overpressure is the size of the load. This approximation (7) applies for both loading and unloading when the fluid flow is negligible. However, in the last case, with a permeability of k 0 = 1 ⋅ 10 −17 m 2 , fluid flow is noticeable, and underpressure develops at the present time as seen from Fig. 6c. (3) (1) (a) (b) 0.01 -Permeability exponent ( n 0 ) 3 -Surface permeability (highest) ( k 0 ) 1.65 ⋅ 10 −17 m 2 Surface permeability (lowest) ( k 0 ) 1.65 ⋅ 10 −21 m 2

Method
This section presents a poroelastic model of pressure buildup by glacial loading. It is assumed, as in the preceding section, that the glacier is frozen to the seafloor and that the load is the entire weight of the ice. The poroelastic pressure increase from a load L(t), uniformly distributed on the surface, becomes proportional to the load as assuming undrained conditions, where the factor f p is Skempton's coefficent (Wang 2000), and where is the Biot coefficient, S is the storage coefficient and M is the constrained modulus. The derivation of the pressure increase (8) is given in Appendix B. It is assumed that the strain is zero in the horizontal plane and that the pressure is initially hydrostatic. The Skempton's coefficent f p can be expressed as a function of the Biot coefficient and other poroelastic parameters, using where K is the drained bulk modulus, K f is the fluid modulus and is the Poisson ratio. The increase in lateral compressive stress is also proportional to the load where the coefficient of proportionality is denoted by f h . The product S can be written as the following function of the Biot coefficient See Appendix B for a derivation of the expressions (10) and (11). In contrast to the model C, where the sediments have soil properties, the poroelastic porosity fully recovers during unloading.

Results for the aquifer
The factor f p is plotted as a function of the Biot coefficient in Fig. 7, which shows that f p is close to 1 over a wide range of in the case of K∕K f ≪ 1 . The overpressure is then close to the surface pressure of the load. In this case, the matrix is much more compressible than the fluid; therefore, the fluid carries most of the load. In the other regime, where K∕K f ≫ 1 , the factor f p is close to zero. In this regime, the matrix is much less compressible than the fluid, and it is the matrix that carries most of the load.
The condition for pressure build-up to be comparable to the load is a bulk modulus K that is much less than K f = 2.5 GPa , a typical value for water. No direct measurements of the bulk modulus for the Utsira sand seem to have been made, probably because the sand has a weak rock frame. The Utsira sand/sandstone is fine-grained, weakly consolidated, highly porous (30-40%) and permeable ( 1 ⋅ 10 −12 m 2 to 3 ⋅ 10 −12 m 2 ) (Arts et al. 2008;Subagjo 2018). Elenius et al. (2018) used the value K = 0.29 GPa in their geomechanical simulations of CO 2 -storage in the Utsira Formation. The fact that the Utsira Sandstone is (5) (6) (7) (8) weakly consolidated, with a high porosity, supports the idea that its drained bulk modulus could be much less than the fluid modulus. Therefore, it was concluded that the Utsira sand/sandstone would respond with a poroelastic pressure increase close to that of the load. Figure 8 shows the factor f h as a function of the Biot coefficient for different ratios of the drained bulk modulus over the fluid modulus ( K∕K f ). The factor f h is approximately 1 for the case of ≈ 1 when the drained bulk modulus is much less than the fluid modulus ( K∕K f ≪ 1 ). These two conditions could be fulfilled for the Utsira caprock.

Results for the caprock
Geomechanical measurements are available for the Nordland Shale. Zweigel and Heill (2003) made two measurements of Young's modulus for samples at the base of the Nordland Shale in the Sleipner Field, finding E = 0.2 GPa and E = 0.3 GPa , with corresponding Poisson ratios of = 0.25 and = 0.18 . Since E ≪ K f for these measurements, and since the bulk modulus is always less than Young's modulus ( K < E ), this implies that K ≪ K f in the caprock. Harrington et al. (2009) reported similar results for samples taken from the base of the Nordland Shale, with values in the range of K = 0.05 GPa to K = 0.5 GPa.
No measurements of the Biot coefficients seem to exist for the Nordland Shale. Zweigel and Heill (2003) assumed that = 1 for the caprock. There have been reports of Biot coefficients for clay/shale with values of less than 0.5 (Luo et al. 2015); however, these values are for more deeply buried rocks than the Nordland Shale, and with much less porosity. It is not unreasonable that the Biot coefficient could tend towards 1 for the Utsira caprock. Therefore, the overpressure increase in the aquifer, and the lateral stress increase at the base of the caprock, would both be nearly equal to the pressure from the glacial load.

Method
The fluid pressure in an aquifer must exceed the least lateral compressive stress at the base of the caprock to create a chimney. This condition is difficult to achieve if the lateral stress increases by the same amount as the aquifer pressure. It is now shown that situations exist in which the aquifer pressure could increase by nearly the size of the load, while the lateral stress in the caprock remains almost unchanged.
As already discussed, the compressive stress in a low-permeability caprock depends on the spatial distribution of the surface load. If the load is absent over an area directly above a point in the caprock, then it is possible that the increase in the lateral compressive stress is almost absent at this point as well. On the other hand, the fluid pressure in a permeable aquifer depends less on the details of the spatial distribution than on the spatial average of the load.
Finding the depth of a stress-free zone beneath an unloaded part of the surface is not straightforward. To get an idea, we studied a simple situation with two parallel loads separated by an unloaded band, as shown in Fig. 9. The situation illustrated in Fig. 9 was simulated using a 3-D K/Kf=0.01 K/Kf=0.1  . Figure 10 shows a vertical cross-section through the 3D model. The base of the model is 58 km x 58 km. The surface load covers most of the surface, except for the 4.7-km-wide band placed symmetrically around the y-axis. The subsurface consists of a 700-m-thick shale layer, which rests above a 50-m-thick aquifer of sand, as shown in Fig. 10a. This situation could be similar to that of the Sleipner Field in the North Sea, assuming it was partially covered by two parallel ice sheets with an open space between them.

Results
The two parallel loads in Fig. 10a increase linearly from 0 Pa to 10 MPa over a time span of 1000 years. The loading gives a linear overpressure increase with time, as shown in Fig. 11. Equation (8) for the fluid pressure gives an excellent match against the numerically computed fluid pressure using the spatial average of the surface load. The overpressure is almost the same everywhere in the sealed aquifer because lateral overpressure gradients from the loading dissipate faster than the rate of loading. The dissipation of pore pressure gradients is a rapid process when permeability is high. The poroelastic parameters for the numerical model are shown in Table 3, and the properties of the pore fluid can be found in Table 4. Figure 10b shows the fluid pressure in a vertical crosssection through the model. It is seen that the low-permeability caprock has an increased fluid pressure beneath the load, but almost no increase in the fluid pressure where the surface load is absent. This is because the permeability in the caprock is sufficiently low for the fluid to be nearly immobile on the time scale of loading. Figure 10c illustrates the stress increase in the x-direction ( xx ) from the surface load. The caprock has almost no stress increase beneath the band of zero surface load.
The lateral stress increase is almost absent all the way to the base of the caprock. A zone exists at the base of the caprock where there is a transition from nearly zero lateral stress increase to a nonzero stress increase in the aquifer. Despite this thin zone at the base of the caprock, the preferred places for the generation of leakage structures, such as chimneys, may be where the surface loads are small or absent.

Method
It is not a straightforward matter to obtain analytical expressions for the stress increase beneath non-uniform surface loads. One exception is the stress beneath two parallel loads, as shown in Fig. 9. In this case, the stress increase turns out to be proportional to the bulk strain, as shown in Appendix C. It can then be assumed that yy = 0 . The undrained pressure increase is always proportional to the bulk strain, and the bulk strain can serve as a proxy for lateral stress. Figures 10b and 10c show that the places with a near-zero increase in lateral stress are nearly identical to the places with no pressure increase. It is possible to compute the bulk strain produced by the two parallel ice sheets in Fig. 9 by using the Boussinesq solution, as shown in Appendix C. The bulk strain becomes where and where L is the surface pressure from the two loads (see Appendix C). The parameter G u is the shear modulus obtained from the undrained bulk modulus. The two ice sheets are separated by a distance of 2a, and are placed symmetrically around the y-axis. The function (13) for d(x, z) explains how the bulk strain depends on position (x) and depth (z). Expression (12) shows that the bulk strain has the maximum value kk = (1 − 2 ) L∕G u for z → ∞ , and that it approaches 0 at shallow depths beneath the band ( |x| < a and z → 0). Figure 12 shows a plot of d(x, z) for different depths measured by the length scale a. It can be seen that, for depths much less than a ( z ≪ a ), there is a weak pressure build-up in the caprock beneath the band where the load is absent. In

Results
the case of an aquifer such as Utsira, with depths down to 700 m, the load would have to be absent over a band much wider than 1400 m. The distance separating the two loads in Fig. 10 is 2a ≈ 6 km , which is sufficient to fulfil the condition z ≪ a all the way to the top of the aquifer. It is possible to estimate the minimum load needed for the aquifer fluid pressure to exceed the lateral stress at the base of the caprock. The maximum increase in overpressure in the aquifer is the spatial average of the surface load, while it is possible for the caprock to avoid a lateral stress increase in areas where that load is absent. Therefore, the aquifer fluid pressure can be estimated as where z a is the depth to the aquifer and L is the spatial average of the surface load. It is assumed that the aquifer was initially hydrostatic. The initial least compressive stress in the caprock is a fraction f a of the overburden The condition for chimney formation is an aquifer fluid pressure greater than the least lateral compressive stress at the base of the caprock, p f > h , which implies that the load must be larger than With f a = 1 , b = 2300 kg m −3 , f = 1000 kg m −3 , z a = 700 m and g = 10 m s −2 , the spatial average of the load must be greater than the minimum L min = 9.1 MPa . A factor f a = 0.7 , which is not unrealistic, gives a minimum value of L min = 6.4 MPa for the load. The estimate L min = 9.1 MPa indicates that a glacier with an average (14) p f = f gz a +L,  Fig. 11 Pressure increase in the aquifer plotted as a function of time thickness of at least 1 km would be needed to create an overpressure sufficient to exceed the initial lateral stress at the base of the caprock, thereby producing chimneys.

Discussion
Rapid deposition of glacial sediments creates overpressure in sedimentary basins. Glacial loading produces additional overpressure. The condition for chimney formation is fluid pressure in the aquifer exceeding the least compressive horizontal stress at the base of its caprock. It appears that among the six studied models for overpressure build-up, the last two, models E and F, with glacial loads that partially cover the surface, are those which provide the best conditions for chimney generation. These are the only models that allow for overpressure build-up in the aquifer while the least compressive stress remains almost unchanged in the caprock. It was estimated that the average ice thickness must be larger than 1 km for the fluid pressure to exceed the compressive stress at the base of the caprock at 700 m.
The first two models, A and B, where overpressure is generated by rapid deposition of sediments, meet the conditions for chimney generation when the caprock permeability is as low as 1 ⋅ 10 −19 m 2 . Under these conditions, the fluid pressure approaches the overburden at the base of the caprock. The few existing permeability measurements from the Utsira caprock indicate that 1 ⋅ 10 −19 m 2 is one order too large. Even if the overpressure by rapid deposition is not high enough for chimney generation, it contributes to the overpressure build-up. An aquifer that is initially overpressured needs less pressure increase by glacial loading before reaching the critical pressure for chimney generation.
Model C generates overpressure by glacial loading of sediments with soil properties. The pressure increase is equal to the pressure from the load when the surface permeability is lower than 1 ⋅ 10 −18 m 2 . In soft soil-like sediments, the increases in overburden with the weight of the glacial load and the stress state could be near isotropic. Therefore, the least compressive stress may be close to the overburden, and the conditions for chimney formation are not met.
Model D generates overpressure in the Utsira formation by glacial loading of poroelastic sediments. It is possible with this model to select poroelastic parameters such that the fluid pressure in the aquifer increases to the same degree as the load pressure, while the least compressive stress in the caprock increases with only half of the pressure from the load. For the aquifer to have pressure build-up close to the size of the load, it needs a drained bulk modulus that is much less than the fluid modulus ( K ≪ K f ) and a Biot coefficient ( ) close 1. This is most likely the case, according to the measurements on the Utsira sand. For the lateral stress in the caprock to be nearly half the load, it needs a drained bulk modulus that is much greater than the fluid modulus ( K ≫ K f ) and a value for the Biot coefficient ( ) that is around 0.5. Geomechanical measurements indicate that the Utsira caprock is too soft and porous for this to be the case.
Gas generation from deeper formations, and its subsequent migration into the Utsira Aquifer, is one process that could have generated overpressure in the aquifer (Kartens and Berndt 2015). Gas accumulations at the top of the Utsira Formation cause overpressure. The low-permeability Nordland Shale caprock prevents further migration towards the seabed. The combined contribution from all processespressure build-up by rapid deposition, loading by the Fennoscandian ice sheet and gas accumulation-could be the reason for the chimneys in the caprock. Gas accumulation in the Utsira Aquifer was not analysed for this study, but it involves the generation and expulsion of hydrocarbons from source rocks, as well as migration through shale layers. It also involves the Skade Aquifer that lies beneath the Utsira Formation, and the way in which these two aquifers are hydrologically connected (Gasda et al. 2017).
Additional models exist that could provide the necessary conditions for chimney formation above the Utsira formation: for instance, models of heterogeneities that naturally focus the stress and the fluid flow in the caprock. It should also be mentioned that the stress state in the caprock is unknown and that visco-plastic effects could be important. It is also possible that chimney formation is related to gas hydrates in shallow sediments (Virs 2015;Waage et al. 2020).

Conclusions
The Pleistocene caprock above the Utsira sand/sandstone aquifer in the North Sea has several chimney structures. The chimneys were revealed by recent and improved seismic imaging. These structures appear to have been formed by hydraulic fracturing, triggered by high overpressure in the Utsira Formation. The conditions for chimney formation require that the fluid pressure in the aquifer exceeds the least horizontal compressive stress at the base of the caprock. Therefore, high overpressure must have developed towards the end of the Pleistocene for these chimneys to have formed by hydraulic fracturing.
Six models, labelled A to F, for overpressure build-up were tested to see whether they could produce conditions for chimney formation: Models A and B generate overpressure by sediment deposition and models C to F generate overpressure by glacial loading. The latter models assume that the glaciers are frozen to the surface.
The two models A and B have different rheologies, but both models show that the sediment permeability must be as low as 1 ⋅ 10 −20 m 2 for the fluid pressure build-up by rapid deposition to approach the overburden. This value appears to be one order of magnitude too low compared with the current permeability measurements from the caprock.
It is shown that model C, which is glacial loading of a caprock with soil properties, generates fluid pressure approaching the size of the load when the permeability is limited above by 1 ⋅ 10 −18 m 2 . Conditions for chimney formation cannot be obtained with this model, unless the least compressive horizontal stress is less than the overburden.
Model D produces overpressure by glacial loading of poroelastic sediments. It is shown that certain choices of poroelastic properties for the aquifer and for the caprock could provide conditions for chimney formation. The right choice of parameters gives a pressure increase in the aquifer that is close to the size of the load, while the least compressive horizontal stress in the caprock increases with only half of the load. The few geomechanical observations that exist for the Utsira and its caprock indicate that the aquifer is inside the parameter range for chimney formation, but that the caprock is not.
Models E and F produce overpressure by glacial loads that partially cover the sediment surface. It is shown that the overpressure in the aquifer increases with the pressure from the average weight of the glacial load, while the least compressive stress can remain nearly unchanged underneath places where the load is absent. Among the six models studied here, those with a glacial ice cover of variable thickness appear best suited for chimney formation. For the latter models, an estimate shows that the ice thickness must be at least 1000 m for the aquifer pressure to exceed the least compressive stress at the base of the caprock.
where s is the sediment matrix density. Note that expression (21) does not involve the porosity. The final 1-D pressure equation thus becomes The 3-D version of the pressure equation has the flow term replaced by ∇((k∕ )∇p) . In this version, sediment compaction is restricted to the vertical direction, while the fluid flow is in 3-D. The pressure Eq. (22) can also be used to model pressure build-up by glacial loading by replacing the righthand side with the rate of change of the load where L(t) is the surface pressure of the load at time t.

Gibson's solution of the pressure equation
The pressure Eq. (22) has an exact solution for a basin in a state of deposition at a constant net rate , where = t . The fluid compressibility c f is assumed to be zero. Furthermore, it is assumed that the porosity is given by the function (1) of effective stress, and that the permeability is given by the function (2)

Poroelastic mechanical equilibrium
The mechanical equilibrium is expressed by where (1) ij is the full stress tensor, b is the bulk density and g is the gravitational acceleration (Wang 2000). The Kronecker delta is 1, if the two indices are equal ( i = j ), otherwise it is 0 ( i ≠ j ). Einstein's summation convention is used in the expression (28), which means that there is a summation over every pair of equal indices. The full stress (1) ij is decomposed into the initial stress (0) ij and the poroelastic stress change ij The initial stress field satisfies the mechanical equilibrium (28), which implies that the poroelastic stress change is given by the mechanical equilibrium without any body forces The poroelastic effective stress ( ′ ij ), which is the stress added to the pore pressure (p), multiplied by Biot's coefficient ( ), causes linear elastic deformations in a poroelastic medium where ij is the strain tensor, is Lamé's first parameter and G is the shear modulus (Wang 2000). The effective stress Definition (31) uses the sign convention that compressive stress is negative. The strain is obtained from the displacement field as where u i is the displacement in the directions i = x, y, z . The spatial directions can also be denoted by x 1 = x , x 2 = y and x 3 = z . The displacement field in the three spatial directions are also denoted, as u = u 1 , v = u 2 and w = u 3 , respectively.

Poroelastic fluid flow
Poroelastic mechanical equilibrium is coupled with fluid flow through the pressure equation (28) (1) ij (1) where S is the storativity, k is a scalar permeability and kk is the bulk strain (Wang 2000). In the case of undrained conditions (where the fluid does not move in the pore space), the pressure Eq. (33) gives which can be integrated to assuming that both the pressure p and the bulk strain are zero at time 0. Inserting the stress-strain relation (31) into the equilibrium Eq. (30) gives the equation which can be expressed with respect to the displacement field by use of the definition of strain (32). In undrained conditions, Eq. (35) can be used to replace the pressure p by the bulk strain kk , which gives an equation in only the strain Therefore, the mechanical equilibrium for the undrained pressure response is expressed by the Eq. (37) from standard linear elasticity, by use of the undrained Lamé parameter Similarly, the bulk modulus for drained conditions, K = + 2G∕3 , becomes K u = u + 2G∕3 = K + 2 ∕S for undrained conditions. It is often advantageous to use the shear modulus rather than the bulk modulus, and the undrained version of the shear modulus is proportional to the undrained bulk modulus in the same way that the shear modulus is proportional to the drained bulk modulus, G u = (3∕2)(1 − 2 )K u ∕(1 + ) . The undrained fluid pressure response to a surface load is obtained from the bulk strain, Eq. (35), where the bulk strain is found by solving Eq. (37) for standard linear elastic deformations using the undrained Lamé parameter u .
show that the double integral becomes 2 . The contribution to the bulk strain from the band from x = −a to x = a is where the double integral becomes The difference between kk = A − B gives the bulk strain (12). The Boussinesq solution also gives the vertical strain as zz = 1 2 kk , (Bower 2009). The strain in the y-direction is approximated by yy = 0 because the two parallel loads in Fig. 9 are constant in the y-direction. The assumption yy = 0 implies that xx = 1 2 kk . The components of the normal strain xx = zz = 1 2 kk and yy = 0 give the stress increase as xx = ( + G + 2 ∕S) kk , and the effective stress increase as � xx = ( + G) kk . Hence, the fluid pressure, and also the stress increase in the x-direction, are proportional to the bulk strain.