Impact of uncertainties associated with the choice of the yield stress on the prediction of subsurface reservoir compaction: A field study

The compaction of subsurface chalk during hydrocarbon production leads to technical and societal issues related to well trajectory design, wellbore collapse, seafloor subsidence, and seismic events. Predicting reservoir deformation over time requires well-calibrated constitutive equations to capture the viscoplastic response of chalk according to several parameters including porosity, age, water saturation, and rate of applied stress or strain. The experimental data used to extrapolate trends between e.g. , porosity and yield stress, show a non- negligible data scattering, thereby raising questions about the reliability of the predicted compaction. The present study aims at assessing how the uncertainty associated with the determination of the yield stress man- ifests itself at the field scale. To this purpose, one-dimensional simulations are performed on four producing hydrocarbon reservoirs from the Danish North Sea and quality-checked against seafloor subsidence measure- ments. Two simulations are run per case study by considering in the constitutive equations the initial (conser-vative approach) and final (optimistic approach) yield stress values, which delimit the transition between the elastic and plastic regime of chalk. The results indicate that the accumulated strain and the seafloor subsidence vary by a factor of up to 73% and 50%, respectively between the conservative and optimistic approaches. These differences in the compaction behaviour are influenced by three parameters, the initial porosity, the virgin stress condition, and the stress history of the reservoir. This study demonstrates that the uncertainty analysis of compaction simulation studies is essential for estimation of compaction drive of hydrocarbon fields and to identify the locations potentially subject to failure in the context of hydrocarbon production and field aban- donment as well as CO 2 and thermal storage projects.


Introduction
The high compressibility and reactivity of chalk in the presence of water entail significant field-scale deformations in producing reservoirs that in turn have important economic and societal impacts. 1 For instance, the platform rebuild in the Tyra field due to seafloor subsidence has caused the temporary cessation of hydrocarbon activities. 2,3 In the Groningen gas field, earthquakes have forced authorities to drastically reduce production for several years. 4 Wellbore instability and change in rock physics are other technical issues related to compaction that jeopardise the economic viability of a field and impede the application of advanced reservoir monitoring and stimulation technics such as 4-D seismic, hydraulic fracturing, and radial jet drilling. 5,6 Geomechanical models for chalk have been developed during the last three decades to predict the distribution and amount of strain in reservoirs so that suitable field development plans can be implemented in a safe manner. [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] To this purpose, a yield function is commonly utilised to reconstruct the yield surface of chalk in mean-deviatoric stress (p ′ m -q) plot and is integrated in the constitutive models. The yield surface consists of a shear failure line and end cap that delimit the three main regimes of deformation in the stress space, namely the elastic, shear and plastic regime. The representativeness of the yield surface to capture the deformation behaviour of reservoir rocks is thus of primary importance dictating the quality of the simulation outcomes. The locus of the yield surface in stress space is however not fixed and depends on several factors that are introduced here. The early constitutive models were mainly focused on transferring the laboratory data collected at the laboratory strain rate (0.1%/hr) to the geological rate of deformation (2.10 -10 %/hr) 8,[22][23][24] and on capturing the strength softening of chalk when water replaces oil in the pore space, a phenomenon referred to as the water weakening effect. [25][26][27][28] Advanced experimental studies have then demonstrated the temperature-dependent role of the ion additional electrostatic repulsive force between calcite grains promoting grain-to-grain slippage, a loss in rock cohesion, and pore collapse. 32,33 Thus, more sophisticated constitutive models were built over the past years in an attempt to capture the effects of fluid-rock interactions on chalk geomechanics at in situ conditions and to ultimately enhance the accuracy of predicted deformation. 34,35 Sensitivity and uncertainty analysis studies of the outcomes of geomechanical models (e.g. the amount and distribution of viscoplastic strain and the timing at which it occurs in subsurface) that are related to the accuracy of the experimental data integrated in the constitutive equations have, however, received less attention. [36][37][38][39][40][41] These few works have emphasised the predominant influence of the value of i) the stress at the onset of plastic deformation i. e., the hydrostatic pore collapse or yield stress, ii) the hardening coefficient controlling the shape of the stress-strain curve in the plastic domain, iii) the in situ stress path, and to a lesser extent iv) the elastic modulus and Poisson's ratio. For instance, rising the yield stress value of a carbonate rock in a gas reservoir (Malaysia) by 50% reduces the total seafloor subsidence by 20%. 39 The potential errors of the model outcomes mainly derive from the extrapolated functions present in the constitutive equations and used to capture the trends observed in laboratory. The experimental data, for instance the porosity-dependency of the hydrostatic yield stress for subsurface chalk, 19,42,43 show typically a non-negligible data scattering, thereby raising questions about the representativeness of reported correlation function and thus, about the reliability of the predicted compaction. Such data scattering and associated uncertainty are well-known and caused by the compilation of laboratory data from various hydrocarbon fields. The depositional and diagenetic history specific to each reservoir modifies the texture, mineralogy, and degree of cementation of rocks. Previous laboratory studies have demonstrated the impacts of these three lithological parameters on the geomechanical properties of chalk with a similar porosity, thereby explaining the observed data scattering. 29,30,[44][45][46][47] An alternative is to collect data in one single reservoir to develop a customised constitutive model. However, such dataset is not available in the literature and compiling information from various fields remains the most practical approach.
Identifying the most suitable methodology to quantify the model uncertainty associated with the extrapolated trends derived from experimental data requires an assessment on how these data were collected. Amongst the multiple geomechanical properties needed for compaction simulations, the hydrostatic yield stress value is one of the most important variables as it dictates the stress conditions at which the rock starts to deform plastically. The transition from the elastic to plastic regime is not instantaneous in chalk taking place at one given finite stress value. Contrary, under hydrostatic conditions, it takes place over a stress interval delimited by the initial (σ h0.in. ) and final (σ h0.fin. ) hydrostatic yield stresses and referred to as the transition phase (Fig. 1).
During triaxial tests, σ h0.in. defines the first deviation of the measured stress-strain curve from the interpreted linear elastic trend, whereas σ h0.fin. represents the stress at the onset of linear plastic deformation. The hydrostatic yield stress value assigned to a tested specimen is referred to as the representative, pragmatic, or typical yield stress (σ h0.rep. ) and is located within the transition phase. The method used to estimate σ h0.rep.
varies amongst experimental studies (Fig. 1). While 43 have calculated the yield stress when the ratio Δσ/Δε in the elastic domain (i.e. the elastic modulus) is reduced by 50%, 48 have assigned the yield stress value that corresponds to the strain at which the interpolated linear elastic and plastic trends intercept. The σ h0.rep. value has also been estimated using the orthogonal projection of the intersection point between the linear elastic and plastic trends to the tangent of stress-strain curve. 42 Thus, the estimation of σ h0.rep. in the transition zone is arbitrary and depends on the methodology chosen.
In this context, the present study aims at assessing quantitatively how the uncertainty associated with the determination of the yield stress of subsurface chalk manifests itself at the field scale. To this purpose, one conceptual model and four hydrocarbon fields from the southern part of the Danish North Sea (Dan, Halfdan, Gorm, and Kraka fields) are investigated. Two 1-D compaction simulations are conducted per case study by considering alternatively the σ h0.in. and σ h0.fin. values in the equations of a strain-rate dependent constitutive model. Note that the computed strain of a model using the σ h0.rep. value estimated from one of the three methods shown in Fig. 1 lies within the strain bounds defined by these two simulations. While the conceptual models assist the interpretation of the field data, seafloor subsidence data collected in the studied fields are used to ensure that the simulation outcomes are realistic. Then, the amount and distribution of strain, and its evolution with time are systematically compared for each case study.

