Analysis of Caprock Tightness for CO₂ Enhanced Oil Recovery and Sequestration: Case Study of a Depleted Oil and Gas Reservoir in Dolomite, Poland

Abstract

This study addresses the problem of geological structure tightness for the purposes of enhanced oil recovery with CO₂ sequestration. For the first time in the history of Polish geological survey the advanced methods, practical assumptions, and quantitative results of detailed simulations were applied to study the geological structure of a domestic oil reservoir as a potential candidate for a combined enhanced oil recovery and CO₂ sequestration project. An analysis of the structure sequestration capacity and its tightness was performed using numerical methods that combined geomechanical and reservoir fluid flow modelling with a standard two-way coupling procedure. By applying the correlation between the geomechanical state and transport properties of the caprock, threshold pressure variations were determined to be a key factor affecting the sealing properties of the reservoir–caprock boundary. In addition to the estimation of the sequestration capacity of the structure, the process of CO₂ leakage from the reservoir to the caprock was simulated for scenarios exceeding the threshold pressure limit of the reservoir–caprock boundary. The long-term simulations resulted in a comprehensive assessment of the total amount of CO₂ leakage as a function of time and the leaked CO₂ distribution within the caprock.

1. Introduction

Because CO₂ emissions are an increasing problem, many strategies have been developed to reduce its emissions to the atmosphere, including mineralisation [1,2], CO₂ storage at the bottom of the oceans [3], accelerated weathering [4,5], and subsurface geological storage [6]. Carbon dioxide capture at the source, followed by its long-term storage in exploited hydrocarbon reservoirs, is one of the most practical ways to reduce the CO₂ concentration in the atmosphere.

Reservoir structures suitable for geological CO₂ storage are most frequently water-bearing formations, depleted hydrocarbon deposits, or unminable methane-rich coal seams.

Geological formations that are selected for CO₂ storage must meet the appropriate criteria related to the minimum and maximum depth of the structure. The reservoir rock parameters, including thickness, permeability, porosity, and fracture characteristics, are also critical. However, the most important criterion is the sealing quality of the rock, that is, the overburden’s integrity and thickness play key roles in assessing the safety of long-term CO₂ geological storage.

Most of these criteria are met for depleted oil and gas fields because hydrocarbon production from oil and gas fields implies favourable reservoir properties, including porosity, permeability, and reservoir thickness [7].

During the geological storage of CO₂, most of the injected CO₂ remains in the free (liquid, gas, or supercritical) phase. This phase poses the greatest threat to the seal integrity. Under reservoir conditions, which are at elevated pressures and high temperatures, CO₂ occurs in a dense gas phase, yet its density is lower than that of the formation water. Driven by the buoyant force that occurs as a result of the difference in the densities of the CO₂ and reservoir water, the free fraction of CO₂ tends to migrate upward [8,9,10]. The presence of free CO₂ in the upper parts of reservoir rocks poses a threat of escape to overlying formations and potentially into the atmosphere, in the case of the caprock integrity loss.

Therefore, it is crucial to maintain a tight sealing layer and not inject CO₂ into the formation beyond a level that guarantees the integrity and stability of the seal and overall safety of the operation. The long-term tightness of the overburden may be influenced by the effects of CO₂ injection, as operations change the stress and strain field in subsequent years of sequestration.

The aim of this study was to analyse the tightness of the caprock above the main dolomite (Ca2). The main dolomite structure with partially depleted hydrocarbon accumulation located in the Gorzów Block, Poland, is considered a potential formation for carbon dioxide storage.

Numerous research papers address the caprock tightness problem in the potential site for CO₂ geological storage. Edlmann et al. (2013) [11] postulated the existence of a critical threshold of fracture aperture size controlling the CO₂ flow along the shale samples. Li et al. (2006) [12] examined the occurrence of the volume flow and measured the effective gas permeability for selected post-failure evaporite beds samples caused by CO₂ injection. Zivar et al. (2019) [13] examined the effect of stress magnitude and stress history on porosity and permeability values of anhydride and carbonate rocks, while Hangx et al. (2009) [14] investigated the impact of CO₂ on the mechanical strength of Zechstein Anhydrite, which seals many potential CO₂ storage sites in Central and Eastern Europe.

Under the assumptions of the study, CO₂ was injected into the reservoir rock for both enhanced oil recovery (EOR) and geological CO₂ sequestration. To meet the research objectives, several numerical methods were used, involving the coupling of geomechanic and dynamic modelling of reservoir fluid flow in the main dolomite formation, into which CO₂ was injected, and the surrounding rocks, especially the overburden.

Geological Setting

The study area is located on the Gorzów Block adjacent to Szczecin Trough in the north, with the Mid-Polish Swell to the east and the Fore-Sudetic Monocline to the south (Figure 1). This area has a regional sequence of tectonic disturbances and related uplift of Permian-Mesozoic sediments [15]. These elevated tectonic blocks were accompanied by extensive volcanic rock cover, as well as a depressions series of clastic deposits in the Lower Rotliegend. In particular, the erosive outliers of the Zechstein basement had a significant impact on the development of the overlying Zechstein--Mesozoic complex [16].

In the analysed part of the South Permian Basin, the Zechstein Sea entered a morphologically diverse depositional surface. Considerable denivelations contributed to the division of the basin into shallow and deep-water zones. In the elevated areas, platforms and micro platforms of sulphate deposits of the Werra cyclothem developed and became covered with the main dolomite platforms, which formed as a result of the transgressive cycle during the Stassfurt cyclothem [16]. Sedimentological studies of the main dolomite deposits revealed the existence of different sedimentation environments resulting from considerable bathymetric differences [17,18,19]. These include: the platform and microplatform zones (including barrier and platform flat), slope and toe of the slope, and basin floor [19].

Optimal reservoir properties were found both in the shallow-water platform formations and in the formations developed at the toe of the slope in a deep-water environment [18,19,20].

From a petroleum exploration and potential storage volume perspective, the most important property was the development and preservation of significant secondary porosity in the main dolomite rock (in some areas over 30%), which was caused by the complete or partial dissolution of granular carbonate components. The presence of this high porosity may be related to the dolomitisation and recrystallisation of the primary calcareous rock matrix [20,21].

The main dolomite, which is the reservoir rock and target formation for carbon dioxide storage, is covered with a sequence of thick evaporates, thin carbonates, and very thin layers of shales of the PZ2 (Stassfurt), PZ3 (Leine), and PZ4 (Aller) Zechstein cycles [22,23]. The location of the study area with the lithostratigraphic profile is shown in Figure 1.

