Numerical forward modelling of the overpressure build-up in the Cenozoic section of the Central Graben in the North Sea

One of the most widespread hypotheses for the origin of the present-day overpressure in the shale Post-Chalk section in the North Sea is the very rapid sedimentation from Neogene to present day. We tested this hypothesis by the means of numerical forward finite elements modelling and successfully simulated the overpressure build-up during the Cenozoic filling of the North Sea with relatively simple model set-up. Our model shows that overpressure of approximately 28% above hydrostatic developed in the Neogene, while during the Quaternary, it reached up to 36% above hydrostatic. At present day, the predicted onset of overpressure starts at about 800–1000 m below seafloor, while the maximum (magnitude about 1.36 sg, i.e. 36% above the normal hydrostatic pressure) is at approximately 2100 m. This overpressure profile fits reasonably well with data from wells drilled in the Central Graben. The exact magnitude of the overpressure depends on the used assumptions, the model set-up and the values of the input parameters. Especially the dynamic interaction between high sedimentation rates, clay permeability and low horizontal pressure gradient seems to be a key factor in the efficiency of dewatering of saturated clays during burial. The results indicate that, the assumption of horizontal stress isotropy results in nearly no horizontal fluid flow, despite the same magnitude for the vertical and the horizontal permeability. In these conditions, the vertical permeability plays much bigger role than the horizontal one in the effective de-watering of the sediments during burial. Further investigation is needed to explore the role of horizontal pressure gradient in fluid migration in passive sedimentary basins.


Introduction
Important oil provinces such as the North Sea in Europe and the Gulf of Mexico in North America consist of basins where thick sequences of clay-rich sediments (mudstones, shales) were deposited very quickly during the Cenozoic; both regions are characterized by significant overpressure (pore pressure above the hydrostatic) in these formations (Dickinson 1953;Vejbaek 2008).
The mechanism of the origin of the significant overpressure in the Cenozoic section of the Central Graben (in the North Sea) is vividly debated (Osborne and Swarbrick 1997;Holm 1998;Swarbrick et al. 2002;Yang and Aplin 2007;Vejbaek 2008). Most of the authors agree that vertical effective stress and permeability play both a significant role in porosity reduction and in the build-up of abnormal pressures, but the detailed mechanism and the complex dynamic interaction between these factors is not yet fully understood.
In offshore an environment, porous material (clay, sand, etc.), deposited at the bottom of the sea, is fully saturated and the water contained in the pore space is squeezed out when new sediments are deposited on top of the existing ones. This leads to compaction, i.e. reduction of pore space. When the pore water is expelled efficiently, the fluid pressure inside the pores of the sediments remains hydrostatic. Key parameters relating to pore pressure are the grain size and shape, which relate to permeability, and the sedimentary influx to the basin (Catuneanu 2002;Catuneanu et al. 2011). Intensive sediment influxes result in rapid burial and increase of the applied stress due to the weight of the overlaying material. The speed of the de-watering depends (among other factors) 1 3 on the permeability and the applied stress. The rate of deposition may exceed the speed of de-watering of the underlying formations, thus the fluid remains trapped in the pore space, preventing compaction and generating overpressure (pressure larger than the hydrostatic). This mechanism is sometimes called "undercompaction" or "compaction disequilibrium" and is believed to be the main cause of overpressure in sedimentary basins (Osborne and Swarbrick 1997;Swarbrick et al. 2002;Yang and Aplin 2007;Vejbaek 2008).
Thus, both present day pore pressure (called also formation pressure), stress state (magnitude and orientation of vertical and horizontal stresses), strain and rock properties (porosity, permeability, stiffness etc.), of the rock formations are the result of the geological evolution of the sedimentary basin.
The purpose of the present study is to study the dynamic evolution in pore pressure and mechanical properties of the Cenozoic clay-rich sediments in the Central Graben of the North Sea as result of the Post-Chalk fill in of the basin. The focus is, in particular, to study the build-up of abnormal formation pressure in the thick clay/shale formations deposited during the period from Eocene to present day. The sedimentary process is simulated by the means of forward finite elements modelling. Such a numerical model allows to test and, to some extend quantify (taking into account uncertainties), different hypothesis on, for example, permeability magnitude, rates of deposition, sea depth, etc. in the efforts to understand better the evolution of mechanical and fluid properties of the present day sedimentary formations during the process of basin formation.