Geological setting and production scheme
Located in the southern part of the Danish Central Graben, the four studied hydrocarbon fields include the Maastrichtian Tor Fm. that is composed of five units (M1 to M5 units) and the overlying Danian Ekofisk Fm. consisting of the Danian Porous (D1) and Danian Tight (D2) units. 49 Each field displays specific characteristics, which are described below (Fig. 2). The main Dan field is a 5-km in diameter anticline structure formed by salt tectonics. A major, NE-SW oriented, fault divides the field into two roughly equal-sized blocks, the block A in the NW area and the southeastern block B. This field was discovered in 1971 at a depth of 1850 m and put on stream in 1972. From beginning of production in 1972 to 1988, the main recovery mechanism was pressure depletion. 50 A water injection pilot was first initiated between 1988 and 1991 and implemented at the field scale in 1995. 51 The Maastrichtian M1 unit in the uppermost part of the Tor Fm. represents the main reservoir unit with a second reservoir located in the D1 unit at the top of Fig. 1. Typical stress-strain curve of chalk during hydrostatic loading showing i) the successive four phases during the transition from the elastic to plastic regime and ii) the initial (green dot) and final yield stress (red dot). A zoom at the transition phase displays the linear trend and bulk modulus in the elastic (K e , green dashed line) and plastic domain (K p , red dashed line) as well as the three methods used to estimate the representative yield stress σ h0.rep.,1 , σ h0.rep.,2 , and σ h0.rep.,3 based on 48,43, and 42 respectively. Note that K e /2 represents half of K e and is used by 43 to calculate σ h0.rep. . Note also that the non-linear trend of the stress-strain curve at the onset of loading is caused by the closure of micro-cracks in the rock fabric. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) the Ekofisk Fm. Oil-bearing intervals consist of autochthonous and clean chalk with less than 5 wt% insoluble content. 52 Impure chalk deposits from the D2 unit displaying poor reservoir properties are present in the four fields studied.
The main Halfdan field located in the northwestern area of the Dan field was discovered in late 1998 at a depth of 2100 m and put on stream in 1999. The field extends over ca. 7 km along a NW-SE direction (Fig. 2). No structural trap is present, and hydrocarbon reserves accumulate as a combination of stratigraphic and dynamic traps in which oil and gas flow slowly at a geological time scale towards the Dan field due to the low permeability of the rock. 53 The field is initially developed by kms-long parallel, sub-horizontal wells alternating water injectors and producers to stimulate production. 54 The reservoir rock consists of the M1 unit with a water saturation above 40% and the D1 unit. The Maastrichtian and Danian chalk shows similar characteristics in term of lithology and stratigraphic subdivision with the deposits of the Dan field. Matrix porosity and permeability average 25-30% and 0.5-2 mD, respectively.
The Gorm field is a 3.5-km in diameter anticline structure discovered in 1971 at a depth of 2100 m (Fig. 2). As in the Dan field, a major NNE-SSW fault divides the structure into two blocks. Upper Cretaceous and Paleogene, syn-depositional faulting led to significant thickness variations across the field. The production started in 1981 with pressure depletion as the main mechanism for recovery until 1989. 50 A water injection program was initiated in 1989 and full field scale waterflooding of the reservoir was achieved in 1995. The reservoir consists of fractured clean chalk in the M1 unit with a subordinate reservoir in the D1 unit. Matrix permeability is less than 5 mD and porosity is as high as 45% in large part of the field and deteriorates down flank.
The Kraka field consists of a 8 × 5 km, anticline structure that formed by salt tectonic 55 (Fig. 2). The field was discovered in 1966 at 1800 m depth and production starts in 1991. The main recovery mechanism is pressure depletion and no waterflooding program has been implemented in the field. Contrarily to the other studied fields, the main reservoir is located in the D1 unit of the Ekofisk Fm. This unit consists of clean pelagic chalk with an average porosity ranging between 25% and 40% and with the best reservoir properties identified toward the crest and southeastern flank of the anticline structure. 56,57 A subordinate reservoir unit is found in the M1 unit of the Maastrichtian Tor Fm that consists of clean pelagic chalk.

