A 3D model for chimney formation in sedimentary basins

This paper introduces a 3D model for chimney formations in tight rocks in sedimentary basins. This is an adaption of a model for hydraulic fracturing in an anisotropic stress field by fluid injection (fracking). The model assumes that a chimney formation is triggered and sourced by overpressure build-up in permeable units, such as reservoirs or aquifers. Cells in the numerical models fracture when the fluid pressure exceeds the least compressive stress and a random rock strength. Chimney growth is represented by chains of cells (branches) that emanate from the base of the cap rock. The branches have an enhanced permeability during ascension, because the fluid pressure in the fracture network is greater than the least compressive stress. When the branches reach the hydrostatic surface, the fluid pressure drops below the fracture pressure and the fracture network closes. The reservoir is drained by the branches in the closed fracture network that reaches the seafloor. The model produces pipe-like structures and chimneys as accumulations of branches that reach the surface. The degree of random rock strength controls how pipe-like the chimneys become. Chimney formation stops when the rate of fluid leakage through the chimneys surpasses the production of excess fluid by the overpressure-building process. A ‘‘low’’ permeability of the chimney branches produces wide chimneys with many branches, and a ‘‘high’’ permeability gives narrow chimneys made of just a few branches. The model is demonstrated in a setup that could be relevant for the chimneys observed in the cap rock over the Utsira aquifer in the North Sea. By using the proposed model, the permeability of such chimneys is estimated to be of the order of 10 μD.