Figure 1. Lithostratigraphic profile in the reference borehole (left) and the location of the study area on the background of the map of the main dolomite (Ca2) in Poland (right) [24].

Figure 1. Lithostratigraphic profile in the reference borehole (left) and the location of the study area on the background of the map of the main dolomite (Ca2) in Poland (right) [24].

The thickness of the Zechstein evaporates overlying the reservoir rock varies from 253 m on top of the platform flat and reaches 840 m in the basin. Developed over geologic time on top of the reservoir rock sequence, these hardly permeable or impermeable evaporates constitute the seal for hydrocarbon accumulation and are potentially a good barrier for carbon dioxide stored in the main dolomite reservoir rock.

2. Methods

Performing fluid flow modelling coupled with geomechanics is essential for the safe storage of CO₂ or other media in subsurface rock formations. In particular, analysing the mechanical response of rocks to activities related to geological CO₂ storage concerns the reservoir rock itself and the neighbouring rocks, especially the overburden.

As a result of injection, the increased pressure of the fluid filling the pore space in the rock causes a change in the stress field. The disturbance of the stress equilibrium may lead to a potential risk of loss of tightness of the sealing rocks, creating a potential pathway for CO₂ leakage to the overlying strata. To evaluate the tightness of the caprock, we employed numerical methods that combined geomechanic and reservoir fluid modelling of the reservoir and surrounding rocks using the Schlumberger software (Houston, TS, USA)—Petrel™ platform, VISAGE™ geomechanical simulator, and ECLIPSE™ reservoir simulator. The standard approach to this evaluation includes a two-way coupling of both types of simulations realised by iterative, alternating runs of the simulations, as shown schematically in Figure 2.

For the geomechanical analysis, parameter distributions describing the geomechanical and petrophysical rock properties were developed, and boundary conditions were determined. The distributions of the pore pressure necessary for the geomechanical simulations were combined with the hydrostatic pressure distribution in the overburden and the reservoir pressure in the main dolomite (Ca2) developed for certain time steps in the history of hydrocarbon production and CO₂ injection. Upgraded transport properties were determined according to the resultant geomechanical state of the structure and were used in the dynamical model to simulate upgraded pressure and saturation distributions. These results were input to the geomechanical model and thus closed the iteration loop of the two-way coupling simulation procedure.

Relevant Dataset

To analyse the tightness of the geological structure’s caprock, we used three main data types:

• Seismic data, which comprises the result of a 3D seismic interpretation of the study area, the structural surface of the top of the reservoir rock (Ca2) in the depth domain, and the map of the thickness of the reservoir rock based on seismic data;
• Well log and lab data;
• Lithostratigraphic profiles and well log data from 10 of the 27 boreholes drilled in the study area with the results of laboratory measurements of petrophysical and static geomechanical parameters performed on the core material;
• Reservoir engineering data (reservoir fluid saturation distribution, pressure distribution, and reservoir fluid thermodynamic (PVT) properties);
• Hydrodynamic well tests (multi-rate and pressure build-up tests);
• Production data (reservoir fluid production rates and totals, and bottom-hole and well-head pressures).

3. Geological Model

To meet the study’s objectives, we first developed a structural model of the area to define the 3D space, which was subsequently parameterised. Ultimately, we obtained the spatial distributions of the petrophysical parameters necessary to determine the weight of the overburden rocks via geomechanic modelling and during the modelling of the reservoir fluid flow.

3.1. 3D Structural Geological Model

The structural 3D model of the study area presents the geometry of the entire study area profile, specifically the reservoir rock (Ca2), but also, for the purposes of geomechanic modelling, the structure of the overburden rocks reaching up to the ground surface, as well as the side- and under-burden. The developed spatial structural model was used to determine the 3D space for parameterisation.

Using the structural map of the main dolomite Ca2 top, the map of the Ca2 thickness, and the stratigraphic borehole data, we developed a map of the bottom of the main dolomite.

Further, using available interpretations of seismic horizons in the depth domain and boreholes stratigraphic markers, as well as tools available in the Petrel software, we constructed structural maps of the remaining stratigraphic levels. These maps began with the Zechstein limestone (Ca1), through the lithostratigraphic units of the Werra (A1D, Na1, and A1G), Stassfurt (Ca2, A2, Na2, and A2G), Leine (I3 and A3), and Aller cyclothem sequence (A1D, Na1, and A1G) occurring in the area of the study. They are followed by the sediments of the Triassic, Jurassic, Cretaceous, and finally, the Paleogene, Neogene, and Quaternary sediments. The longitudinal extent of the structural model is approximately 13.8 km, and the latitudinal length is approximately 14.5 km.

The horizontal resolution of the 3D grid was 100 m × 100 m. The vertical resolution of the model varied depending on the individual zones. To assess the overburden’s tightness, it was necessary to detail the zones that could play a key role in the possible leakage of the stored CO₂, which included the rocks located in the vicinity of the main dolomite with minimal permeability. Accordingly, the highest resolution was used in the zones critical for further investigations, such as the reservoir rock (the main dolomite Ca2) and its closest sealing layer, the basal anhydrite (A2). The mean vertical resolutions in these intervals were 1.88 m and 4.85 m, respectively.

The final geometry of the 3D grid used in the geomechanical simulation of the VISAGE™ simulator (Schlumberger) considered the rocks surrounding the reservoir rock. The geomechanical model was limited from the top at the ground surface, and at the bottom, which is defined at a depth of approximately 41 km. At the sides of the reservoir scale model, a volume of rocks with a length of approximately 45 km on each side was added, resulting in the final latitudinal extent of approximately 103.8 km and a longitudinal extension of 104 km.

As shown in Figure 3, the geometry of the obtained spatial structural model of the main dolomite reservoir rock and surrounding rocks, including the overburden, were identified.

This developed structural model was then parameterised, resulting in spatial models of the petrophysical and geomechanical parameters.

3.2. 3D Modelling of Petrophysical Properties

Because the modelling workflow included a geomechanic modelling task that required the definition of rock properties for the entire geological profile, the structural model had to be further developed. Property modelling was carried out in two different model scales: (1) a model covering the geological profile up to the ground level, which was populated with total porosity and rock density values; and (2) a detailed 3D model of the interval of interest covering the Ca2 main dolomite reservoir rock, in which a higher vertical resolution was applied.