Subsidence data
Seafloor subsidence data were collected in platforms located in the four studied fields over a period ranging between 5 years in the Halfdan field and 28 years in the Gorm field (Fig. 2). Vertical movements of the platforms were monitored by a GNSS (global navigation satellite system), which consists of satellites transmitting positioning and timing data from a GNSS emitter located on the platform to a receiver located onshore. The elevation data refer to mean sea level. Importantly, the GNSS monitoring surveys were initiated 4 to 16 years after start of production. Consequently, the total seafloor subsidence modelled prior to the start of GNSS survey is added to the measured subsidence values for comparison purposes between observed and simulated data.

Stress conditions
Initial pore pressure and its evolution with time were calculated every year from a total of 691, 663, 355 and 136 repeat formation tests (RFT) data collected in the Dan, Halfdan, Gorm and Kraka fields, respectively. The RFT data were screened to reflect pore pressure in the central part of the field, close to the platforms monitored by GNSS. A margin of error at 90% confidence level was estimated for each annual pore pressure estimation. During calculation, gas, oil and water pressure gradients of 2.17 kPa/m, 6.79 to 7.47 kPa/m, and 9.84 to 10.2 kPa/m, respectively, were used depending on the field studied. The gas-oil and oil-water contacts were assumed horizontal during stress calculation. Previous studies have clearly shown that the fluid contacts are tilted toward the southeast and associated with a regional lateral gradient ranging between 35 and 100 kPa/km. 53,[58][59][60][61] Given that the crest of the Kraka, Gorm and Dan fields is 3 to 5 km in diameter and the 7 km-long Halfdan field, the tilted fluid contacts imply a maximum of 105-700 kPa pressure variation across the fields. The latter values are negligible compared to the reservoir pore pressure data exceeding 18 MPa.
The methodology to reconstruct the change of pore pressure with time is divided in two steps. During pressure depletion, an exponential function is favoured to capture the pressure decline. Once a waterflooding programme is implemented, RFT data display multiple periods of pressure stabilisation, re-pressurisation, and depletion. Such an irregular pattern combined with the inherent uncertainty associated with each in situ pressure measurement is best fit using a smooth spline function. Furthermore, the calculation of the effective vertical stress is based on the Tergazhi equation that subtracts the pore pressure from the total vertical stress assuming a Biot coefficient equal to 1. The total vertical stress is calculated by assuming a 2.05 g/cm 3 , 2.11 g/cm 3 and 1.05 g/cm 3 density for the unconsolidated and consolidated, shaledominated overburden and seawater, respectively.