Introduction
The rapid increase in the atmospheric CO 2 concentrations over the last 100 years is likely the reason for current climate change (Bryant, 1997). More than 30 Gt of anthropogenic CO 2 has been emitted globally every year since 2000 (International Energy Agency, 2016). The storage of CO 2 in aquifers or depleted oil and gas reservoirs is considered to be a promising way for reducing CO 2 in the atmosphere. A particular concern with subsurface CO 2 storage, however, is the integrity of the sealing rocks above the reservoir (Bachu, 2008;Benson and Cole, 2008;Bickle, 2009).
Chimneys and pipe structures are vertical fluid-flow escape pathways through sealing rock layers. Such structures appear to be common features in sedimentary basins, based on improved seismic imaging. Most of our knowledge about seismic pipes and chimneys is from interpretations of seismic data. The seismic anomalies that distinguish chimneys from their surrounding sediments have been thoroughly characterised, as for instance by Kartens and Berndt (2015) and Cartwright and Santamarina (2015). The anomalies appear as dimmed or distorted reflections inside the pipe structure compared with the layered reflections on the outside of the pipe. The terms ''pipe'' and ''chimney'' are often used synonymously, although Andresen (2012) used ''pipe'' E-mail address: Magnus.Wangen@ife.no.
for vertical columnar structures with a diameter of less than 300 m and ''chimney'' for vertical structures that are wider and have a more complex shape. Pipes and chimneys have been observed in a number of places around the world. For example, they have been mapped, using seismic imaging, in the cap rock of the Utsira Formation in the North Sea (Kartens and Berndt, 2015;Kartens et al., 2017), in the Faeroe-Shetland Basin offshore UK and in the Namibia Basin (Cartwright et al., 2007;Moss and Cartwright, 2010) and in the Niger Delta (Løseth et al., 2009(Løseth et al., , 2011). An overview of important datasets has been provided by Cartwright and Santamarina (2015).
Chimneys have been interpreted as localised porous channels for fluid flow, and can often be traced to a reservoir formation (Løseth et al., 2009(Løseth et al., , 2011Kartens and Berndt, 2015;Kartens et al., 2017;Räss et al., 2019). These vertical structures serve as leakage pathways through the seal for reservoir fluids, and are therefore important with respect to the seal integrity of the reservoir units used for storing CO 2 . Chimneys form as a result of hydraulic fracturing, when the reservoir pressure exceeds the least compressive stress (Løseth et al., 2009(Løseth et al., , 2011Cartwright and Santamarina, 2015). An existing pipe may have enhanced permeability, and the reservoir overpressure can dissipate by fluid leakage through the pipe (Løseth et al., 2009(Løseth et al., , 2011. Since chimney formation can be driven by high reservoir fluid pressure, they are closely linked to processes of overpressure build-up. Still, only a few attempts to model pipes and chimneys have been made. Some models have been concerned with the possible importance of existing chimneys in CO 2 leakage from aquifers and reservoirs (Kartens et al., 2017;Tasianas et al., 2016). These studies used standard reservoir simulators to model gas leakage scenarios through known chimneys and to constrain their permeability.
Only a few models have dealt with actual chimney formation. One particular concept that has been studied as a possible mechanism for the generation of pipes and chimneys involve so-called ''porosity waves'' (Appold and Nunn, 2002;Räss et al., 2015Räss et al., , 2018. Porosity waves are viscous, porous deformations of areas of increased porosity, which are driven upwards by buoyancy. The porosity-wave process can transport fluids through low-permeability layers in the form of localised ascending regions of high porosity. They form spontaneously, and are self-propagating, high-porosity channels in low-permeability rocks, such as shale. This type of model was initially suggested for magma migration by Mckenzie (1984), although Mckenzie (1987) proposed the same model for the compaction of sediments. Iyer et al. (2017) recently developed a finite-element model for hydrothermal venting in sedimentary basins driven by heat from magmatic intrusions. This approach is similar to the model introduced here, with fractured elements occurring where the fluid pressure exceeds the least compressive stress. These elements have their permeability enhanced, which in turn allows fluids to ascend through rocks that were initially virtually impermeable.
This paper proposes a model for the formation of chimneys and pipe-like structures that builds on a recent model for hydraulic fracturing and damage of low-permeability rocks by fracking operations (Wangen, 2017(Wangen, , 2019. The rock fails when the pore fluid pressure exceeds the least compressive stress and a random rock strength. Chimneys and pipes form by a hydraulic fracturing process that is sourced from an overpressured reservoir (Løseth et al., 2009(Løseth et al., , 2011Kartens and Berndt, 2015;Kartens et al., 2017).
A topic related to chimneys is Pockmarks, which are crater-like depressions observed on the seabed. They are found in a large variety of geological settings on continental margins (Greinert et al., 2010;Kocherla et al., 2015). They may have formed during methane release or directly by the melting of gas hydrates in the shallow sediments (Mazzini et al., 2017). Pockmarks have been mapped at the seafloor and their positions are analysed with statistical methods (Hammer, 2009;Cartwright et al., 2011;Mazzini et al., 2017). Mazzini et al. (2017) found no clustering of the pockmarks.
This article is organised as follows. Firstly, it is explained how the model was developed from one of fracking, which again is based on invasion percolation. Then, how the model deals with in-situ stress anisotropy, how permeability changes in the model, and how the fluid pressure is used for fracture and damage propagation are demonstrated. A reference case is presented, in which the pressure build-up could be from glacial loading. Chimney development based on a reference case is discussed, and how the chimney structure depends on the chimney permeability is demonstrated, before the chimney development is related to the volume balance of the pore fluid.

Model based on concepts from invasion percolation
The model presented here for chimney formation is based on 2D and 3D models for the simulation of hydraulic fracturing, damage and microseismicity (Wangen, 2017(Wangen, , 2019. These models build on a percolation model for hydraulic fracturing, and the associated microseismicity, developed by Norris et al. (2014Norris et al. ( , 2015bNorris et al. ( ,a, 2016. Percolation theory is the study of the connectivity of clusters, and is normally studied numerically on regular grids (Feder, 1989;Stauffer and Aharony, 1992;Sahimi, 1994). Invasion percolation was introduced by Wilkinson and Willemsen (1983) as a model for the slow displacement of a wetting fluid (e.g. water) into a porous medium saturated with a non-wetting fluid (e.g. oil). The displacement process is controlled by the capillary entry pressure of the pore throats. The invading fluid enters the pores with the lowest capillary thresholds, leaving behind a cluster of pores filled with the invading fluid. The invasion process is intermittent, and takes place in bursts. The dynamics of the invasion process and the burst sizes have been studied by, among others, Furuberg et al. (1988Furuberg et al. ( , 1996, Måløy et al. (1992) and Aker et al. (1998). It is the burst dynamics that makes invasion percolation suited to the study of the microseismicity of hydraulic fracturing.
A related model to that of the invasion percolation for the hydraulic fracturing of tight rock should be mentioned. Miller and Nur (2000) suggested one based on ideas of self-organised criticality to model fluid generation and expulsion in tight rocks. This model involves burst dynamics and cluster generation.
The model described here for chimney formation makes use of the fluid pressure in determining which cell will break next. The nearest neighbour cells are connected by transmissibilities (also called bonds; see Fig. 1a), which break when the fluid pressure exceeds the least compressive stress and a random bond strength. The random bond strength has a similar function to the random capillary entry pressure in models of fluid flow in porous media. The bonds in this model break in an intermittent manner. When a bond breaks, the intact cell connected by the bond breaks as well. The newly broken cell becomes permeable and is invaded by fluid. The neighbours of the newly broken cell may also then break, potentially leading to an avalanche of breaking events.

Stress anisotropy
Chimney formation is driven by the fluid pressure. The fluid pressure is determined by solving a pressure equation that is obtained using the finite volume method (see Appendix). Therefore, the pressure is represented at the centre of each cell. A transmissibility (or bond) accounts for the permeability of the rock matrix between the cell centres, the distance between the cell centres and the area of the common interface between the cells (see Appendix).
The least compressive stress is different on bonds with different spatial directions, as shown in Fig. 1b, but otherwise is the same at the same depth. The components of the stress state in the principal system are the least horizontal stress ℎ , the maximum horizontal stress and the vertical stress . The vertical stress is taken to be the greatest principal stress, where: and the numerical grid is aligned with the direction of the principal stress. Both the vertical and horizontal stresses increase with depth, and it is assumed that the vertical stress in the rock is the same as the weight of the overburden: where is the depth measured from the seafloor, is the bulk density of the sediments, is the water density, ℎ is the water depth and is the gravitational acceleration. The densities and the water depth are assumed to be constant, for the sake of simplicity. The effective stress in the horizontal plane is taken to be proportional to the vertical effective stress: where ℎ and are coefficients. A prime denotes the effective stress, which is defined as the stress minus the pore fluid pressure multiplied by the Biot coefficient. The fluid pressure is taken to be hydrostatic in the overburden. In the case of isotropic stress in the horizontal plane, ℎ = , the coefficients ℎ and are equal, and all bonds are subjected to the same least effective compressive stress, ′ ℎ = ℎ ′ .

Chimney permeability
An important rock property is the permeability of the chimney, which results from hydraulic fracturing and damage. The term ''damage'' is used to describe a dense and pervasive microfracture network. It is based on the observation of hydraulic fracturing (fracking) in tight petroleum-rich shales, where as much as 7.5 ⋅ 10 3 m 3 to 11 ⋅ 10 3 m 3 of fluid is injected into the shale over a short time span (typically 6 h) (Montgomery, 2013;Norris et al., 2016). This fine and dense network of microfractures enables the injection of large quantities of fluids. The microfractures provide permeable pore space in the rock when they are open. The microfracture network in oil-and gas-rich shales have to be fine and pervasive to enable the production of large quantities of petroleum from an otherwise tight rock.
The chimney most likely has a different permeability when it is rising towards the seafloor than after it has reached the seafloor. When the chimney ascends towards the seafloor, it has a critical fluid pressure that keeps the fractures open, because the fluid pressure is greater than the least compressive stress. When the chimney reaches the seafloor, the fluid pressure at the top of the chimney becomes hydrostatic. A stationary overpressure profile develops, in which there is a linear decrease in overpressure from the root of the chimney to the surface, and the chimney fractures close. Cells in a closed fracture network are assumed to be more permeable than intact cells with nearly zero permeability. Therefore, the chimney stays permeable after it has formed.
The permeability in a chimney propagating towards the surface is called the open-chimney permeability ( ) and the permeability after it has reached the surface is the closed-chimney permeability ( ). However, a precise value for the open-chimney permeability is not needed, as long as it is sufficiently large to keep the fluid pressure greater than the least compressive stress for the ascending chimney. The closed-chimney permeability, which is less constrained, turns out to be important in relation to the width of the chimney.

Fluid pressure and rock damage
Chimney formation in this model is driven by overpressure build-up. The overpressure is obtained from solution of the pressure equation: where is time, is the porosity, is the system compressibility, is the rock permeability at position , is the fluid viscosity and is a source term.
A broken bond changes the associated cell from an intact cap rock cell to a damaged chimney cell. The chimney cell has an enhanced permeability that allows fluid to invade the newly damaged cell. Fluid flow into a new cell changes the fluid pressure locally, and therefore it is necessary to resolve the pressure equation. It would be too timeconsuming to solve this for the fluid pressure in the entire numerical domain each time a bond and cell is damaged. To save computation time, the cells in the cap rock are assumed to be impermeable, except for the nearest neighbours to the reservoir or chimney cells. Therefore, the numerical solution mainly covers the permeable cells -those in the reservoir and the chimney. The nearest neighbour cells to the reservoir and chimney cells are assumed to be low-permeability cells with zero overpressure. These intact cap rock cells are assigned zero overpressure by the boundary conditions. The numerical pressure equation requires the permeability of the bond that connects two neighbour cells, and the bond permeability is obtained as the harmonic mean of the two neighbour-cell permeabilities.
Breaking a bond changes the bond permeability, which makes the model nonlinear. Therefore, the stiffness matrix in the linear equation system must be recomputed each time a bond breaks. The discrete pressure equation not only has to be resolved after each broking bond and cell event, but also the stiffness matrix has to be regenerated.

Reference case for the Utsira Formation
Here, a case is presented that could be relevant to the chimneys observed in the cap rock of the Utsira Formation in the North Sea, close to the CO 2 injection site (Kartens and Berndt, 2015;Kartens et al., 2017). There are two questions related to these chimneys -how they formed and what the present-day permeability might be.
The Utsira Formation in the North Sea comprises a more than 400 km long Pliocene sand (Helland, 2011). The average depth from the sealevel to the Utsira Formation is 900 m, and its thickness ranges from 50 m to 350 m. The reservoir permeability is roughly 1 D (1 ⋅ 10 −12 m 2 ) (Helland, 2011). It has served as a storage site for the CO 2 captured from the natural gas produced in the Sleipner Field, with about 14 Mt of supercritical CO 2 having been injected between 1996 and 2013 (Cavanagh and Haszeldine, 2014). The pressure and temperature conditions at the injection site are just within the region for supercritical CO 2 (Helland, 2011). The Utsira Formation reservoir is covered by a layer of fine-grained Quaternary sediments, mostly clay. Above the cap rock, there is seawater, ranging from 100 m to 200 m deep (Helland, 2011).  The geometry of the simulation model is shown in Fig. 2. The numerical representation covers a lateral extent of 2 km × 2 km and a height of 1 km. A flat seafloor is at the top of the model, at = 0. The top of the reservoir is sin-shaped, and oscillates between 700 m and 900 m deep. The colours in Fig. 2 indicate the depth of the top of the reservoir unit. The cap rock lies between the top of the reservoir and the seafloor. The base of the reservoir, at 1000 m, is also the base of the model.

Pressure build-up in the reservoir
It is assumed that chimneys are formed as a result of pressure buildup in the reservoir unit. The overpressure in the numerical model is created by a source term in the reservoir cells. Glacial loading has been suggested as a possible process for the pressure build-up that triggered the chimneys observed in the cap rock (Kartens and Berndt, 2015). The entire North Sea was covered by ice during the Last Glacial Maximum (Lambeck et al., 2000;Mangerud et al., 2011). Whilst there are no data for the palaeo-ice thickness in the North Sea, models have indicated that the ice was considerably thicker than 1000 m over the Utsira Field during the Last Glacial Maximum (Siegert et al., 2001;Kleman et al., 2013).
Instead of modelling the possible pressure build-up directly from the ice loading, which is a large topic in itself, a simple source term is applied in the pressure equation. The overpressure build-up is not modelled from zero, but from an initial value close to the fracture pressure of the cap rock. The critical fluid pressure necessary to break the weakest bond in the cap rock is denoted . The initial overpressure is set to 1 = 1 , with 1 = 0.9. The pressure build-up takes place over a time span of 0 , and the reservoir overpressure would have ended at 2 = 2 , with 2 = 1.2 in the case of no chimney formation. The pressure equation (13), with negligible pressure gradients in the reservoir, is: where 0 is the source term in the reservoir unit. A linear pressure increase from 1 = 1 to 2 = 2 over the time span 0 gives the source term: The source term (6) is used in pressure equations (4) to simulate the geological time interval during which the chimney formation took place. Because the source term does not model one specific process for pressure build-up, it can represent the effects of a number of different overpressure-generating processes. The simple source term (6) makes it possible to study chimney formation without the complications of a full model for the specific process of pressure generation on a geological time scale. In particular, the complicated geological history of repeated glacial loading and unloading over several hundreds of thousands of years is not the topic of this article. An eventual pressure build-up in the tight cap rock from, for instance, glacial loading can be ignored. The cap rock is taken to be virtually impermeable.

Chimney formation
The chimney formation for the geometry of Fig. 2 is shown in Fig. 3. Table 1 provides a listing of the simulation data used. Fig. 3 shows the development of the chimney at three times, = 100, = 200 and = 400 years, where the time span is 0 = 400 years. Chimney formation starts with a single chain of cells (branch) that ascends towards the surface. The chains grow through the cap rock from the highest parts of the reservoir rock because the lateral compressive stress on the bonds is the least there. Even though the bonds have random bond strengths, the bond strength is, in this case, not great enough to dominate the difference in compressive stress between neighbour cells in the vertical direction. The random bond strength is a means of representing the effect of the often strong heterogeneity of the rock. An important observation is that the effective compressive stress is proportional to the depth from the seafloor, but the random bond strength is independent of depth. The maximum bond strength is = 1 MPa, = 1 MPa and = 0.5 MPa, and the bonds are weaker in the vertical than in the horizontal direction, which favours the vertical propagation of the branches.
After the first chains of cells reach the surface, new chains nucleate and propagate towards the surface. The permeability in the chimneys is important. Before a chain of cells reaches the surface, they must have high permeability to ensure that the fluid pressure is greater than the least compressive stress in the entire chain. A fluid pressure that exceeds the least compressive stress opens the fractured pore space, which again produces an enhanced permeability. When the chains reach the hydrostatic seafloor, the fluid pressure drops below the least compressive stress and the cells in the branch close. The permeability of the hydraulically-fractured but closed cells is still higher than the low permeability of the intact sealing rock matrix. Therefore, the reservoir overpressure starts to leak through the chimney structure that has reached the surface. Fig. 3 shows the overpressure in the reservoir and chimney cells. There is a linear pressure decrease through the chimney, from the base at the overpressured reservoir to the hydrostatic surface.
The permeability of the reservoir and the cap rock are well constrained. The reservoir permeability is = 10 −12 m 2 and the cap rock is virtually impermeable, with a permeability estimated to be less than 10 −19 m 2 , (Helland, 2011). The permeability of the hydraulicallyfractured and damaged chimney cells is taken to be constant because of the lack of a better permeability model. In the simulation case shown in Fig. 3, the open hydraulically-fractured cells have a permeability of = 10 −12 m 2 and when they close, the permeability drops to = 10 −17 m 2 . A larger open permeability than = 10 −12 m 2 will not change the results noticeably, because this permeability is sufficiently large for the pressure drop inside a rising chimney path to be negligible. Decreasing the permeability will at some point lead to pressure drop inside the rising chimney that will eventually stop the chimney for ascending towards the surface. The chains that have reached the surface, and that drain the reservoir of excess fluid, do not refracture. They typically experience a linear overpressure decrease from the reservoir to the surface. On the other hand, new chains of cells continue to nucleate and grow towards the surface.

Branches in a chimney
This model of chimney formation produces a pipe-like structure, in terms of fracture branches emanating from the top of the reservoir. These branches are numerically represented as chains of damaged (fractured) cells. The model differs from the percolation invasion models of Norris et al. (2014Norris et al. ( , 2015b and Wangen (2017Wangen ( , 2019 in not having a loopless fracture network. In other words, the branches in the fracture network are allowed to connect. The reason for not connecting branches in the fracking models is that the pressure in two branches that are close is almost the same. This is not the case here, where there are usually different pressures in a branch (chain of cells) that is ascending and a branch that has reached the surface. An ascending branch is short-circuited once it connects to a branch that is already in contact with the surface, and it stops.
The permeability of damaged cells in the chimney ( ), once they are connected to the surface, is an important parameter. In the case where this permeability is high, there is an implication that overpressure build-up stops, since more fluid leaks from the reservoir than is generated by the source term. On the other hand, low permeability in the damaged cells in the chimney implies that the pipe-like structures will not drain the reservoir sufficiently to stop pressure build-up, and so the process of chimney formation will continue. These effects are examined in Fig. 4, which shows chimney formation for the three different values of closed-chimney permeability, = 10 −16 m 2 , = 10 −15 m 2 and = 10 −14 m 2 . With increasing , it can be seen that progressively fewer branches are needed to drain the pressure build-up in the reservoir. In Fig. 4a, it is shown that a closed-chimney permeability as high as = 10 −14 m 2 will produce only one branch in the chimney. Decreasing by one order of magnitude gives two branches. When is decreased by yet another order of magnitude, is sufficiently lowered to avoid a too rapid leakage from the reservoir, and several branches develop. The value = 10 −17 m 2 seems to be reasonable for a case such as Utsira (see Fig. 3), where several chimneys have been observed with widths of ca. 300 m.
It is possible to estimate the number of non-interacting and leaking branches that would be necessary to drain the reservoir at the same rate as the source term (6) fills the reservoir. This condition can be expressed as: where is the bulk volume of the reservoir, is the -cross-section of a cell and 0 is the distance from the shallowest part of the reservoir to the seafloor. The left side of Eq. (7) is the time-rate of volume production of excess fluid in the reservoir and the right side is the time-rate of volume leakage through the chimney. Inserting the source term (6) into Eq. (7) gives the number of straight vertical chains of cells reaching the surface: It can be seen that the number of vertical branches is inversely proportional to the permeability of a closed chimney ( ). If < 1 implies that the numerical resolution with an -cross-section is too large for the given permeability , or vice versa. On the other hand, if is larger than the number of cells that cover the surface, the permeability is too low for the chimney to drain the reservoir at the same rate as the pressure build-up.
The data from Table 1 indicates that from Eq. (8) is = 70. The number of fractured branches in Fig. 3c that reached the surface is 64. Even though the estimate is not exact, it is still quite useful as an estimation of the size of the chimney.
The Eq. (8) shows that the area of the chimney cross section, = , does not depend on the grid size. The final chimney cross section is controlled by the overpressure generating process. The number of the chimney paths increases with increasingly finer grid size. If there is a lower or upper limits on the chimney paths are difficult to tell from the modelling. More observations are needed to improve this aspect of the modelling.
The random rock strength adds disorder to the branches that constitute the chimneys. Increasing the random rock strength causes the chimney branches to be dominated by the random rock strength and, to a lesser degree, be controlled by the compressive stress that decrease with depth. A reduction in the disorder produces increasingly straighter chimneys, as demonstrated in Fig. 5. These three cases have isotropic random bond strengths of 0 = 0.1 MPa, 0 = 0.2 MPa and 0 = 0.5 MPa, respectively.
The branches that constitute the chimney start from the shallowest parts, at the base of the cap rock, because the seafloor is flat and there is no lateral dependence of the effective compressive stress. In real cases, there are most likely lateral differences in the compressive stress, which may be important for where chimneys prefer to develop. In the case of a nearly flat base of the cap rock, chimneys will develop locally in the cap rock where the compressive stress is the least. Fig. 6 shows the pore fluid volume balance of the permeable parts of the system -the reservoir and chimneys. The green line represents the volume of excess fluid that is generated in the reservoir unit and that drives the overpressure build-up. The red curve indicates the amount of reservoir fluid that has leaked out. This is the sum of what has leaked to the surface and what has leaked into the cap rock. The blue curve illustrates the amount of excess fluid that remains in the reservoir. The black dots represent the sum of the red and blue curves at each time step of the simulation. It can be seen that the volume conservation is excellent.

Volume balance for the reservoir and chimneys
Nearly all the fluid generated in the reservoir stays in the reservoir until the first branch of the damaged cells reaches the surface. That takes place at approximately 80 years, at which point the leakage increases and the volume of the excess fluid in the reservoir starts to decrease. The leakage of excess fluid continues to increase with an increasing number of branches in the chimneys. The chimney-forming process ends when the pressure build-up stops because the leakage through the chimneys is greater than the generation of excess fluid in the reservoir.
The average reservoir pressure is shown in Fig. 7 for pressure buildup over the time spans 0 = 200 years, 0 = 400 years and 0 = 800 years.
It increases until the first chimney branches reach the seafloor. From then on, there is a slight decrease in the average reservoir overpressure.

Conclusions
A physical model has been proposed here for the formation of the chimneys and pipe-like structures commonly observed in tight rocks in sedimentary basins. The chimney model is an extension of a model for hydraulic fracturing of tight rock by fluid injection, and builds on concepts of invasion percolation. It was assumed that chimneys are fluid leakage structures that are triggered by overpressure buildup in aquifers or reservoirs. The dynamics of the chimney formation is therefore closely linked to overpressure generation. A specific case study was presented, which could be relevant for the chimneys observed in the cap rock of the Utsira Formation in the North Sea, where chimney formation may have been driven by overpressure build-up. The pressure build-up could have been by glacial loading. A solution procedure for the pressure equation was suggested, which restricts the numerical solution to the permeable part of the system (the reservoir and the chimney). A source term was devised for the pressure equation, which allows for the study of the actual time interval during which chimney formation takes place. Chimney formation starts when the overpressure exceeds the least compressive stress and a random rock strength. This model produces chimneys that develop as branches of fractured cells emanating from the base of the cap rock. A chimney accumulates fracture branches, thus growing in width over time. Chimney growth ends when fluid leakage through the fracture branches keeps pace with the production of excess fluid, the process that is responsible for the pressure build-up. The fracture branches that constitute the chimneys are of two types. The first type, ascending branches, have a pore fluid pressure that exceeds the least compressive stress, and these branches have high permeability because the fracture network is open. The second type of branches are those that have reached the hydrostatic surface. These have a closed fracture network because the fluid pressure is less than the least compressive stress, and so the permeability is reduced. It was shown how this closed-chimney permeability controls how fast the reservoir will be drained by the chimneys, and how wide the chimneys become. An estimate was provided for the maximum number of branches in a chimney, which is a useful result. The case study suggested that the permeability of the branches in a chimney are of the order of 10 −17 m 2 . Whilst the model was demonstrated with an overpressured reservoir as the source for the chimneys, it could also apply to pipe-like structures driven by melting gas hydrates. The layered structure of sediments of different composition and age plays a role when considering the random strength of the rocks. This paper explores only the simplest choices.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
Parts of this work were supported by the Research Council of Norway through Project 280567 ''CO2-Paths''. The author is grateful to Christian Berndt and Jens Karstens for their constructive comments on the manuscript.

Appendix A. Pressure equation
The pressure equation is derived from conservation of fluid mass in a similar way as done by Wangen (2019), where mass conservation is expressed as: and where is the fluid density, is the porosity, is the Darcy flux and ( ) is source term for fluid generation. The source term is zero, except for the cells in the reservoir unit. Expression (9) can be transformed to a pressure equation by introduction of the effective compressibility: where fluid pressure is , and introduction of Darcy's law: The fluid pressure is replaced by the overpressure which is the fluid pressure minus the hydrostatic fluid pressure ℎ . The equation for overpressure becomes where the Darcy flux is written without the gravity term. The permeability is dependent on whether the position is in the reservoir, the chimneys or elsewhere in the cap rock. The pressure equation (13) is solved numerically with the finite volume method, where the pressure is represented at the centre at each cell in the grid. The numerical pressure equation becomes after volume integration over each cell. The volume integral is replaced by a surface integral by use of Gauss' theorem. The subscript denotes cell and its properties, such as its compressibility , porosity , source term and cell volume . The properties of the connection between two neighbour cells and is denoted with the double subscript . These properties are the distance between cell centres , the interface area , and the permeability between cells .

Appendix B. Computer code availability
The computer code used in this study is available at the web-address I (see Table 2). The code makes use of a library (the ABC-library), which can be downloaded from web-address II in Table 2. The chimney-code and the ABC-library are provided under GNU General Public License. The chimney code makes also use of the linear equation solver cholmod from SuiteSparse-4.2.1, which can be downloaded from web-address III in Table 2.