The interval within the geological profile, which is of main concern in this study, was modelled with increased vertical resolution to enable a more detailed modelling and interpretation of the geomechanical and fluid flow processes. Figure 4 shows the division of the stratigraphic and lithological units into layers. The reservoir interval of the Ca2 was divided into 6 layers, and the neighbouring sediments above and below were divided into varying numbers of layers depending on their average thickness in order to achieve a vertical resolution of approximately 5 m within the Ca2 and approximately 25 m in the salt rock formations, with anhydrite intervals of several meters thick.

3.2.1. Density and Porosity Models of Entire Geological Profile

The property modelling was based on the following data sources: geophysical wellbore logging and interpretation along the entire wells’ profiles, the results of the laboratory measurements conducted for the main dolomite reservoir rock, which were used to constrain petrophysical interpretation, and the 3D seismic data (previously transformed into seismic attributes form), which were applied as a secondary input in the 3D property modelling of the Ca2 reservoir porosity.

For the global grid population with porosity and density values, eight borehole profiles were used. The analysis included the definition of density and porosity variation ranges within each stratigraphic unit, as well as the modelling of semivariograms. Further, 3D density and porosity modelling were accomplished using a stochastic algorithm (Gaussian random function simulation), in which the modelling process was iterated 20 times to produce 20 equally probable realisations of the porosity and density 3D models for the entire profile. Finally, the arithmetic averages of the 20 realisations were calculated for geomechanic modelling.

3.2.2. High-Resolution 3D Petrophysical Model of Target Reservoir

The static modelling workflow was focused on the detailed modelling of the petrophysical properties of the Ca2 main dolomite, which is a potential CO₂ storage reservoir.

For this reason, extended datasets were used, including petrophysical interpretation along eight boreholes calibrated with dense laboratory data, and 3D seismic volumes of simultaneous inversion and seismic attributes.

Porosity reflects the volume of pore space within the reservoir, which, in this case, was the oil-saturated dolomite interval. Meanwhile, permeability defines the reservoir’s capability to conduct fluid flow (oil, natural gas, water, and CO₂). Thus, the 3D models of porosity and permeability are of the highest importance for dynamic modelling the fluid flow within a storage reservoir.

The porosity model was elaborated using both primary wellbore data and secondary seismic attributes previously transformed into seismic properties, which exhibit higher correlation coefficients with porosity along the wellbore profiles. The Gaussian random function simulation algorithm was applied in a co-kriging version. To develop the permeability model, porosity vs. permeability relationship was combined with the borehole profiles of permeability. Figure 5A shows the outcome of porosity modelling within the Ca2 reservoir, while Figure 5B shows the permeability 3D distribution.

4. Geomechanic Model

During hydrocarbon production, the reservoir pressure declines, and the effective stresses increase. Conversely, during CO₂ injection, to utilize the reservoir pressure as a tertiary oil recovery method (EOR) or for CO₂ geological storage, there is an increase in pore pressure, which causes a decrease in the effective stresses. This effective stress principle was first proposed by Therzagi (1945) [25], which was later applied to rocks to understand their ultimate strength and ductility [26].

The change in stress within the reservoir and surrounding rocks can be expressed with an equation for isotropic rocks as follows [27]:

$${\mathsf{\sigma }}_{\mathrm{h}}=\frac{\mathsf{\nu }}{1-\mathsf{\nu } }\left({\mathsf{\sigma }}_{\mathrm{v}}-\mathsf{\alpha }\mathrm{P}\right)+\left[\frac{\mathrm{E}}{1-{\mathsf{\nu }}^2 }\right]{\mathsf{\epsilon }}_{\mathrm{h}}+\left[\frac{\mathrm{E}\mathsf{\nu }}{1-{\mathsf{\nu }}^2 }\right]{\mathsf{\epsilon }}_{\mathrm{H}}$$

(1)

where:

–

${\mathsf{\sigma }}_{\mathrm{h}}$ and ${\mathsf{\sigma }}_{\mathrm{v}}$ are the minimum horizontal and vertical stresses, respectively;

–

$\mathsf{\alpha }$ is the Biot’s coefficient;

–

$\mathrm{P}$ is the pore pressure;

–

$\mathsf{\nu }$ is the Poisson’s ratio;

–

$\mathrm{E}$ is the Young’s modulus;

–

${\mathsf{\sigma }}_{\mathrm{h}}$ and ${\mathsf{\epsilon }}_{\mathrm{H}}$ are the strains in the direction of the minimum and maximum horizontal stresses, respectively.

In geomechanic modelling, a set of 3D models of geomechanical parameters were developed, and then, the boundary conditions were determined to calculate the stress field and predict the geomechanical behaviour of the rock during the period of hydrocarbon production and CO₂ injection into the reservoir rock.

4.1. Modelling of Geomechanical Properties

To model the distribution of geomechanical properties in 3D space, the Young’s modulus, Poisson’s ratio, unconfined compressive strength (UCS), and tensile strength were modelled in 1D for 10 boreholes before they were modelled in the geological model 3D space.

4.1.1. 1D Modelling of Elastic and Strength Properties

Herein, the models of the static elastic properties of the rocks in 1D, including Young’s modulus and Poisson’s ratio, were conducted using the well logs from 10 boreholes, the results of laboratory measurements of the static Young’s modulus parameter, and the Poisson’s ratio obtained from the core material from these boreholes. Based on the well log data, the dynamic Young’s modulus and Poisson’s ratio were calculated, wherein the relationships between the borehole velocity of the compressional vp, shear wave vs, and density ρ were used as follows [28,29,30]:

$${\mathrm{v}}_{\mathrm{dyn}}={\mathrm{v}}_{\mathrm{p}}{}^2-\frac{{\mathrm{v}}_{\mathrm{s}}{}^2}{2}\left({\mathrm{v}}_{\mathrm{p}}{}^2-{\mathrm{v}}_{\mathrm{s}}{}^2\right),$$

(2)

$${\mathrm{E}}_{\mathrm{dyn}}={\mathsf{\rho } \mathrm{v}}_{\mathrm{s}}{}^2\left[\left(3{\mathrm{v}}_{\mathrm{p}}{}^2-4{\mathrm{v}}_{\mathrm{s}}{}^2\right)\left({\mathrm{v}}_{\mathrm{p}}{}^2-{\mathrm{v}}_{\mathrm{s}}{}^2\right)\right],$$

(3)