Geomechanical model
A strain-rate dependent constitutive model that is a generalisation of an empirical formulation first designed to simulate sandstones (de Waal,  1986) and successfully applied to carbonate and chalk rocks 11,62,63 captures the reservoir deformation. The yield function F that models the yield surface of both Danian and Maastrichtian chalk is expressed as follow: where σ hy is the hydrostatic yield stress at each simulation step and at the current strain rate, P t the tensile strength, p ′ m the mean effective stress, φ the angle of internal friction, m a material parameter determining the shape of the end cap, θ the Lode's angle and is equal to − 30 • , and J 2 the 2nd invariant of the deviatoric stress q with J 2 = q ̅̅ 3 √ for isotropic radial stresses. An associated flow rule is assumed while estimating the viscoplastic strain.
Eq. (1) is a modified Mohr-Coulomb yield function capturing simultaneously the shear failure line determined by the friction angle and tensile strength and the end cap thanks to the hydrostatic yield stress and m parameter. 10 The condition F<0 implies isotropic linear elastic deformation and, when the stress conditions intersect the yield surface and, thus satisfy the condition F = 0, shear or plastic deformation occurs depending on the locus of the intersection point either along the shear failure line or the end cap, respectively. The value of σ hy in Eq (1) is dependent on porosity, water saturation as well as the hardening effect and the deformation rate as follow: where, σ h (Φ, S w ) is the hydrostatic yield stress of chalk at the initial porosity (Φ) and current water saturation (S w ) and at the laboratory vol is the volumetric viscoplastic strain, and ε vp vol the volumetric viscoplastic strain rate. Note that in the strain-rate dependent model, viscoplastic strain represents the contribution of both plastic and creep strain 15 . The term ε geol. is the geological strain rate equal to 2.10 -10 %/hr. The parameter χ is referred to as the hardening coefficient with λ, к and e 0 the two compression coefficients and the initial void ratio, respectively. The hardening coefficient reflects the ability of a rock to sustain an increasing amount of stress beyond pore collapse. Besides, the deformation response of chalk varies according to the strain rate applied to the rock and the softest behaviour is observed at low strain rate. The b-factor dictates the rate-dependent mechanical behaviour of chalk and reflects the creep parameter. 11,22,63 Note that, during unloading caused by, for instance, pressure support from the waterflooding of the reservoir, the modelled uplift of the reservoir top is caused by the elastic relaxation of chalk. For each study window, the 1-D model framework consists of 200 m by side cells stacked vertically from the top of the reservoir to the base of the underburden. A vertical well located below each platform monitored by a GNSS survey (Fig. 2) was selected and the initial density porosity and water saturation were estimated from the well log data. In the Halfdan field, however, a vertical well was missing beneath the monitored platform and an artificial one was created from the 3-D static model of the field built by. 64 The model framework of the studied fields is divided into two main units. Firstly, the reservoir consists of a stratigraphic interval from the Ekofisk Fm. to the base of the M1 unit of the Tor Fm. Using the finest resolution available from the well log data, the cell thickness is 0.15 m for the Kraka and Dan fields, 0.3 m for the Gorm field and 0.6-0.9 m for the Halfdan field. Secondly, the underburden consists of five to eight cells from the M2 to M5 units that are composed of fully water-saturated, low porosity chalk that deforms elastically during simulations. Furthermore, while the base of the model is fixed, uniaxial boundary conditions are assumed such that no strain is allowed to develop on the sides of the model. The effective horizontal stress is assumed isotropic and estimated from the effective vertical stress by using a k 0 ratio (ratio of vertical to lateral earth pressures) equal to 0.7.

Model calibration
The input functions used to calibrate the geomechanical model (Table 1) were determined using data collected exclusively on reservoirs from the Norwegian and Danish sectors of the North Sea such as the Ekofisk, Eldfisk, Tor, Valhall, Hod, Dan, South Arne, and Tyra fields. Firstly, the analysis of the stress-strain curves of fully oil-saturated specimens during triaxial tests at laboratory strain rate (0.1%/hr) allowed estimating the Young's modulus and Poisson ratio 65 as well as the hydrostatic yield stresses (Fig. 3), friction angle, and rock cohesion. 42 The tensile strength was determined by 62 using core data from the Valhall field. Then, the tensile, shear, and plastic properties were inserted in the constitutive model (Eqs (1) and (2)) to reconstruct in p ′ m -q plots the yield surface of chalk. 42 During this process, the value of the m parameter (Eq. (1)), which controls the shape of the end cap was adjusted to fit the experimental data. The water weakening effect was determined by comparing the strength of fully oil-and water-saturated, Maastrichtian specimens measured during 12 hydrostatic compaction tests. 42 Besides, the parameter χ that defines the hardening effect of chalk (Eq. (2); Fig. 1) was estimated while back-analysing seven uniaxial compaction tests carried out on fully oil-and water-saturated specimens of Danian and Maastrichtian age. 65 Finally, transferring the knowledge acquired on chalk geomechanics at the laboratory strain rate to the field depletion rate was carried out by analysing data collected by a formation subsidence monitoring tool (FSMT). FSMT measured periodically along a well in a compacting field from the Danish North Sea the distance between radioactive bullets implanted at known depths. The repetitive measurements allowed estimating the evolution of in situ strain with time. A geomechanical model that included the entire stratigraphic column monitored by FSMT was run during a seven year-long period and the b-factor that dictates the rate-dependent mechanical behaviour of chalk (Eq. (2)) was adjusted for each cell. 65 Once the calibration was completed, the b-factor values were plotted against porosity and the resulting input function was inserted in the constitutive model (Table 1).

Simulation scenarios
Assessing the effect of the yield stress value on the mechanical behaviour of chalk reservoirs is carried out in two steps. Firstly, a conceptual 1-D model is built and composed of four stacked cells simulating the compaction of oil-saturated, Maastrichtian chalk with a porosity of 30%, 35%, 40%, and 45% underlaid by an underburden interval. The rock is subject to a linear increase of the effective vertical stress from 7 MPa to 28 MPa over a 16 yr-long period. Two realisations are run in the Table 1 Porosity-dependent functions defining the geomechanical properties of Danian and Maastrichtian chalk and inserted in Eqs (1) and (2). Φ = porosity, E = Young's modulus, ν = Poisson's ratio, and the ratio К between the hydrostatic yield stress of chalk fully saturated by water and oil represents the water weakening function lowering the yield stress as water saturation (S w ) increases.  (1) and (2) where, ε vp is the normalised difference in vertical viscoplastic strain between the conservative and optimistic approach. ε vp varies between 1 (i.e., the vertical viscoplastic strain ε vp computed in the optimistic model is negligible compared with that of the conservative model) and 0 (i.e., ε vp computed in the optimistic and conservative models are equal). ε vp,cons. and ε vp,opt. are the vertical viscoplastic strain computed in the conservative and optimistic model, respectively. Note that the normalised difference in total strain including both the elastic and viscoplastic strain is also calculated and referred to as ε. Secondly, a conservative and optimistic simulation is run in 1-D along the selected wells for each of the four studied fields (Fig. 2) and the outcomes are compared using ε vp and ε.