Geological background
The North Sea, situated between the British islands, Scandinavia, NW Germany, the Netherlands, Belgium and NW France, connects to the North Atlantic through the Norwegian Sea and the English Channel. The Central Graben in the North Sea is part of an aborted trilete rift system (Armour et al. 2004;Ziegler 1992). The rifting was a response to long-running extensional tectonics in the North Atlantic realm, spanning close to 300 million years (Ziegler and Cloetingh 2004). The main rifting events occurred during the Late Carboniferous-Early Permian, the Late Permian-Triassic and peaked during the Late Jurassic (Armour et al. 2004. Following the opening of the North Atlantic during the late Paleocene-earliest Eocene (Mosar et al. 2002;Talwani and Eldholm 1977), a mid-ocean ridge was established in the North Atlantic. This mid-ocean ridge was a source of intra-cratonic stress in most of north Europe (Mosar et al. 2002).
The infill history of a sedimentary basin is controlled by a complex interplay between tectonic evolution and climate. The Cenozoic filling of the North Sea Basin has been investigated intensively since the 1960′s (Gibbard and Lewin 2016). The Post-Chalk sequences deposited in the Danish sector of the North Sea, originated primarily from Scandinavia, Scotland and the Baltic region (Konradi 2005). The most intensive sediment flux from the Fennoscandian shield into the North Sea occurred during the Late Pliocene-Pleistocene and some estimates suggests a maximum of 1000 m Quaternary out of 3500 m thick Cenozoic sediments in the North Sea Basin, north of Dogger Bank (Gibbard and Lewin 2016).
The Central Graben, a part of the Norwegian-Danish basin in the North Sea, is the main depocentre of the Post-Chalk clay-rich sediments. Generally, the Cenozoic succession thins towards both the east and the west, but many factors such as local topography, salt diapirs or domes, etc. have local control of sediments thickness (Evans 2003;Glennie and Underhill 1998).
The deposition in the North Sea Basin was continuously marine throughout the Cenozoic (Schiøler et al. 2007). The two main Cenozoic units in the Norwegian-Danish Basin, the Nordland and the Hordaland Group, consist mainly of shales, with smectite as the dominant clay mineral (Nielsen and Rasmussen 2015). The natural variation of the eustatic sea level during the Cenozoic is believed to be in the order of few hundreds of meters for most of the Cenozoic (Zachos et al. 2001). For example, Rasmussen (2004) reports sea level changes of up to ca. 100 m during the Oligo-Miocene. The present day sea depth varies between 40 m in the south of the Danish Central Graben to 80 m in the north.

Model setup, theory and input data
In this study, we model the evolution in time (from Eocene to present-day) of the poro-elasto-plastic properties of clayrich sediment, from its deposition on the sea floor to a consolidated state, corresponding to approximately 2100 m of burial.
The evolution of material and fluid properties, for example, porosity reduction and overpressure build-up over geological time as a function of permeability and sedimentation rate was investigated by the means of forward finite-element modelling.
The litho-stratigraphic information used to construct the model was derived from the final well report of the exploration well North Jens-1 drilled in the Central Graben. The Post-Chalk sequence at the location is approximately 2000-2100 m thick. It is considered that such thickness is representative for the Cenozoic overburden in the Danish sector of the Central Graben. Seismic data suggest flat topography (i.e. nearly horizontal layers) from seabed to approximately 2000 m depth.

3
The Post-Chalk depositional process was simplified to 3 main stages: Stage 1, corresponding roughly to the period from Eocene to Oligocene (Paleogene), with duration of approximately 33 Million Years, MY; Stage 2, corresponding to the Neogene, from Miocene to Pliocene (approximate duration 21 MY) and Stage 3, Quaternary (Recent), with duration of approximately 2 MY.
Since the focus of our study was the clay-rich section of the Cenozoic sediments, the Paleocene was not included in the present modelling, because, in the Central Graben, it consist mainly of chalk (the Danian Ekofisk formation) and sand-rich (Selandian-Thanetian) sediments.
Obviously, the depositional history of the North Sea Basin is far more complex than what is outlined above. Several phases of inversion and compression have been suggested during the Cenozoic in the North Sea, primarily focused around the Paleocene (Nielsen et al. 2005;Japsen et al. 2014) and during the early Miocene (Rasmussen 2004). As mentioned earlier, the Paleocene is not a part of our model. In addition, it is believed that the magnitude of the Miocene events did not affect significantly the depositional process in the Central Graben, and, at the present, were not included in the simulation of the deposition.
The offshore depositional environment is modelled, assuming a water column of 100 m at the top of the sediments at each stage, which assures a realistic representation of deposition in offshore conditions and fully saturated layers at all times. Although the sea depth in the North Sea was changing during the Cenozoic, we consider a sea depth of 100 m at each stage as a suitable simplification for the purpose of this study.
The Cenozoic shale section is modelled as a simple, 2D column, 500 m wide and approximately 2100 m tall (at present day). A width of "only" 500 m was chosen because it is reasonable to assume that the geological, structural and litho-stratigraphic conditions will be the same across the domain.
The sediments are considered to be at their maximum depth of burial at any time, i.e. no erosion, uplifting or any other tectonic events were included in the simulation. This setup was chosen, since it is considered to be representative for the two major Cenozoic shale groups, the Nordland and the Hordaland Group, in the Central Graben of the North Sea.
In our model, the new sediments are deposited on the top of a pre-existing structure, the so-called base or underburden, represented by a rectangular domain with horizontal surface, with the same width as the overlying column. The base domain can be considered as an analogue of the formations underlying the Post-Chalk sediments, such as the Paleocene deposits (e.g. chalk of Danian age) present in most of the central North Sea. It should be noted that modelling of the chalk deposition and properties is not a focus of this study.
Deposition is modelled as a discrete process, the sediments are deposited as 5 horizontal layers per 1 MY (i.e. one layer is deposited every 200,000 years). After the deposition of each new layer, the undelaying ones are merged (deposition by aggregation). In other words, bedding is not considered in accordance with the assumption with homogeneous and isotropic material. The only horizontal boundaries preserved during the simulation are the boundaries separating the three main geological periods: Stage1, Stage2 and Stage3 in Fig. 1.
The following simplifications and assumptions were used: • All depo layers, in all three stages, are composed by the same material, which is homogeneous and isotropic. We used this simplification, since the main focus is to inves- Fig. 1 Simplified representation of the depositional process, not to scale. Paleogene here refers to the time from Eocene to Oligocene tigate the overpressure build-up as function of the speed of sedimentation and permeability. • The sediments are at their maximum depth of burial at any given time, i.e. no erosion, uplift, or other tectonic events are considered; • The sea depth is 100 m at each stage; • The deposited material is fully saturated with a single phase fluid (water); • The maximum stress is vertical (the weight of the sediments) and the horizontal stress is isotropic and caused by the weight of the overlaying sediments; • No temperature and/or chemical effects (such as diagenesis and/or water interaction with the material) are considered; Information regarding the geological age and the depths at which the respective tops were found along the well paths were derived from final well report of exploration wells drilled in the Central Graben in the North Sea, see Table 1. The depths are rounded to the nearest 10-th or 5-th meter.
The simplified model of the Cenozoic (three main geological periods with respective sediment thicknesses and duration in time) and the mechanical properties (density, porosity etc.) of the materials composing the Cenozoic and the base are reported in the Results Section.

Finite elements model
The dynamic, forward numerical modelling of the evolution in material and fluid properties during the build-up of the 2D sediment column in geological time, outlined in Fig. 1, is performed with the finite element software Elfen. The numerical model is semi-coupled: the mechanical field (Eq. 18) is solved explicitly; the fluid field (seepage, Eq. 19) is solved implicitly and the solutions are coupled at regular time steps. The equations, used in the finite elements modelling, are reported below in the "Theory" Section.

Meshing
At the beginning of the simulation (time 0), before the start of deposition, the finite elements mesh consist of 22 triangular elements and 20 nodes. After the deposition of each layer, new elements are added and the mesh is updated. At the end of the simulation, the mesh consists of 296 elements and 180 nodes.

Initial and boundary conditions
The initial and boundary conditions, the model size and finite elements mesh at the final stage are illustrated in Fig. 2.
At the beginning of the simulation, the pre-existing base is allowed to settle under gravity. Gravity is applied at the top of the sediments, acting downwards. The domain is fixed at the sides and at the bottom, i.e. no deformation is allowed at the bottom of the base and across the vertical sides, while the material can deform inside the domain. Deformation is assessed assuming plain strain conditions (uniaxial compaction). These conditions are considered representative for passive sedimentary basins (such as the North Sea, the Gulf of Mexico etc.), where the major stress is caused only by the weight of the deposited material.
The deposited material is modelled as porous medium fully saturated with water, which can flow both vertically and horizontally within the domain. The fluid cannot drain trough the bottom and no additional fluid supply from outside sources is permitted.

Theory
The modelling simulates the Post-Chalk deposition of a "fresh clay", i.e. very soft, fully saturated material that has not been previously consolidated or compacted, at the bottom of the North Sea from the Eocene to present day (duration approximately 56 million of years).
The material is represented by a modified Cam Clay Model (Roscoe and Burland 1968;Crook et al. 2003Crook et al. , 2006aPeric and Crook 2004;Thornton and Crook 2014).
Only the most basic equations are reported in this section. Detailed description of the computational framework and methods used in the software Elfen can be found, for example, in Thornton and Crook (2014), Crook et al. (2006a) and Peric and Crook (2004).
The bulk volume, V tot , of a porous material is the sum of the volume of the solids, V solids , and the volume of the pore space, V pore : Material properties such as porosity, , specific volume, , and voids ratio, e, are determined by these volumes as follow: and are related as: During mechanical compaction (i.e. no temperature and chemical effects considered), the volume of the solids, V solids is constant, thus the reduction of the total volume is due to reduction of the pore space: ΔV = ΔV pore .
The total volumetric plastic strain (i.e. the irreversible change in volume) during burial is defined as: where V o is the initial bulk (total) volume. From Eqs. 1-4, the volumetric plastic strain can be expressed in terms of porosity or specific volume or voids ratio: The subscript "o" refers to the initial state (at depth/stress/ time 0).
The volumetric elastic strain (i.e. the recoverable change in volume) is described as a function of the mean effective stress, σ′: where K is the bulk modulus. The total volumetric strain is then: The change of the bulk modulus with the mean effective stress P′ is defined as (Thornton and Crook 2014): where φ o is the initial porosity, K o is the initial bulk modulus, P co is the initial pre-consolidation pressure, A is a weighting factor, κ and λ are material constants.
If only density data are available, the porosity is estimated from the bulk density and grain density as follows: In the cases where only compressional sonic log from well data is available, the density can be estimated by the empirical relationship, proposed by Gardner et al. (1974): where V p is the p-wave velocity.

Vertical, horizontal and effective stress
The vertical stress is caused by the weight of the sediments and is estimated as: where ρ i and T i are the density and the thickness of the i-th sediment layer.
In the absence of active tectonic forces, and assuming that the minimum and the maximum horizontal stresses are equal, the horizontal stress is derived from the vertical as: where C eff is the so-called "matrix stress coefficient", or effective stress ratio (Matthews and Kelly 1967), related to the Poisson ratio, μ.
The mean effective stress is defined as: where α is the Biot's coefficient.

Pore (formation) pressure
The hydrostatic pressure (called also normal pressure) of a column of water at given depth, H, is calculated as: where ρ water is the water density, g is the Earth's acceleration (9.81 m/s 2 ) and H is the depth.
In undrained conditions, the pore pressure, called also formation pressure (i.e. the fluid pressure inside the pores of a saturated porous material) changes as function of the change in the mean effective stress, σ′ (Fjaer et al. 2008): K f is the fluid bulk modulus, K fr is the frame bulk modulus and α is the Biot's coefficient.
The formation pressure can be also presented in terms of "Overpressure" (OVP), i.e. the ratio between the estimated formation pressure and the normal (hydrostatic) water pressure, P n at given depth, H: Similar expression is used to compare the pressure P w exercised by a drilling fluid with density ρ inside the wellbore, to the normal pressure P n resulting by a water column The overpressure (OVP) is a dimensionless quantity, often expressed in "sg" (specific gravity). Overpressure of 1 means that the pore pressure, P p , is equal to the hydrostatic pressure, P n . A value of, for example 1.3, means that the pore pressure (P p ) is exceeding the hydrostatic pressure by 30%.

Seepage
In the present simulation, we considered only a mono-phase formation fluid (water), and the fluid flow, Q, at given time and given depth is modelled by the Darcy's law: This unit is often called "specific gravity".
where k is the matrix permeability (independent of viscosity), A is the cross-section area, L is the length of the section, ΔP is the pressure difference at the two ends of the section, and μ is the fluid viscosity.

Evolution in time of a coupled mechanical-seepage field
The mechanical field, representing a porous medium, is expressed as linear momentum balance: The Darcy flow over geological time for a single fluid is represented with a transient equilibrium equation: where: σ′ = σ − αP is the effective stress; P f is the fluid pressure; ρ s and ρ f are the solid and the fluid density, respectively; φ is the porosity; k(φ) is the porosity-dependent permeability; g is the Earth's acceleration; a s is the acceleration of the solid phase; K f is the fluid bulk modulus; K fr is the frame bulk modulus; α is the Biot's coefficient and ε v is the volumetric strain.
In the finite-elements simulation, the mechanical field is solved explicitly, while the seepage field is solved by implicit time integration schemes and the two fields are coupled at given time intervals.

Input data
The main material properties used to model compaction during burial are the porosity in conjunction with permeability, as outlined in the Theory section.

Density and porosity
Porosity and density data are typically obtained from cores and/or petrophysical logs, acquired in wells. Until the advent of shale gas and shale oil exploration in recent years, the main efforts were focused on studying reservoir rock properties, there is a significant lack of data (core samples and/ or petrophysical logs e.g. density, porosity, sonic, pressure measurements) from non-reservoir formations such as a clay-rich overburden. This is also true for the Post-Chalk shale sequences in the Central Graben, especially in the shallow section, from seabed down to approximately 1000 m depth.
From now on, all the depths mentioned in the text are given in meters below seabed.
Since our model geometry and depositional history were constructed mainly based on lithostratigraphic information from well North Jens-1, we used the available petrophysical logs (Gamma Ray, Compressional Sonic, Density) from this well to derive information about formation porosity. Unfortunately, the porosity indicator logs (such as sonic, density) were acquired only starting from approximately 1200 m below seabed. We examined (to the extent to which data were available) petrophysical logs collected in other deep wells drilled in the Central Graben, and we found that a well, drilled in the Norwegian sector of the Central Graben, well 1/2-1, acquired sonic and/or density logs from approximately 540 m below seabed. The data from the Danish well (North Jens 1) were provided by Danish Hydrocarbon Research and Technology Centre, while the data from the Norwegian well (1/2-1) were provided by the Norwegian Petroleum Directorate. As mentioned earlier, the two major Cenozoic units, the Nordland and the Hordaland Group are composed mainly of smectite rich shales (Nielsen et al. 2015) at both well locations. The total thickness of the Cenozoic section, composed by the Nordland and the Hordaland Groups, is approximately 2900 m at the location of the Norwegian well (well 1/2-1), and approximately 2100 m at the Danish well site. The difference is due to the fact, that the former is situated centrally in the Cenozoic depocentre of the Norwegian-Danish basin, whereas the latter is located somewhat more marginally (Evans et al. 2003;Amour et al. 2004;Glennie and Underhill 1998) with respect to the depocentre.
In order to fill in the data gap between seabed and 500-1000 m depth, we reverted to data collected by the Ocean Drilling Program (ODP) project.
Since we are studying the evolution of the burial of freshly deposited clay in a marine environment, we needed well, as deep as possible, penetrating a column dominated by Cenozoic marine muddy sediments at their maximum depth of burial and with limited amounts of quartz or carbonate content.
We therefore identified among the ODP campaigns, an offshore area in the Arctic with no deposition of coarse clastic material and with very limited calcite content due to carbonate producing organisms. This area is located in the Labrador Sea between Greenland and Canada, and its depositional history in the last 10 million years (Arthur et al. 1985) can be considered as a suitable analog for the Late Cenozoic evolution in the Central Graben.
We selected site 646 of ODP Leg 105 campaign (Arthur et al. 1985), where a borehole was drilled to a total depth of 767 m below sea floor (mbsf), and cores were taken from almost the entire section. The core analysis showed lithologies with a low quartz or carbonate content (ca. 5-20%, but mostly < 10%), and predominantly clay, with only a few intervals dominated by silt (Arthur et al. 1985). It is therefore reasonable to assume that the clay content in the sediments from the ODP 646 site is similar to that of the clay rich overburden deposited in the North Sea during the Cenozoic. The p-wave velocity data from site 646 showed also very good agreement with velocity data from North Sea clays (Japsen 2018).
Thus, the data collected at the site 646 of ODP Leg 105 campaign provided us with the needed information from seabed to approximately 800 m depth. It should be noted that the porosity and density measurements from the ODP site were not corrected for compaction nor sea depth (Arthur et al. 1985).
Combining the data from the ODP site 646 and the two Central Graben wells allowed us to create a composite density-depth and porosity-depth profiles for the Cenozoic clay-rich formations from seabed to approximately 2100 m below sea floor. However, the petrophysical logs acquired in the wells drilled in the Central Graben (North Sea), provide much less details and accuracy on, for example, clay content and age, compared to the ODP data.
While clay content, density and porosity of the sediments at the ODP site 646 were measured from cores with great details (Arthur et al. 1985), no similar detailed information was available in the well reports of the Central Graben wells. In order to identify clay-rich intervals at the locations of the North Sea wells, we used the gamma ray log to distinguish between clay/shale and non-clay formations. Only formations with GR > 50 API were considered to be clay/ shale rich and used to compare the formation densities of the two areas, the ODP site 646 and the Central Graben. This method allows mainly qualitative classification; it was not used to derive quantitative clay content. We are fully aware of the uncertainties related to this approach and the respective uncertainty in identifying clay/shale lithology.
Where sonic, but not density log data was available, the density was estimated by Eq. 10. Finally, a composite density profile from seabed to approximately 2100 m depth was created by combining the three data sets.
The porosity of the Central Graben clay/shales was estimated from the bulk density by Eq. 9. A composite porositydepth profile combining the ODP with the Central Graben wells was created in the same way as the composite density log. The composite density and porosity datasets are shown in the Results section. We then used this composite porositydepth profile to obtain the changes in voids ratio/specific volume during burial, with increasing mean effective stress (Eqs. 1-8).

Permeability
Permeability in fine-grained materials as clay, mudstone and shale is notoriously difficult to measure. In addition, significant variations in the clay mineralogy are observed in the North Sea, both stratigraphically and laterally (Nielsen et al. 2015;Anell et al. 2012). Therefore, it is not possible to define a universally correct clay type for our modelling. However, as already mentioned, smectite is the dominant clay mineral in most of the central North Sea for a significant portion of the Cenozoic section (Nielsen et al. 2015). Mondol et al. (2008) performed a number of laboratory measurements on clay mixtures with porosity between 0.1 and 0.55, ranging from pure smectite to pure kaolinite, under increasing vertical loading. Increased vertical loading is a reasonable approximation of the process of burial during the filling in of a sedimentary basin. Neuzil (2019) reports matrix permeabilities from 0.1 μDarcy to 0.1 Darcy for high porosity (0.6-0.8) clay-rich sedimentary formations.
We created a permeability-porosity profile for the clayrich depositional material, combining the Mondol et al. (2008) data for pure smectite in the porosity range 0.1 to 0.5 and the upper boundary of the permeability range, reported by Neuzil (2019) for high porosities (0.6-0.8).
Yang and Aplin (2007) measured the permeability (among other parameters) of 30 samples from deeply buried and, therefore, already consolidated mudstones. We considered these data suitable to model the permeability of the material composing the pre-existing underburden upon which the Cenozoic column is deposited. Yang and Aplin (2007) report that the vertical permeability of mud rocks can be up to 10 times smaller than the horizontal one, due to alignment of the plate-shaped clay particles orthogonally to the vertical stress field. Our modelling does not include directly the shape and the alignment of the particles, thus we used the same values for the permeability in the horizontal (X) and the vertical (Y) direction. This is also in accordance with the assumption of homogeneous and isotropic material, i.e. all properties are equal in all space directions.
The final permeability-porosity profile, together with the values reported by Yang and Aplin (2007), Mondol et al. (2008) and the approximate permeability range for clayrich sediments outlined by Neuzil (2019) are presented in Fig. 3. The red dots mark the values assigned to the deposition material, while the purple dots (Yang and Aplin, 2007) show the values used to characterize the pre-existing material in our model.

Results
As mentioned in the previous chapter, porosity and bulk density data for freshly deposited clay-rich finely grained sediments, were derived from the ODP site 646 in the Labrador Sea.
Samples were collected from seabed to approximately 800 m depth and were analyzed in great detail (Arthur et al. 1985). Thus, good quality measurements are available on both bulk density, porosity and the composition, in terms of sand, clay and silt content.
At this site, the age of the sediments is ranging from Pleistocene (Quaternary) at seabed to Upper Miocene, approximately 8 million years ago (Ma) at depth of approximately 800 m below seabed (Arthur et al. 1985).
The ODP data are shown in Fig. 4. All depths are given in meters below seabed. The bulk density is plotted to the left, the lithological composition, in terms of sand, silt and clay content, is given in the centre of the figure and porosity is given to the right.
The most dominant lithological components are the finegrained clay and silt, accounting together for more than 80% of the material (Arthur et al. 1985). The clay is the most dominant component, typically, above 50%, especially below 400 m depth. The silt content is high, as well, but it decreases with depth and below 400 m is typically less than 50%.
The sand content is generally low. In the uppermost 100 m, the sand content is less than 15%, with just a couple of measurements of 20-25%. The sandiest interval is between 180 and 260 m depth, where a couple of values reached 40 to 50%. Below 260 m and down to approximately 800 m the sand content is very close to 0%.
The porosity of newly (recent) deposited clay at seabed is very high, up to about 75-80%, which is in agreement with values reported by literature (i.e. Tucker 2001). As Fig. 2 shows, there is a relatively well-defined trend of porosity decrease from 80% at seabed to 40% at 800 m, while, between 130 and 350 m, there is an interval with slightly higher porosities/lower bulk densities, deviating a bit from the general trend. As outlined above, the sand content alone Fig. 3 Permeability versus porosity. The red dots indicate the values used to characterize the newly deposited material. The purple dots mark the values used to characterize the pre-existing base (underburden). The dashed lines outline the approximate permeability range for clay-rich sediments, summarized by Neuzil (2019) 1 3 cannot explain this anomaly, since only few samples had a sand content more that 10%. The anomalous interval corresponds to the period of the Late Pliocene (Arthur et al. 1985). However, detailed analysis of the porosity-depth trend of the ODP site is beyond the scope of our study. The lithological composition supports the assumption that the upper Neogene-Holocene sediments of the last 8 million years at the ODP site 646 in the Labrador Sea are composed almost entirely of fine-grained formations, mostly clay and silt.
The shale intervals were identified by the Gamma Ray logs, as described in the Methods section. Data related to detailed lithological composition, as in the case of the ODP site in the Labrador Sea, was not available.
The shale density, obtained from petrophysical logs from the two Central Graben wells, is shown in Fig. 5. At the Danish site (left side of Fig. 5), data were available only below approximately 1200 m, while at the Norwegian site (right side), the petrophysical logs were acquired starting from appr. 550 m depth. The figure illustrates well the Fig. 5 Shale density, from petrophysical logs, at the locations of the Central Graben wells. Left: data from the Danish (DK) well. Right: data from the Norwegian (NO) well. The depth is in meters below seabed. Here, Paleogene refers to the time from Eocene to Oligocene general lack of data in the shallow section (seabed to appr. 500-1000 m depth) of the Cenozoic in the Central Graben.
The density data from the two deep Central Graben wells, well 1/2-1 and North Jens-1, are both between 1.9 and 2.1 g/ cc in the depth interval 1200-1600 m (which corresponds roughly to the Neogene period at the two well sites) and then diverge from each other.
This can be expected, since the distance between the two wells is approximately 170 km and it is very natural to have differences in the local conditions, such as lithology, mineralogy, topography etc., since, as mentioned previously, the two wells are placed differently with respect to the axis of the Central Graben (Evans 2003;Amour et al. 2004;Glennie and Underhill 1998). Also the present day thicknesses of the Quaternary, Neogene and Paleogene sections (taking account for the uncertainties in the determination of depth and geological age) are different, as shown in Fig. 5.
The porosities at the two Central Graben sites were estimated from the density by Eq. 9.
The bulk density and porosity data from the ODP site and the Central Graben data derived from petrophysical logs from well 1/2-1 (Norway) and well North Jens-1 (Denmark) are plotted together on Fig. 6. The ODP site data are given in red, the North Jens-1 data are plotted in blue and the well 1/2-1 data are displayed in black. The figure shows that there is a relatively good match between the different dataset in the overlapping depth intervals (500-800 m and 1200-1600 m). Figure 6 shows clear trends with depth/burial, increasing for the density and decreasing for the porosity. The freshly (most recent) sediments at seabed, are characterized by low bulk densities and very high porosities-up to 80% (or 0.8).
The similarity of the density and porosity data allowed as to construct a composite profile from seabed to approximately 2100 m depth, as follows: from seabed to appr. 600 m-ODP data, from 600 to 1200 m-data from well 1/2-1 and from 1200 to 2100 m-data from well North Jens-1. This composite profile was used in the further modelling.
It should be noted again, that while the density and the porosity from the ODP borehole are measured on core material, the data from the Central Graben deep wells are derived from borehole logs, with an associated increase in uncertainty. In addition, it must be mentioned, that the present-day sea depth is rather different at the three locations. Figure 7 shows the final (composite) porosity-depth profile, used in the modelling, together with a hypothetical average trend (with upper and lower limit) of decreasing porosities from seabed to appr. 2100 m depth, with an upper and a lower boundary.
As it can be seen on Fig. 7, the porosity is rather high at seabed (between 0.70 and 0.80) and decreases to 0.25-0.35 at approximately 2000 m depth. The largest decrease from 0.7-0.8 to approximately 0.40, however, occurs in the uppermost section (0-800 m). The large spread in the porosities, especially in the shallow section, illustrated by the upper and lower bounds of the decreasing porosity trend, show clearly the large uncertainties related to establishing a unique normal compaction trend. Figure 7 shows that the majority of the data points between approximately 850 and 1600 m depth deviate visibly to the right (i.e. towards higher porosity values) of the average trend line. If consider the upper bound of the trend line, again there is an interval between 1200 and 1600 m with porosities higher than the hypothetic upper trend. Below 1600 m the observed data are closer to the trend line. This is most likely an expression of under-compaction, i.e. preservation of porosity, related to inefficient dewatering, and consequently, the pressure of the fluid trapped in the pores could be larger than the hydrostatic.
The trend curves, plotted in Fig. 7, are, of course, affected by uncertainty (which is difficult to quantify), given also the different degree of data accuracy in the three wells. As mentioned above, the ODP data are direct measurements on core material, while the Central Graben wells data were derived from petrophysical logs. It is obvious that these hypothetical trend curves represent the present-day state of compaction in the subsurface, as derived from well data.
No further details with respect to compaction trends were investigated, because the main focus of the present study was to study the dynamic interplay between the mechanical and fluid properties during the burial of fully saturated clayrich sediments over the last approximately 56 MY. In this respect, the purpose of the here presented modelling was not to derive a specific normal compaction trend, i.e. an empirical relation between depth or effective stress and porosity and provide values for the parameters in such a relation.
Mean effective stress, voids ratio, specific volume and volumetric plastic strain were estimated as given in the Theoretical section (Eqs. 1-12).
The material properties used to describe the two materials, the newly (fresh) deposited sediments and the material composing the pre-existing base are summarized in Table 2.
The information (main geological stages, duration in time of the given geological period and the corresponding present day sediment thickness) used to model the depositional process in time was based mostly on the data from the well North Jens-1 in the Danish Sector of the Central Graben. The 2D-column (sketched in Fig. 2) is composed of two major domains: a pre-existing base (or underburden) upon which the Cenozoic clay-rich sediments were deposited; the geometry is summarized in Figs. 1 and 2 in the Model Section. The specific numerical values assigned to the geological age and the respective thicknesses are approximate, since in the well report the age was reported qualitatively, as for example "Early Pliocene" or "Mid to Early Miocene" etc. Table 3 shows the duration (in MY) of each stage, together with its present day-thickness (rounded to the 1 3 closest 10-integer) and the number of discrete layers used in the numerical simulation of the depositional process.
Simplified depositional rates for each geological period were derived from the respective present-day thicknesses Fig. 6 Comparison of density and porosity data from the Labrador (ODP site 646) and the two Central Graben wells: well 1/2-1, Norway (NO) and well North Jens-1, Denmark (DK). Depth in meters below sea floor and time duration in million years (MY). As mentioned earlier, the Paleogene is restricted to Eocene and Oligocene (i.e. the post-chalk sequences).
A graphic illustration of the build-up of the sedimentary column in time is shown in Fig. 8.
As it can be seen in Table 3, the deposition rate was steadily increasing from the Paleocene (Eocene-Oligocene) to present, and the Quaternary rate was approximately 9 times higher than that in the Paleogene and 3 times higher than the Neogene's rate. It must be emphasized again, that the depositional rates etc. are derived from the present day thicknesses. In basin analysis, it is a common practice to use iterative back stripping in order to derive the original thicknesses at specific time and location of the major stratigraphic units. In our study, we used a preliminary forward "dry-run" to quantify roughly    The aim of our study was to produce a generic dynamic model for the origin and the evolution of overpressure in the North Sea, not to replicate in great details the original thicknesses of given sedimentary unit.
To illustrate the build-up of overpressure, a 1D plot of the formation pressure Pp versus depth in the centre of the domain was extracted at the end of each main geological period/stage. The result for the isotropic permeability case, where the vertical and the horizontal permeabilities are the same, PermX = PermY is shown on Fig. 9. The plot shows the overpressure (OVP) along a hypothetical vertical line placed in the middle of the domain with 20 m vertical resolution. The overpressure is defined (Eq. 16) as the ratio between the actual pore pressure, P P , at a given depth and the corresponding hydrostatic (or normal) pressure P n . Thus, an overpressure of 1.5 means that the P P exceeds the P n by 50%. The results for the three case scenarios are presented below. As the Fig. 9 illustrates, the formation pressure was hydrostatic during the Paleogene deposition. Overpressure started to develop during the Neogene deposition, starting at about 5-600 m below seabed and reaching approximately 1.28 (i.e. exceeding the hydrostatic by 28%) at the base of the Paleogene section by the end of the period. The overpressure continued to build-up also during the Quaternary, adding additional about 8% and reaching approximately 1.36% at the end of the simulation. Overpressure developed also in the pre-existing base.
Since no direct measurements of the present-day Cenozoic formation pressure were available, the estimated Quaternary formation pressure was compared to the static pressure of the drilling fluid, used to drill the North Jens 1 (the rightmost side of Fig. 9). The drilling fluid pressure, i.e. the pressure exercised inside the well of a column of drilling fluid (called often mud) with given density (the so called mud weight, MW) can expressed in terms of specific gravity and compared to overpressure (Eq. 16). As a standard procedure in well planning, the drilling fluid pressure is carefully selected to be above the expected formation pressure in order to prevent formation fluid flowing in the well. This, when lacking detailed information or measurements of the actual formation pressure, the drilling fluid pressure can be used as indicator for the expected formation pressure. This is, of course, just an indicator, and should not be confused with the actual formation pressure.
As the figure shows, the drilling fluid pressure was increased from 1.1 to about 1.44 sg at 1150-1200 m below seabed, probably to accommodate expected increase in the formation pressure. This fits rather well with the pore pressure increase, predicted by our model. The drilling pressure was increased to about 1.6 sg below 2300 m, probably in anticipation of additional increase in formation pressure, due to expected hydrocarbon accumulation. Our modelling is focused on the Cenozoic overburden only, and hydrocarbon reservoirs were not considered.
The post-drill analysis of the well Fasan-1 (Jensen 2004), located in the Danish Central Graben, not far away from well North Jens-1, showed that the pore pressure starts to increase at approximately 900-1000 m depth, reaching 1.30-1.35 sg between 1600 and 2300 m depth. These values fit well with the magnitude of the present-day overpressure predicted by our model.
Thus, our approach, with relatively simple representation of the depositional process, succeeded to capture the profile of overpressure build-up in the clay rich Cenozoic Section of the Central Graben. It must be kept in mind that the exact numerical values are the result of the assumptions, the geometry, the model set-up and the material properties used in the forward finite element simulation of the depositional process. Figure 10 shows the estimated porosity, permeability, fluid velocity and overpressure at Present. Please, note that the permeability is plotted versus depth, not versus porosity. The estimated porosity (black curve, at the leftmost panel of Fig. 10), fits reasonably well the observed porosity (red dots) in the very shallow section (between seabed and approximately 200 m) and below 1250 m, while the match is not very good in the middle section (200-1200 m).
One possible explanation for the poor prediction in the middle section (200-1200 m) could be the present-day sea depth. Here we model it as 100 m, which is representative for the North Sea, while the porosity in the interval 0-600 m is derived from the ODP site cores, where the present-day sea depth is more than 2000 m (Arthur et al. 1985). In addition, the deposition history and sedimentation rates at the ODP site (the Labrador Sea) are not exactly the same as the ones in the North Sea. Further modelling efforts and details are needed to improve the prediction of the porosity in the section between 200 and 1200 m depth.
The onset of Overpressure (between 600 and 700 m) corresponds to permeability dropping below 1 microDarcy and vertical fluid velocity dropping below 50 m/MY. The maximum values of the estimated overpressure are reached at about 2100 m, where both the vertical and horizontal fluid velocity approach 0.
Despite the fact that permeability is the same in both vertical and horizontal direction, the vertical fluid velocity is much higher than the horizontal in the section from seabed down to about 1200-1500 m. In the very shallow section, 0-200 m the vertical fluid velocity is between 130 and 100 m/MY and decreases to about 10-20 m/MY at 1500 m depth. The horizontal velocity is close to 0 for the entire section from seabed to bottom domain. In other words, dewatering is much faster and efficient in the vertical direction, then horizontally. This is related to the fact that, fluid flow 1 3 is driven by the pressure gradient ΔP. Since the horizontal pressure gradient is nearly 0 (due to the assumption of horizontal stress isotropy), the horizontal fluid flow will be close to 0 no matter how large the permeability is. On the other hand, the vertical pressure gradient is much larger due to the weight of the sediments resulting in higher fluid velocity. The vertical fluid velocity is approaching the horizontal one for permeabilities below 0.1 microDarcy.
The difference in the fluid velocity is not an artefact of boundary conditions, since this condition requires only nondeformation on the vertical sides of the domain, and does not impose any restrictions on fluid flow across the vertical sides of the sedimentary column. Besides that, the horizontal location of the vertical profile corresponds to the middle of the modelled sedimentary column. This significant difference in the fluid flow in vertical and horizontal direction, despite the isotropic permeability, might have important implications in understanding in more details the role of horizontal stresses and the paths of fluid migration in basin analysis. More investigations are recommended, exploring also different stress regimes, such as non-isotropic horizontal stresses and/or compressional and strike-slip regimes, where the maximum stress is horizontal. Our results indicate that under the used assumptions, values of input parameters and model set-up, despite the same vertical and horizontal permeability the nearly absent horizontal stress gradient prevents effective de-watering of saturated sediments during burial and this contributes to the generation of overpressure in the Cenozoic clay-rich section of the Central Graben. The magnitude of this overpressure is governed by the magnitude of the permeability and the pressure gradient.

Discussions
In our study, we used porosity data derived from the shallow Ocean Drilling Program-cores to fill-in the data gap since the ODP data had a reasonable match with the available well log derived porosities from two different well sites in the North Sea (Fig. 5). It must be pointed out, that the sea depth of the North Sea differs significantly from that of the Labrador Sea (where the ODP site is located). In our simulation, we used a sea depth of 100 m, representative for the North Sea, but not for the Labrador Sea (more than 2000 m). In addition, the deposition history and sedimentation rates are not exactly the same at the two locations. The most prolific porosity reduction seems to take place during the first few hundreds of meters of burial (Fig. 7), where the porosity data are derived from the ODP core measurements, which are not corrected for deposition rate, compaction or sea depth. This, in conjunction with the assumption of horizontal stress isotropy (resulting in very low horizontal pressure gradient) can probably partly explain why our model fails to predict the porosity in the upper section (0-1000 m). Further analysis is needed to address this issue.
Our model uses a porosity-permeability relation (Fig. 3), based on laboratory measurements (Mondol et al. 2008) for a smectite dominated clay lithology, which seems to be the prevalent type in the Cenozoic overburden in the North Sea (Nielsen et al 2015). The poro-perm profile was extrapolated for porosities above 60% within the range outlined by Neuzil (2019), who reports values of approximately 3 mD for clays with a porosity of 80%. The large uncertainties in permeability lead to uncertainty in the predicted pore pressure magnitude.
Our simplified representation of the Cenozoic era does not include all details of the depositional history of the North Sea Basin, such as phases of inversion and compression in the Paleocene, which is anyway not part of the current modelling, (Nielsen et al. 2005;Japsen et al. 2014) and during the early Miocene (Rasmussen 2004). Furthermore, sedimentation rates were kept constant in each major stage. After experimenting with various other scenarios we chose to model the sedimentation process as depositing one discrete layer at each 200,000 years. Too large number of depo layers resulted in solution convergence problems; too few did not represent adequately the deposition process. This simplified depositional model might contribute further to the mismatch between the predicted and observed porosity in the shallow section (approximately, between 200 and 1200 m below sea floor, see Fig. 10).
Although the North Sea Basin is not located at an active tectonic plate margin, several tectonic phases have caused deformation, which could lead to leakage of the overpressure. Furthermore, the overpressure itself may cause deformation such as intra-formational faulting (Cartwright and Lonergan 1996) and fluid expulsion features (Andresen and Clausen 2014). Such features will therefore reduce the actual overpressure, implying that our results could represent a theoretical maximum for overpressure in a thick clay overburden such as the North Sea. These features were not included in the present model.
Despite all simplifications used, our model is able to predict the build-up of the present-day overpressure in the clay rich Cenozoic overburden at the studied location in the North Sea.
It must be noted, that this study focuses solely on the mechanical compaction, and processes such as thermal induced clay mineral transition from smectite to illite, tectonic effects such as fracturing, folding, erosion and uplift, as well as changes in sea depth were not included into this model. As Mondol et al. (2008) demonstrated, significant differences in permeability exists between various clay minerals. Thermal transition takes place at depths of around 2-3 km (Tucker 2001), and since our model mostly deals with ca. 2000 m of clay, this effect was not considered of relevance in the present model. Neither are capillary and/or diffusion effects included in the present modelling. However, our numerical forward approach allows incorporating the above-mentioned features in future modelling.

Conclusions
The results presented here reflect a simplified 2D finite elements forward modelling of the overpressure build-up in the Central Graben caused by the deposition of clay-rich sediments during the last 56 million years from Eocene to present. The results show that the observed present day overpressure in the clay rich overburden originated in the Neogene (Miocene-Pliocene), and increased further during the Quaternary.
The predicted present-day formation pressure is close to normal (i.e. hydrostatic) down to approximately 5-800 m below seabed, a depth interval corresponding to the presentday Quaternary-Upper Neogene litho-stratigraphic section. At approximately 1000 m depth (mid-Neogene section), the pore-pressure exceeds by approximately 10% the hydrostatic and reaches the highest value (approximately 36% above the normal) between 1500 and 2000 m below seabed, in the Eocene section of the litho-stratigraphic column. In our model, the depth at which the overpressure starts to increase (Fig. 9) is corresponding to the observed overpressure below the Mid Miocene Unconformity (Japsen 1999(Japsen , 2018Jensen 2004).
In the absence of direct measurements of the pore pressure in the clay-rich overburden, the estimated pore pressure was compared to the pressure of the drilling fluid used to drill the well North Jens-1. The comparison had shown that the pore pressure predicted by our model matches the profile of the pressure of the used drilling fluid. Our prediction fits well also with the post-drill analysis of the pore pressure in another Central Graben well, Fasan-1 (Jensen 2004). Thus, we believe that, our model captures reasonably well the observed present-day pore pressure profile of the clay-rich Cenozoic overburden in the Central Graben.
We modelled the Cenozoic as composed by three main stages: stage 1-Eocene to Oligocene, stage 2-Neogene and stage 3-Quaternary, with average post-compaction sedimentation rates of 18, 57 and 150 m/MY, respectively.
The time evolution of the formation pressure showed normal pore pressure during the slow Eocene-Oligocene deposition (ca. 33 MY). Overpressure started to develop during the more intensive Neogene sedimentation, reaching approximately 28% above the hydrostatic at the end of the period. The very rapid Quaternary deposition contributed further to the overpressure build-up, leading to pore pressure about 36% above the hydrostatic.
Thus, the Neogene deposition (Miocene-Pliocene, duration of appr. 21 million years) accounts for around 78% of the predicted present day overpressure, while additional about 22% can be attributed to the intensive Quaternary deposition during the last 2 million years.
The model highlights the dynamic interplay between the sedimentation rate, permeability and pressure gradients in the development of the overpressure.
The results show that the predicted overpressure starts to build up with permeabilities below 1 μD. A test with very high permeability magnitudes (data not reported) shows that when the permeability is above 1 μD, abnormal pressure does not develop despite the rapid Neogene and Quaternary deposition.
Under the assumption of isotropic horizontal stress, the vertical permeability plays an important role in the effective de-watering of the sediments and the generation of overpressure, while the horizontal permeability seems to have very little impact. This is in agreement with the Darcy's fluid flow, which depends on both permeability and pressure gradient. It can be reasoned that, in the case of isotropic horizontal stress, the horizontal stress gradient is approaching 0, resulting in extremely inefficient horizontal drainage and leading to overpressures even for high horizontal permeability. These results show that further efforts are needed to investigate also the importance of horizontal stress anisotropy in conjunction to clay permeability as a controlling parameter for generating the overpressure.
Our results and the specific numerical values of the estimated overpressure are, as in any other subsurface modelling, affected by uncertainties, which are difficult to quantify. These uncertainties are coming from the used assumptions; the chosen numerical values for various input parameters, the model set-up, into some extend also the size and geometry of the finite elements mesh, etc.
Despite the inherited uncertainties, we believe that our numerical forward modelling approach can be used as a useful numerical tool to test and quantify different hypotheses on the effect of, for example, depositional history, sedimentation rates, tectonic events, sea level changes, stress anisotropy, etc. on fluid flow and overpressure build up in sedimentary basins.
Funding Danish Hydrocarbon Research and Technology Centre (DHRTC). Grant: LOCRETA TRD.1E.01. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Declaration
Conflicts of interest On behalf of all the co-authors, the corresponding author states that there is no conflict of interest.
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/.