where ${\mathrm{v}}_{\mathrm{dyn}}$ denotes the dynamic Poisson’s ratio, ${\mathrm{E}}_{\mathrm{dyn}}$ represents the dynamic Young’s modulus, ${\mathrm{v}}_{\mathrm{p}}$ is the velocity of the compressional wave, ${\mathrm{v}}_{\mathrm{s}}$ is the velocity of the shear wave, and $\mathsf{\rho }$ is the rock density.

The static equivalents of the elastic properties were determined by comparing the calculated dynamic properties with the results of laboratory measurements of the static elastic properties. The correlation coefficients were 0.640 and −0.588 for the Young’s modulus and Poisson’s ratio, respectively.

To estimate the static uniaxial compressive strength in the profiles of the analysed boreholes, the compressional wave velocity was compared with the results of uniaxial compressive strength laboratory measurements, which produced a linear relationship with a satisfactory correlation coefficient of 0.682.

Further, to estimate the static tensile strength, T, which was not measured on the core material, we used a simple and well-known relationship proposed by Hoek (1966) [31]:

$$\mathrm{T}=\frac{\mathrm{UCS}}{10}.$$

(4)

As shown in Figure 6, the input data and developed 1D models of static elastic parameters, including the Poisson’s ratio, Young’s modulus, and $\mathrm{UCS}$ (marked with continuous line), were calibrated with the results of their measured static equivalents (marked with dots).

4.1.2. 3D Modelling of Elastic and Strength Properties

The profiles of the static geomechanical parameters from the study area were scaled up to obtain their average values along the boreholes with a defined vertical resolution of the structural model. Then, they were subjected to a geostatistical analysis, which included determining the degree of spatial correlation among the borehole data for each geomechanical parameter.

While calculating the spatial distribution of the geomechanical parameters, the truncated Gaussian simulation algorithm (TGS) with the Co-kriging option was employed, thereby allowing us to use the 3D seismic data as secondary data. The advantage of this approach is that it compensates for missing data and, therefore, unknown relationships, especially in the horizontal direction. As a result, we obtained 3D models of the geomechanical parameters, including the Young’s modulus, Poisson’s ratio, and UCS, in the main dolomite reservoir interval, as shown in Figure 7A–C, respectively. Moreover, in Figure 7D–F, the average values of the static Young’s modulus, Poisson’s ratio, and UCS in the reservoir zone, respectively, are shown.

The spatial distributions of these geomechanical properties exhibited typical relationships, in which there was a negative correlation between the elastic parameters and a positive correlation between the Young’s modulus and UCS. These relationships were also observed in the distribution of the values of these parameters illustrated as histograms. The 3D distribution of the tensile strength was obtained in a similar manner as the 1D models, wherein Hoek’s relationship with UCS was utilized. Other geomechanical parameters, such as friction angle, angle of dilatation, and the Biot’s coefficient, were assigned based on the typical values for appropriate lithologies listed in the literature [32,33,34,35,36]. The properties used as an input in the geomechanic modelling, both as a result of the modelling procedure and assumed based on the literature, are listed in

Table 1. Calculated and estimated values of petrophysical and geomechanical parameters in 3D model of the study area.

Parameter [Unit]	Cenozoic (Clay, Sand, Gravel)	Cretaceous (Clayey Shales)	Jurassic (Sandy Shales)	Triassic (Sandstones)	Zechstein
					Rock Salt	Anhydrite	Reservoir MAIN Dolomite	Limestone	Rotliegend (Underburden)
Young’s modulus [GPa]	0.5	4	5.56	28.5	6.89	52.69	3D model	42.06	46.19
Poisson’s ratio [-]	0.3	0.32	0.19	0.17	0.3	0.25	3D model	0.18	0.3
Density [g/cm3]	3D model	3D model	3D model	3D model	3D model	3D model	3D model	2.75	2.3
Biot’s coefficient [-]	1	1	1	1	0	0.10	0.7	0.8	1
Porosity [%]	3D model	3D model	3D model	3D model	3D model	3D model	3D model	2.99	4
Unconfined compressive strength (UCS) [MPa]	2.8	48	56.98	50.7	27.33	90.3	3D model	14.93	50
Friction angle [°]	30	32	20	59	29.08	64	28.6	22.8	30
Dilatation angle [°]	0	0	0	0	0	0	0	0	0

.

4.2. Boundary Conditions

The boundary conditions applied to the model in order to define the initial stress for the geomechanical simulation were determined as the global tectonic stresses (minimum horizontal stress, ${\mathsf{\sigma }}_{\mathrm{h}}$, and maximum horizontal stress, ${\mathsf{\sigma }}_{\mathrm{H}}$) based on a recent investigation of a tectonic stress field in Poland conducted by Jarosiński (2006) [37]. The characteristics of the assumed stress field are listed in

Table 2. Characteristics of regional principal horizontal stresses in the neighbouring area based on borehole breakouts (Jarosiński, 2006).

Stress Characteristic Parameter	Assigned Value
Gradient of horizontal stress (${\mathsf{\sigma }}_{\mathrm{h}}$) [MPa/m]	0.01707
Gradient of ${\mathsf{\sigma }}_{\mathrm{H}}$ [MPa/m]	0.02134
Azimuth of ${\mathsf{\sigma }}_{\mathrm{H}}$ [°]	6

. Note that vertical stress, ${\mathsf{\sigma }}_{\mathrm{v}}$, is caused by gravity, and was thus determined using the overburden density (3D distribution of rock density).

The developed 3D models of the geomechanical and petrophysical parameters were used as inputs to the geomechanical simulation. These were run for seven time steps (2003, 2020, 2050, 2052, 2053, 2056, and 2100) over the course of the hydrocarbon field exploitation and CO₂ injection, which were defined by the distribution of the reservoir fluid pressure developed in the reservoir fluid flow modelling. The calculated stress field reflecting the geomechanical states of the analysed rocks at specific time steps were then used to evaluate the tightness of the caprock, which could potentially leak the injected and stored CO₂.

5. Dynamic Model

A dynamic model of the analysed structure was constructed using the developed geological model and was supplemented with the following additional components:

–

An initial distribution of reservoir fluids (oil and water) under hydrostatic conditions;

–

Reservoir fluid transport properties (relative permeabilities);

–

Reservoir fluid (oil) thermodynamic model.

5.1. Reservoir Fluid Distributions