Simulation assumptions
Predicting seafloor subsidence based on 1-D simulations requires three key assumptions. A time delay between the onset of reservoir depletion and the onset of seafloor deformation and a subsidence/ compaction ratio (S/C) between the vertical displacement at the seafloor and at the reservoir top need to be assumed. Previous studies on North Sea fields have assessed the time delay and S/C ratio 21,62,67-70 as well as provided methods to estimate them, 71 thereby allowing filling the knowledge gap. Firstly, the present study assumes that seafloor subsidence is initiated when pressure depletion reaches 2 MPa as observed in the Valhall field buried at a depth similar to that of the studied fields.
Secondly, following Geertsma (1973) method that assumes a disc-shaped reservoir and homogeneous properties in the reservoir and the overburden, the S/C in the four studied fields is expected to range between 0.3 and 0.8. For each case study, the value of S/C was adjusted within the 0.3-0.8 interval to fit the seafloor subsidence data and a ratio equal to 0.7, 0.5, 0.4 and 0.3 was assigned to the Dan, Halfdan, Kraka, and Gorm reservoirs, respectively. Thirdly, the onset of change in water saturation during simulation and, thus the start of the water weakening effect is constrained by the well production data. 72 During simulation, the change in water saturation in each cell is assumed to coincide with the onset of water injection with a 1% increase in water saturation per simulation time step (2.4 months). For the Kraka field, no waterflooding recovery scheme was implemented and, thus the decline in the amount of produced oil and associated rise in produced water are used as a proxy to estimate the onset of water saturation increase in the modelled cells.

Conceptual model
As anticipated, changing the yield stress from its initial (conservative model) to final value (optimistic model) shifts the onset of viscoplastic deformation towards high stress conditions. Along the stress-strain and strain-time curves, the transition phase from the elastic to plastic regime is characterised by a progressive acceleration of the accumulated viscoplastic strain until it reaches a steady linear slope in the elasto-plastic phase (Figs. 1, 4a and 4b). After a given amount of compaction, the strain vs stress curve grows exponentially as the result of strain hardening. Considering the cell of the model with a 45% porosity, chalk compaction is first initiated in the conservative model and, thus ε vp value is equal to 1 (Fig. 4c). As the applied stress exceeds σ h0.fin. , irreversible deformation starts occurring in the optimistic model and ε vp slightly decreases to ca. 0.97. The temporary plateau at ca. 0.97 is caused by the incremental increase in viscoplastic strain at the onset of the transition phase of chalk in the optimistic model, whereas in the conservative model strain accumulates at an increasing rate (Fig. 4b and c). The applied stress continues increasing and then the optimistic model experiences the acceleration of the accumulated viscoplastic strain typical of the transition phase, causing an abrupt drop in the ε vp value observed between the year 3 and 5 down to 0.6. After year 5, the difference in strain between both models becomes progressively steady since the geomechanical properties of chalk in the elasto-plastic and strain hardening phases are set to be identical in both the optimistic and conservative approaches. Thus, the contrast in the compaction behaviour of chalk with time in the conservative and optimistic model is characterised by three distinct sections, a high ε vp value close to 1 followed by a drop and stabilisation towards lower values that reflect the shift of the transition phase of chalk between the two simulation approaches. For the cells with a lower porosity, only a section of the reported trend in ε vp with time is observed as chalk becomes stronger and its stress at yielding increases. Note that the temporary plateau observed at ca. 0.97 in the cell with a 45% porosity, is reduced to 0.85 for the 30% porosity cell. This is due to the additional creep strain occurring in the optimistic model during the 9 years period between year 7 and year 16 in Fig. 4b and c.

Pressure depletion
In the block A of the Dan field, pore pressure data cover a period from 1983, 10 yrs after the beginning of production, to 2007 (Fig. 5). After ca. 20 years of production, the pressure depletion stabilises in ca. 1991 at 8 MPa. This depletion, which corresponds to a formation pressure equal to 18 MPa remains relatively constant until 1998. A re-pressurisation phase is then observed reducing depletion by half prior to a recrudescence of  (Fig. 5). After 9 years of production in 1990, the maximum depletion reaches 6 MPa, a value that is maintained relatively constant for 11 years. A re-pressurisation is then observed between 2001 and 2007 at a decreasing rate with time. As in the Dan field, the pore pressure evolution is closely correlated with the implementation of a waterflooding programme. In the Kraka field, no waterflooding programme was implemented and, thus, a continuous depletion is measured from the onset of production in 1991 until 2011. The depletion rate slows down after 6 years of production in 1997 likely due to the decline in the total volume of fluids withdrawn from the reservoir (Fig. 5). Within the Halfdan field, the waterflooding programme was implemented at the beginning of production in 1999. Despite pressure support, pressure depletion occurs continuously to reach ca. 6 MPa. A temporary pressure stabilisation is observed after 6 years of production at a time when water injection reaches its peak.