The presence of three reservoir fluids (water, gas, and oil) was assumed in the simulation model of the analysed reservoir. The water–oil contact was established at a depth of 3280 m b.s.l., while the gas–oil contact was established at a depth of 3178 m b.s.l. The J-Leverett function was used to generate the initial water saturation distributions that matched those obtained from the geophysical measurements. This approach enabled us to obtain the initial dynamic equilibrium by simultaneously reconstructing the water saturation values in the wells. The J-Leverett function of water saturation, S_w, considers the dependence of capillary pressure, ${\mathrm{P}}_{\mathrm{cow}}$, on the parameters of the reservoir rock and has the form:

$$\mathrm{J}\left({\mathrm{S}}_{\mathrm{w}}\right)=\sqrt{\frac{\mathrm{k}}{\mathsf{\varphi }} } \frac{{\mathrm{P}}_{\mathrm{cow}}\left({\mathrm{S}}_{\mathrm{w}}\right)}{{\mathsf{\sigma }}_{\mathrm{ow}}\cos \left({\mathsf{\theta }}_{\mathrm{ow}}\right)},$$

(5)

where $\mathrm{k}$ and $\mathsf{\varphi }$ denote the permeability and porosity of the reservoir rock, respectively, and ${\mathsf{\sigma }}_{\mathrm{ow}}$ and ${\mathsf{\theta }}_{\mathrm{ow}}$ are the interfacial tension at the oil–water interface and their contact angle, respectively. The correct fit of the water saturation distributions in the wells was obtained using the four-parameter J curve model as follows:

$$\mathrm{J}=\left(\frac{\mathrm{A}}{{\mathrm{S}}_{\mathrm{r}}^{\mathrm{n}}+\mathrm{B}}\right)\left(1-{\mathrm{S}}_{\mathrm{r}}^{\mathrm{p}}\right), {\mathrm{S}}_{\mathrm{r}}=\frac{{\mathrm{S}}_{\mathrm{w}}-{\mathrm{S}}_{\mathrm{wc}}}{1-{\mathrm{S}}_{\mathrm{wc}}},$$

(6)

where ${\mathrm{S}}_{\mathrm{wc}}$ is the connate water saturation, in which the correlation parameters were found to be: $\mathrm{A}$ = 0.0094, $\mathrm{B}$ = 0.0657, $\mathrm{n}$= 2, $\mathrm{p}$ = 10, and ${\mathrm{S}}_{\mathrm{wc}}$ = 0.0035, while ${\mathsf{\sigma }}_{\mathrm{ow}}=64.65$ dyne/cm and ${\mathsf{\theta }}_{\mathrm{ow}}=0°$.

The water saturation distributions were adjusted for every well. As a result of this analysis, it was found that the model correctly reproduced the water saturation variation with depth in most wells. As a result, the average water saturation in the model generated corresponds with high accuracy to the average water saturation generated using the interpolation method in the geological model.

5.2. Transport Properties

In general, transport properties in porous media include various flow mechanisms ranging from continuous flow of viscous type (Darcy flow) through slippage effects to Knudsen and molecular diffusion. Similarly, thermodynamics associated with transport phenomena may obey classical laws of locally equilibration state for large porosity systems and follow non-equilibrium type for small porosity system where molecule collisions with the walls of the system dominates over inter-particle ones. In contrast to low permeability, low porosity rocks [38,39] (tight and shale formations), conventional rocks such as dolomites included in the analysed case reveal standard viscous flow characterized by permeability tensors and equilibrium thermodynamic variables like pressure drop and shear stress. The other type of rocks in the analysed model refers to caprock anhydrite. Although it is an extremally low permeability rock, its main feature used in this study concerns the dependence of transport properties upon geomechanical state that was established in terms of effective permeability [13]. Consequently, in order to model effects of gas transport across the anhydrite caprock, a conventional approach was applied adopting stress-dependent permeability of a single porosity micro-fractured system.

An important part of the transport described by permeabilities is played by relative permeabilities of various reservoir fluids. Owing to the lack of direct measurements of the relative permeability curves, k_r, and reduced fluid saturation, S_r, for the analysed reservoir, a power model was adopted according to the relationship: k_r = (S_r)ⁿ, wherein reduced saturation was determined as follows:

$${\mathrm{S}}_{\mathrm{r}}=\frac{\mathrm{S}-{\mathrm{S}}_{\min }}{{\mathrm{S}}_{\max }-{\mathrm{S}}_{\min }},$$

(7)

where S_min and S_max refer to the minimum and maximum available saturations, respectively. The model used the following parameter values: for water (k_r = k_rw, S = S_w): irreducible water saturation S_min = S_wcr = 0.0528, maximum water saturation, S_max = 1; for oil in the oil–water system (with no gas k_r = k_row, S = S_o): S_min = S_orw = 0.4917, S_max = S_omax = 1 − S_wco = 0.9964; for oil in the oil–gas system (at connate water saturation, k_r = k_rog, S = S_o) S_min = S_orw + S_wco, S_max = S_omax; and for gas in the oil–gas/water–gas system (k_r = k_rg, S = S_g) S_min = S_gcr = 0.1, S_max = S_gmax = 1 − S_wco = 0.9964. Both the exponent of the relationship and the other parameters were adjusted during the model calibration procedure.

5.3. Reservoir Fluid Model

As the simulation model of the analysed reservoir (structure) is compositional, a thermodynamic model of the reservoir fluid (oil and gas) was constructed and calibrated based on the Soave–Redlich–Kwong (SRK) equation of state (EoS) using a standard software tool (PVTSim software, [40]). SRK equation of state is a conventional equation that relates temperature, pressure and volume of fluid and uses two parameters taking into account interaction of fluid molecules and their volumes. Those parameters are expressed in terms of several physical quantities, such as critical pressure, temperature, and volume. Other quantities include eccentricity factors [41], molar masses and boiling point temperatures. All these quantities are well defined for single component fluids such as light hydrocarbons (up to C6) and non-hydrocarbon components (N₂, CO₂, H₂S). For heavier hydrocarbons they are determined by empirical correlations with densities and molar masses [42]. In order to effectively perform compositional reservoir simulations, the number of fluid components is typically reduced by grouping or lumping the heavier ones. Effective parameters of such pseudo-components are calculated as averages of the real components with weights proportional to their concentrations and molar masses. In addition, the equation of state for multicomponent mixtures includes interaction effects of different polar molecules by so called mixing rules expressed via binary coefficients [43]. The viscosities of gas and liquid phases of the analysed fluid were determined using Lohrenz–Bray–Clark correlation [44] of viscosity and density and the results of the SRK equation of state for the densities. The SRK equation of state with parameters described above is typically sufficient to predict fluid properties. However, its results can be improved by the method of EoS regression [45] to experimental data.