Seafloor subsidence
The simulation outcomes are compared with the seafloor subsidence data collected on the platforms of each field (Fig. 2). The range of measured seafloor subsidence is between 0.8 m in the Kraka field and 0.07 m in the Halfdan field (Fig. 6). Overall, the GNSS data fall between the seafloor deformation predicted by the conservative and optimistic models, which consider σ h0.in. and σ h0.fin. as the yield stress value in the constitutive equations (Eqs (1) and (2)), respectively. In fact, the observed subsidence trends tend to mainly follow the deformation behaviour captured by the conservative approach. The results also show that, relatively low at the beginning of simulation, the difference in the amount of vertical displacement of the seabed between the conservative and optimistic approach rises significantly as deformation accumulates in the reservoir. The simulated seafloor subsidence using σ h0.in. can be twice that of using σ h0.fin. after a 20 yr-long simulation time as observed in the Dan and Kraka fields. Obviously, this contrast in deformation results from the fact that chalk with a yield stress equal to σ h0.in. is weaker and such compacts more. The additional viscoplastic strain accumulates progressively in the reservoir, thereby promoting the divergence between both models with time.
Furthermore, the successive phases of subsidence and uplift of the seafloor are properly modelled with a vertical displacement rate that is similar to the measured values. In more details, the maximum subsidence rate is observed and simulated in the Dan field and ranges between − 9 cm/yr to − 14 cm/yr, while the rate decreases to − 0.18 cm/yr in the Halfdan field. Positive values indicating uplift are observed in the Gorm and Dan fields reaching temporary up to 3 cm/yr. The positive and negative movements of the seafloor appear closely related to the periods of pore pressure decline and re-pressurisation that are specific to each case study (Fig. 5). The calibrated geomechanical model can thus predict the deformation range of the seabed taking place in different chalk reservoirs characterised by a specific geological setting and exploitation strategy.

Total strain distribution
The deformation in the reservoir is not distributed evenly and strain mainly accumulates within the most porous i.e., the mechanically weakest, chalk intervals. Within the four studied fields, the D1 and M1 units consist of porous rocks compared with the remaining units, the D2 and M2 to M5 units ( Table 2). The highest total vertical strain that exceeds 6% in the Kraka and Gorm fields occurs within these two porous units (D1 and M1). In comparison, deformation in the D2 and M2 to M5 units can be considered negligible with strain values below 0.06%. The most pronounced deformation in these low porosity units is estimated in the conservative model of the Kraka field and is equal to 0.35% total vertical strain. Nevertheless, this strain value represents 5% of the total deformation measured in the Ekofisk Fm. Therefore, the compaction is primary located within the D1 and M1 units, which are further investigated below.
The modelling results show that, within these two latter units, the contrast in the total amount of strain between the conservative and optimistic models and approximated by ε varies among units and fields without any distinguishable trend at a first glance (Table 2). Firstly, in comparison with the optimistic approach, the conservative approach over-estimates the total strain occurring in the D1 and M1 units by a highly variable ε factor spanning from 0.26 to 0.73 and from 0.38 to 0.53, respectively. Thus, the choice of the yield stress value ranging between σ h0.in. and σ h0.fin. has significant impacts on the total computed vertical strain that are not dependent on the stratigraphic formation. Secondly, in the Kraka field, both the conservative and optimistic models give a relatively similar amount of deformation in the Ekofisk Fm. (ε = 0.26), whereas a higher difference (ε = 0.53) is observed in the    Table 2). An opposite observation is reported in the Gorm field. The convergence (ε close to 0) and divergence (ε close to 1) between the outcomes of the conservative and optimistic models are therefore not consistent between the studied reservoirs either. In the following, an assessment of compaction with time is carried out to further characterise the relationships between the assigned yield stress values in Eqs (1) and (2) and the deformation behaviour of the reservoirs.

Initial stress condition
The normalised viscoplastic strain value ε vp (Eq. (3)) at the onset of compaction simulation can be grouped into two categories. In the Dan and Halfdan fields as well as in the D1 unit of the Gorm field, the value of ε vp is high and ranges between 0.8 and 1, denoting a great divergence in the compaction prediction between the conservative and optimistic models (Fig. 7). Note that the absolute strain values at the beginning of simulation are nevertheless low (<0.1%) as shown in Fig. 8. On the other hand, a lower ε vp value spanning between 0.53 and 0.68 and indicating a higher similarity between the models is estimated in the Kraka field and in the M1 unit of the Gorm field. Here the question is what the simulation parameters controlling this contrast between studied fields are.
The initial effective vertical stress is 10.5 MPa for the Kraka and Dan fields, whereas its reaches 12.3 MPa for the Gorm and Halfdan fields.
The ε vp value at the onset of simulation is thus not directly dictated by the virgin stress conditions that can hinder or favour compaction. The stratigraphic intervals associated with a low ε vp value are however systematically characterised by a high mean density porosity. The most porous D1 unit is located in the Kraka field with a mean porosity equal to 38%, whereas the remaining fields display a porosity below 34.9% (Table 2). Similarly, the most porous M1 unit is located in the Gorm and Kraka fields with a mean porosity equal to 34.8% and 32.5%, respectively, whereas that of the Halfdan and Dan fields does not exceed 30.8%. These petrophysical characteristics of the Kraka field and the M1 unit of the Gorm field lead to mechanically weak units that start collapsing at the onset of simulation under virgin stress conditions. As shown in Fig. 8, viscoplastic strain is recorded early in both the conservative and optimistic models of the two latter fields during the first years of simulation, thereby resulting in a low ε vp value. On the other hand, within the Dan and Halfdan fields, and the D1 unit of the Gorm field, chalk is overall tighter and mechanically stronger. The virgin stress triggers viscoplastic deformation first in the conservative model while the mechanical response of chalk from the optimistic model is primarily elastic (Fig. 8). Consequently, the ε vp value is close to 1 at the beginning of simulation (Fig. 7). To conclude, the variability in the computed ε vp values at the onset of simulation between fields and units, or in other words, the degree of divergence in the predicted compaction between the optimistic and conservative models, results from the porositydependent mechanical properties of chalk, which dictate whether the virgin stress induces viscoplastic deformation.

Stress history
The ε vp is monitored throughout the simulation period. Interestingly, during pressure stabilisation i.e., during creep, the ε vp values stay relatively constant independently of the stratigraphic unit and stress history (Fig. 7). The creep deformation does not seem to contribute significantly to the contrast in compaction prediction between the conservative and optimistic models. Furthermore, the ε vp value systematically decreases with time indicating that the total compaction predicted by both models tends to progressively converge (Fig. 7). As observed in the results from the conceptual models (Fig. 4), this convergence results from the increase in stress with time that promotes the viscoplastic deformation of most of the porous cells present in both models and thereby a similar accumulated viscoplastic strain given that except the yield stress value the other mechanical properties (e.g., Young's modulus, hardening and water weakening effect, and time-dependency) are identical. Nevertheless, the magnitude of the decline in ε vp with time varies among fields and units. For instance, the maximum decline from 0.9 to 0.4 is observed in the D1 unit of the Dan field, whereas it remains relatively constant and high in the D1 unit of the Halfdan and Gorm fields with a slight drop from 0.9 to 0.75 (Fig. 7). This contrast between fields in the D1 unit seems related to the specific stress history applied to each reservoir and to the porosity-dependent mechanical properties of chalk as discussed below.
The lowest pressure depletion that is equal to 6.3 MPa is measured within the Gorm and Halfdan fields, whereas formation pressure declines in the order of 8.4 MPa in the Dan field (Fig. 5). Besides, the D1 unit of the Gorm and Halfdan fields is characterised by a low mean porosity (<33.1%) compared with that of the Dan fields (34.9%) ( Table 2). The combination of tight and mechanically strong rock, and low-pressure depletion hinders the viscoplastic deformation of chalk and this is observed in the D1 unit of the Gorm and Halfdan fields. In these two fields, the acceleration of the accumulated viscoplastic strain with time to a steady slope that characterises the transition from the elastic to plastic regime ( Fig. 4a and b) is less pronounced in the optimistic model than in the conservative model (Fig. 8). It indicates that the majority of the collapsing cells remains at the onset of the transition phase, whereas in the conservative model the softer mechanical properties assigned to the Danian chalk promote their deformation until the elasto-plastic phase. As a consequence, the ε vp values remain close to 1 throughout the simulation time (Fig. 7). Contrary to the observations reported in the Halfdan and Gorm fields, the acceleration of the accumulated viscoplastic strain with time in the D1 unit is observed in both the optimistic and conservative models of the Dan field (Fig. 8). As a result, the ε vp value decreases significantly down to 0.4 before the end of the depletion phase (Fig. 7). Within the other intervals such as the M1 unit of all the studied fields and the D1 unit of the Kraka field, the difference in the ε vp value between the onset and end of the simulation runs ranges between the two extreme cases discussed above. The exact value of this difference depends on i) the initial porosity and virgin stress that Table 2 Mean porosity (Φ) and total vertical strain (i.e., elastic and viscoplastic strain) computed by the conservative (ε cons. ) and optimistic (ε opt. ) model and for each unit of the Ekofisk Fm. (D1 and D2 units)  dictate the ε vp value at the beginning of simulation and on ii) the stress history that controls whether most of the reservoir chalks in the optimistic model reach the elasto-plastic phase prior to the end of the pressure depletion.