In the analysed case the original fluid model [46] included 25 components determined by chromatographic analysis (gas phase components) and true boiling point analysis (liquid phase components) that were subsequently grouped into eight effective components, as summarized in

Table 3. Composition of the reservoir fluid after component grouping.

Component	% mol
N2	31.588
CO2	0.612
H2S	5.085
C1	19.353
C2	3.567
C3–C6	11.99
C7–C11	12.27
C12+	15.5

.

The model was calibrated using the regression method to the following experimental data: the pressure at the saturation point, constant mass expansion tests (relative volumes of gas and liquid phases, effective compressibilities), differential depletion tests (oil and gas formation volume factors, gas-in-oil solubility, oil density, gas z-factor, oil and gas viscosities, gas gravity), separator tests (gas–oil ratio, gas gravity, oil formation volume factor), and viscosity tests (oil viscosity).

The following regression parameters of the SRK EoS were used: critical temperatures and pressures, eccentricity and Ω_A factors of 2 grouped components (C₇–C₁₁, C₁₂₊)—the complete set of the EoS parameters and their values are given in

Table 4. Equation of state (EoS) parameters of the reservoir fluid model for the Soave–Redlich–Kwong (SRK) model.

Component	Critical Temperature Tc [K]	Critical Pressure Pc [bar]	Eccentricity Factor ω	Parameter ΩA	Parameter ΩB	Molar Mass M	Boiling Point Tb [K]	Critical Volume Vc	Critical Gas Compressibility Factor Zc	Parachor
N2	126.2	33.9	0.0400	0.4275	0.0866	28.0	77.4	0.090	0.2905	41.0
CO2	304.2	73.8	0.2250	0.4275	0.0866	44.0	194.7	0.094	0.2741	78.0
H2S	373.2	89.4	0.1000	0.4275	0.0866	34.1	213.5	0.099	0.2837	80.1
C1	190.6	46.0	0.0080	0.4275	0.0866	16.0	111.6	0.099	0.2874	77.3
C2	305.4	48.8	0.0098	0.4275	0.0866	30.1	184.6	0.148	0.2847	108.9
C3–C6	453.2	34.7	0.2315	0.4275	0.0866	65.4	296.9	0.299	0.2752	221.2
C7–C11	641.1	27.3	0.3182	0.4221	0.0866	120.6	418.2	0.631	0.3231	347.6
C12+	784.1	17.3	0.4975	0.4221	0.0866	234.6	575.6	1.180	0.3135	626.1

and

Table 5. Fluid model binary coefficients for the Soave–Redlich–Kwong (SRK) equation of state (EoS) model.

	N2	CO2	H2S	C1	C2	C3–C6	C7–C11	C12+
N2	-	-	-	-	-	-	-	-
CO2	−0.0315	-	-	-	-	-	-	-
H2S	0.1696	0.0989	-	-	-	-	-	-
C1	0.0278	0.1200	0.0800	-	-	-	-	-
C2	0.0407	0.1200	0.0852	0.0000	-	-	-	-
C3–C6	0.0808	0.1200	0.0655	0.0000	0.0000	-	-	-
C7–C11	0.0928	0.1006	0.0006	0.0000	0.0000	0.0000	-	-
C12+	0.0928	0.1006	0.006	0.0000	0.0000	0.0000	0.0000	-

. Additional regression parameters were those of Lohrenz–Bray–Clark correlation coefficients—given in

Table 6. Coefficients of the Lohrenz-Bray–Clark viscosity model.

a1	a2	a3	a4	a5
0.4703	−0.1017	0.0585	−0.0408	0.0093

.

6. Model Calibration

The model of the analysed reservoir was calibrated against the production data covering 16 years of operation with 11 producing wells. The data included daily oil, gas, and water production rates from individual wells, bottom-hole pressures, and well test results.

6.1. Calibration Results

Model calibration involved modifying the following model parameters: global and local petrophysical parameters (absolute and relative permeabilities and permeability anisotropies), well productivity indices, and skin-effect coefficients. As a result, a significant increase in the gas–oil-ratio (GOR) during operation was reproduced for wells C-1, A-1, and A-4, which were located directly or in the vicinity of the primary gas cap (Figure 8), and for the wells located relatively far away from the primary gas cap but influenced by the formation of a secondary gas cap.

The necessity to accurately reproduce the measured bottom-hole pressures required a linear characteristic of the relative gas permeability at the top zone of the structure, which suggests the existence of micro-fractures in the reservoir. The critical oil saturation was found to be consistent with the results obtained from the laboratory experiments regarding oil displacement with water. Further, the permeability vertical-to-horizontal anisotropy increased, which resulted in the correct depression in the production wells and reduced the migration of the released gas to the top of the structure.

Acid treatments of wells C-1, C-2k, and A-4 were modelled by modifying the absolute permeability and effective pore volume in their near-wellbore zones. To reproduce the relatively high pressure difference measured between the C-1, C-2k, and C-11k wells, deteriorated petrophysical properties were introduced among these wells. A local reduction in the absolute permeability anisotropy near the A-11H and A-13k wells caused the simulated results of the GOR to become closer to the observed values. In general, the productivity indices and skin-effect coefficients of the wells were accordingly adjusted to the depressions observed during the well tests.

Figure 9 and Figure 10 present the total oil and gas production from the field as a result of the calibrated simulation model versus the measured values, revealing that the model matched the production data with high accuracy. At the well level, Figure 11 shows an example of the adjustment quality of the bottom-hole pressures measured at a production well. Meanwhile, as shown in Figure 12, good GOR fitting in an individual well was achieved. Further, Figure 13 shows an example of the adjustments to the bottom-hole pressure evolution reported during the well production test. Thus, the simulation results are in good agreement with the measured values for all the available historical production data.

6.2. Model Characterisation after Calibration

After the calibration process, the simulation model of the analysed reservoir was characterised by the following parameters: total area: 15 km × 10 km, the nature of the model: single porosity and permeability, lateral dimensions: 376 blocks × 242 blocks, lateral sizes of the blocks: 100 m × 100 m, layer structure: 15 layers, number of active blocks: 20,325, initial contact depths: oil–water: 3282 m b.s.l., gas–oil: 3178 m b.s.l., initial pressure: 430.2 bar (@ 3282 m b.s.l.), reservoir temperature (constant): 126.8 °C, total model pore volume: 70.2 million m³, average values of basic parameters: porosity: 9.6%, horizontal permeability: 34.91 mD, vertical permeability: 1.4 mD, and total thickness of a simulation layer: 2.77 m.

7. Pressure Evolution

The dynamic flow model of the analysed reservoir was used to simulate the reservoir behaviour during enhanced oil recovery with CO₂ injection followed by exclusive CO₂ injection. Production rates of the oil producing wells were assumed at the nominal levels provided by the reservoir operator. In addition, they were constrained by the minimum well-head pressure, maximum water cut, maximum gas oil ratio, minimum economic flow rate, and maximum allowable drawdown. Initially, CO₂ injection was performed using four wells (A-4, A-6H, C-1, and C-4), before injection expanded to other wells, which were converted from producing wells after they ceased production. The injection rates of individual wells were controlled by the total injection rate and well contributions proportional to their injectivities. The injecting wells were constrained by the maximum bottom-hole pressures not exceeding the formation fracturing pressure. The injection phase lasted approximately 36 years. As a result of the simulation of this EOR/CO₂ sequestration project (called the basic scenario), all the significant quantities characterising the project were obtained. Figure 14 shows the total oil production and average reservoir pressure evolution during the project.

The contributions of individual wells to the reservoir oil production are shown in Figure 15 in terms of their total oil production.

The quantitative characteristics of the CO₂ injection process were expressed as the total CO₂ injection of the project and as the total CO₂ injection of individual wells, as shown in Figure 16 and Figure 17, respectively.

Figure 18 shows the effects of EOR with CO₂ injection by comparing the total oil production of this project with total oil production using the primary oil production method. The incremental oil production exceeds 6.3 × 10⁶ Sm³, which is equivalent to approximately 130% of the primary production (4.8 × 10⁶ Sm³).

The simulation results include pressure distributions across the analysed structure extended model at various stages of the EOR/CO₂ sequestration project, as shown in Figure 19.

To present a more detailed variation in the pressure distribution, pressure maps for the top layer of the reservoir are shown in Figure 20 for various stages of the project.

8. Caprock Sealing Analysis

The migration of stored CO₂ into the overlying strata is inhibited by the sealing properties of very low permeability and high capillary forces, which prevent gas penetration into the water-saturated caprock. The parameter governing this sealing ability is the threshold pressure. If the value of the threshold caprock pressure exceeds that of the gas, the caprock pores begin to be penetrated, ultimately leading to the development of interconnected pathways for gas to escape upward [47]. This threshold pressure can be measured using the method first proposed by Thomas et al. (1968) [48], which was later modified using empirical models adopted for a specific lithology. In general, among the different lithologies, anhydrites, carbonates, and shales have high threshold pressure and low permeability, making them potentially good sealing rocks for gas storage.

In this study, to calculate the threshold pressure in the basal anhydrite (A2), a direct sealing layer neighbouring the reservoir rock, we used the empirical relationship as follows:

$${\mathrm{P}}_{\mathrm{th}}=a{ \mathrm{k}}^{\mathrm{b}},$$

(8)

where ${\mathrm{P}}_{\mathrm{th}}$ is the threshold pressure, $\mathrm{k}$ is the permeability, $a$ is the multiplication coefficient, and $\mathrm{b}$ is the exponent. The coefficients $a$ and $\mathrm{b}$ were developed by Ibrahim et al. (1970) [49], who studied seven samples of anhydrites. The permeability in the history of reservoir exploitation and forecasting of CO₂ injection into the main dolomite reservoir rock was estimated based on the stress dependency on permeability, wherein the stress variation was obtained using geomechanical simulations.

Numerous scientific investigations have revealed that both porosity and permeability are affected by dominant loading conditions and the history of stresses within the sedimentary basin. This relationship is widely known and well documented [13,50,51,52]. Based on the laboratory permeability measurements performed by Morrow et al. (1984) for fault gouges [53], Shi and Wang (1988) suggested that the power law describes the nonlinear decline in permeability with increasing effective stress in the low-permeability rocks well in the following form [50]:

$$\mathrm{k}={\mathrm{k}}_0{\left(\frac{\mathsf{\sigma }}{{\mathsf{\sigma }}_0}\right)}^{\mathrm{p}},$$

(9)

where $\mathrm{k}$ denotes the permeability under effective stress $\mathsf{\sigma }$, ${\mathrm{k}}_0$ is the initial permeability at the initial effective stress ${\mathsf{\sigma }}_0$, and $\mathrm{p}$ is the material constant. Regarding effective stress, we used the magnitude of the minimum horizontal stress as it is mostly responsible for the decline in permeability in the mechanism of closing microcracks potentially present in the anhydrite rocks.

As the data regarding the petrophysical characteristics of the Zechstein anhydrites in Poland were not available, we used a typical value of permeability for anhydrites [54] and the material constant p determined by Zivara et al. (2019) [13]. The results of geomechanical simulations performed for the selected time steps provided the distribution of effective stresses evolving throughout the CO₂ injection operation. The calculated permeability was then used to estimate the threshold pressure at selected time steps, ending approximately 80 years after the initial CO₂ injection.

Table 7. Assumed values of parameters used in the calculation of permeability versus effective stress, and coefficients for calculation of threshold pressure in the anhydrite caprock.

Parameter		Assumed Value
Coefficient $a$ (Equation (8))	[49]	2.6 × 10−7
Exponent $\mathrm{b}$ (Equation (8))	[49]	−0.348
Initial permeability ${\mathrm{k}}_0$ [m2] (Equation (9))	[54]	9.6 × 10−21
Material parameter $\mathrm{p}$ (Equation (9))	[13]	0.6288

lists the assumed values of initial permeability, material constant, and the other coefficients employed for the threshold pressure calculation in the basal anhydrite (A2).

As shown in Figure 21, the evolution of pressure in the selected location at the top of the reservoir Ca2 (red line) and in the basal anhydrite A2 (blue line) exhibit a negative correlation between the magnitude of the minimum horizontal stress and reservoir pressure, a positive relationship between the magnitude of the minimum horizontal stress and the threshold pressure and the evolution of the threshold pressure (green line), and a pressure step between the top of the reservoir and basal anhydrite (orange line).

The spatial distributions of the pressure (A), magnitude of the minimum horizontal stress (B), estimated permeability (C), and calculated threshold pressure (D) for the time steps representing the end of hydrocarbon production (2020—left column), during the process of CO₂ injection (2056—middle column), and at the end of CO₂ injection (2100—right column) in the sealing anhydrite (A2) are shown in Figure 22.