Discussion
Although scarce in the literature, the simulation studies carrying out sensitivity analyses on the relative influence of the input data and functions (e.g., rock properties, stress, pore fluid composition, reservoir thickness) on the deformation behaviour of hydrocarbon fields can be grouped into two categories. The first category in which the present study belongs consists of running and then comparing multiple compaction realisations of the same reservoir by successively changing one single parameter at a time. Following this methodology, 39 have performed a series of one-way coupled simulations on a carbonate gas field and have demonstrated that the yield stress value is the most influential parameter. Increasing its value by 50% and 100% reduces total seafloor subsidence by 20% and 43%, respectively. In comparison, 23% and 15% reduction in seafloor subsidence was estimated when doubling the Young's modulus value of the overburden and reservoir rock, respectively. Another study focused on land subsidence caused by hydrocarbon production in southern Louisiana has emphasised the impact of the rock compressibility that appears directly proportional to the estimated surface subsidence, whereas the Poisson's ratio has a negligible impact. 73 Besides, a change of 10% of the radius and depth of a disc-shaped reservoir yields an up to 20% variation in the estimated surface subsidence under shallow conditions and, as the depth of the compacting rock increases, this variation is reduced to less than 2%. The present study also indicates the significant influence of the yield stress values on subsurface deformation. The total strain (Table 2) and viscoplastic strain (Figs. 7 and 8) accumulated in the studied fields decline by a factor of minimum 26% when shifting from the initial to final yield stress value that represents a maximum yield stress increases of 30% and 40% for the Maastrichtian and Danian chalk, respectively (Table 1). However, this study also demonstrates that the influence of the yield stress values on compaction is not straightforward and its magnitude dependent on several factors including the rock types (porosity and age of the rock), the virgin stress, and the stress history. These three latter parameters dictate the time at which pore collapse is initiated and the amount of accumulated strain in each realisation. Such quantitative assessment of the main parameters controlling the compaction behaviour of producing fields represents the first step towards accurate uncertainty analysis of the deformation predicted by geomechanical models. 74 The first category of sensitivity studies is however time-consuming as numerous simulations are required considering that the constitutive model contains many input variables such as the elastic modulus, yield stress, in situ stress path, hardening coefficient, b-factor, and the shape of the yield surface (Eqs (1) and (2)). An alternative is to utilize genetic algorithms, machine learning methods, and other optimising methods to screen the most influential geomechanical properties and quantify their impacts in the context of subsurface deformation. This second category of sensitivity studies estimates the rock properties by optimising the error between simulation outcomes and observations such as, seafloor subsidence and in situ strain measurements. 36,37,40,75,76 Mainly focused on land subsidence related to groundwater and hydrocarbon production, these studies have also highlighted the role of the plastic properties of the rock such as the hardening coefficient as well as the in situ stress path, and to a lesser extent, the elastic modulus and Poisson's ratio on field deformation. Although fast to execute, a too simplistic view of the reservoir is commonly considered for simulations with no lateral heterogeneity in the rock properties. This is due to the deterioration of the optimisation process by adding geological and thus, geomechanical heterogeneity into the models. 38 Importantly, inverse modelling carried by 69 on a Danish North Sea field have demonstrated that a simplistic model composed of one single cell and a more realistic model composed of several cells can both provide a good fit to the data (i.e. seafloor subsidence) given that appropriated properties are assigned to each cell. Thus, a model capable of reproducing deformation observed at the field scale is not a guarantee that the optimised geomechanical properties F. Amour et al. assigned to the cells reflect the true characteristics of the subsurface rock. Both categories of sensitivity studies are therefore complementary in assessing the uncertainty of the outcomes of geomechanical models related to the input variables. While the first category of studies accurately quantifies the influence of each input parameter within a realistic geological framework, the second one rapidly screens the key variables impacting field deformation.

Conclusion
The present study assesses how the uncertainty associated with the determination of the yield stress of chalk samples in the laboratory manifests itself in a compacting field during 1-D simulation. For each of the four Danish North Sea fields investigated, two compaction simulations are carried out by integrating alternatively in the constitutive equations the σ h0.in. (conservative approach) and σ h0.fin. (optimistic approach) values of Danian and Maastrichtian chalk that delimit the stress interval between the elastic and plastic regime. The results show that the seafloor subsidence measured in each field falls between that predicted by the conservative and optimistic approaches, thereby validating the geomechanical model under different geological settings and production schemes. In addition, the differences in the viscoplastic behaviour of the four hydrocarbon fields between the conservative and optimistic approach are controlled by three parameters, i) the age and porosity of chalk, ii) the virgin stress condition, and iii) the stress history. According to these results representative of hydrocarbon fields located in the southern part of the North Sea, the viscoplastic strain and seafloor subsidence computed in the reservoirs vary between both approaches by a factor of 26%-73% and 30%-50%, respectively. The simulation outcomes also indicate that the creep deformation does not seem to contribute significantly to the contrasts in compaction prediction between the conservative and optimistic models.
Furthermore, the porosity distribution in each studied area as well as the in situ stress applied to the rock dictate the timing and the amount of viscoplastic strain occurring in each realisation. A high pressure depletion favours the collapse of most the porous cells and, thus a convergence of the simulation outcomes between the conservative and optimistic models. However, under low effective vertical stress, the cells with the mechanically weakest chalk in the conservative model enter in the plastic regime, whereas the same cells in the optimistic realisation that are mechanically stronger, are still in the elastic domain or in the transition phase between the elastic and plastic regime. A great discrepancy in the predicted reservoir compaction between the conservative and optimistic models is observed in the latter conditions. Sensitivity studies during geomechanical modelling are of primary importance to document the different compaction scenarios that a field can experience. These compaction scenarios can then be used to create probability maps, which identify the most mechanically unstable areas of a field that need to be avoided while drilling and installing infrastructures in the context of hydrocarbon production as well as CO 2 and thermal storage.

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.

Data availability
Data will be made available on request.