These distributions show that the increase in pore pressure in the sealing rock along with the injection of CO₂ in the reservoir zone caused a very slight increase in effective stress and, consequently, a small increase in the permeability, from 9.55 nD to 9.80 nD. This change results in a small decrease in the threshold pressure, from 24.16 bar to 24.05 bar. These changes are small and local, and the threshold pressure remains at a safe level of approximately 24 bars. The most sensitive area is localised within the slope of the platform, near injection wells B-4 and B-5. The increase in permeability (yellow-red colour) and decline in threshold pressure (purple colour) were significant in this area. Another sensitive zone is located in the deep-water zone near injection wells A1, A2, and A-4.

The difference between the pressure in the top layer of the reservoir (Ca2) and the sealing layer of the basal anhydrite (A2) shown in Figure 23 helped to determine where additional CO₂ can be injected without exceeding the threshold pressure in the anhydrite, thereby lowering the risk of CO₂ leakage.

9. Gas Leakage Analysis

In addition to caprock sealing conditions, the quantitative studies of gas leakage to the caprock were performed using the structure simulation model described above. As the caprock maintains its sealing properties up to the determined threshold pressure, the leakage occurs only when the pressure step across the caprock-reservoir boundary exceeds the threshold pressure. In this section, the quantitative results of the leakage to the caprock were determined for several scenarios, characterised by varying amounts of injected/sequestrated CO₂ as listed in

Table 8. Total CO₂ injection/sequestration capacity.

Scenario	Total Injection/Sequestration Capacity [×109 SCm3]
Basic	12.01
P490	12.17
P500	12.26
P510	12.38
P520	12.43

. The varying amounts of total injection/sequestration were realised by the modifications of the well group injection rate. All the other constraints were kept unchanged.

The progress of the CO₂ injection process was presented in Figure 24 and Figure 25 in terms of time evolution for total CO₂ injection and average reservoir pressure, respectively.

The detailed pressure distribution at the reservoir top revealed a pressure step across the reservoir–caprock boundary exceeding the threshold pressure at several locations. As a consequence, CO₂ leaked from the reservoir to the overlying anhydrite layers. The effects of this leakage process are shown in Figure 26 for each of the four scenarios considered. The final (140 years after the injection finish) amounts of the total leakage and their value as a fraction of injected CO₂ are listed in

Table 9. Total leakage of CO₂ to the caprock.

Scenario No.	Total Leakage [×106 SCm3]	Total Leakage as Fraction of Injected CO2 [%]
P490	0.98	0.0080
P500	1.20	0.0098
P510	1.39	0.0113
P520	1.59	0.0128

.

Although the CO₂ leakage rate gradually decreased with time, it continued up to the extended end of the simulation period (140 years of relaxation phase after the injection finish). Note that the leaked CO₂ accumulated in the bottom layer of the anhydrite horizon, wherein a relatively small amount (2.5%) of CO₂ reached the anhydrite top, as shown in Figure 27.

The process of CO₂ leakage from the reservoir to the caprock was very heterogeneous because of the following:

–

inhomogeneous reservoir rock properties,

–

varying depths of the reservoir-caprock boundary,

–

inhomogeneity of the CO₂ injection process.

The detailed distributions of CO₂ saturation at the end of the simulation period found in the top layer of the reservoir and the three layers composing the anhydrite horizon are shown in Figure 28 for the selected scenarios.

10. Summary and Conclusions

The studies described herein address the practical problem of geological structure tightness for the purpose of CO₂ sequestration applied to a geological structure to be assessed as the potential site of a future sequestration project. The studies took the advantage of advanced methods, realistic assumptions, and quantitative results of previous, general studies to practically evaluate the operation of the real geological structure of a domestic oil reservoir including the initial phase of enhanced oil recovery by CO₂ injection and final phase of CO₂ sequestration.

In particular, we performed an analysis of the tightness of the basal anhydrite (A2), which is the sealing layer for the hydrocarbon accumulation in the main dolomite formation (Ca2). To this end, we used numerical methods that combined geomechanical and reservoir fluid flow modelling. Based on a large set of data, geological, structural, and parametric models were developed. To predict the reservoir and surrounding rock behaviour caused by CO₂ injection, dynamic flow and geomechanical models were constructed based on the geological model supplemented by all the other required components. These models were satisfactorily matched with the historical reservoir operation data. Then, the models were used to forecast the reservoir behaviour (evolution of detailed pressure and fluid saturation distributions) under the assumptions of an EOR project with CO₂ injection followed by a CO₂ sequestration phase. The dynamic simulations were combined with the geomechanical simulations using a two-way coupling procedure and utilising correlations between the geomechanical state (stress and strain distribution) and transport properties of the reservoir rock and caprock, including the threshold pressure at the reservoir–caprock boundary.

The results of the procedure determined the sequestration capacity of the structure. In addition, the effects of CO₂ leakage from the reservoir to the caprock were simulated in cases where the reservoir pressure at the structure top (reservoir–caprock boundary) exceeded the limits of the threshold pressure. The long-term simulations resulted in an assessment of the total amount of CO₂ leakage and its distribution within the caprock horizon as a function of time for various total CO₂ injection volumes above the sequestration capacity.

Based on the simulation results, the following conclusions can be drawn:

(1)

General conclusions:

–

The method applied in the studies proves the necessity to employ an extended model of the analysed structure, wherein the geomechanical and dynamical simulations allow precise estimations of the threshold pressure and provide information regarding critical locations at the reservoir–caprock boundary where leakage could occur.

–

The precise determination of the threshold pressure (with inhomogeneous distribution across the reservoir–caprock boundary) and its evolution with time are crucial factors for estimating the sequestration capacity of the structure.

–

The following two correlations are key factors when the sealing properties of the reservoir caprock boundary are evaluated:

• The correlation between the geomechanical state (stress field) and transport properties (permeabilities) of the caprock;
• The correlation between the caprock permeability and threshold pressure at the reservoir-caprock boundary.

(2)

Conclusions specific to the analysed geological structure:

–

The determined threshold pressure revealed the potential CO₂ sequestration capacity of the structure, showing that it could safely store approximately 12 × 10⁹ Sm³ of gas.

–

The relatively low (up to 3.6%) excess over the determined sequestration capacity resulted in a very small total CO₂ leakage (0.13‰ of the sequestrated volume) up to approximately 100 years of relaxation phase after CO₂ injection is complete.

–

Most of the leaked CO₂ accumulates in the bottom part of the caprock.

–

The leakage process does not cease even at the end of the simulated (100 years) relaxation phase.