Integrated Three-Dimensional Structural and Petrophysical Modeling for Assessment of CO₂ Storage Potential in Gas Reservoir

Pakistan is confronted with significant challenges in mitigating its carbon footprint, attributed to escalating carbon dioxide (CO₂) emissions from activities like deforestation and fossil fuel consumption. Despite advancements in renewable energy, fossil fuels continue to dominate Pakistan’s energy landscape. Therefore, the adoption of CO₂ geological sequestration (CGS) technology holds immense promise for emission reduction, especially in sectors such as energy and manufacturing. Pakistan boasts substantial potential for storing atmospheric CO₂, particularly in subsurface oil and gas reservoirs in sedimentary basins. This presents a unique opportunity for effective carbon sequestration, fostering a sustainable ecosystem. To achieve efficient CO₂ sequestration, it is imperative to identify geological formations with suitable mechanisms for CO₂ injection and containment over extended periods to prevent migration and leakage [1]. Thus, conducting a comprehensive examination of the geological characteristics, including the arrangement, layering, and physical properties of rocks, is essential to assess the suitability of reservoirs units for CGS in Pakistan.

Three-dimensional (3D) geological modeling has emerged as a crucial tool for evaluating the feasibility of CO₂ storage in various geological formations [2]. It encompasses structural modeling and property modeling. Structural modeling constructs models utilizing a combination of geometric entities, that is point, line, surface, and body, emphasizing the shape of geological structures and their interrelationships [3]. Property modeling subdivides the 3D space into regular or irregular voxels, focusing on the spatial distribution of physical and chemical properties, which can be easily integrated with geophysical fields [4-6], geochemical fields [7-9], uncertainty analysis [5, 10], finite-element methods [5, 11, 12], structural analysis [12-14], metallogenic analysis [15-17], and reservoir characterization [18-20]. These 3D geological models amalgamate various data sources, for example well logs, seismic data, and core analysis, to create precise representations of subsurface structures and rock properties [21, 22]. Three-dimensional reservoir models can be generated by accurately distributing petrophysical characteristics using geostatistical methods. Three-dimensional geological modeling holds promise in optimizing CO₂ utilization in petroleum industries [23]. By integration of geological, geophysical, and geoengineering data, accurate 3D models can be created. These models are indispensable for assessing suitable CO₂ storage reservoirs and determining storage potentials at specific sites, providing crucial insights into potential, variability, and injection purposes, aiding in informed decision-making regarding site selection and project economics [24]. Some studies have provided empirical evidence on the efficacy of 3D geological modeling in evaluating the potential of various geological formations for CO₂ storage, for example unmineable coal fissures, depleted oil and gas reservoirs, and saline aquifers [24-26]. Apart from assisting in site characterization and selection, these models facilitate the estimation of effective and CO₂ storage potentials, crucial for risk assessment and project planning [2, 27].

Previous research has extensively examined the Kadanwari Gas Field (KGF), situated in the Middle Indus Basin (MIB), with a specific emphasis on its hydrocarbon production and exploration capabilities. Jadoon et al. [28] examined the stratigraphic and structural capture mechanisms that contribute to the accumulation of hydrocarbons in the KGF, with a particular emphasis on the intricate interaction between faults and dip closures [28]. Ali et al. (2022) conducted a detailed petrophysical evaluation of the reservoir quality and fluid properties of the Lower Goru sands in the KGF, providing valuable insights into the reservoir characteristics [29]. The previously mentioned investigations have enhanced knowledge regarding the geological structure and reservoir characteristics of the KGF, thereby establishing the groundwork for subsequent efforts in exploration and development. Despite the prospective advantages associated with CGS attempts in Pakistan, the country has not yet performed CGS projects at the local or regional level. Prior research in the MIB has predominantly focused on the exploration and extraction of hydrocarbons [30, 31]. However, limited research has been conducted on the possibility of CO₂ storage in this area. Kazmi and Jan [32] presented initial approximations of the potential of gas and hydrocarbon fields in Pakistan to store CO₂. Among these fields includes the KGF, which is situated in the MIB. Their findings indicated that the KGF’s short distance to industrial CO₂ sources and its favorable geological conditions made it a possible location for CO₂ storage.

To perform a complete assessment of the potential for CO₂ storage in the Lower Goru sand intervals, more precisely between the G-sand and E-sand zone of reservoir of the KGF, this research attempts to address a significant knowledge gap by constructing detailed 3D static structural and petrophysical models of the reservoir. Using a methodology that integrates well logs and 3D seismic reflection data, this study focuses on CO₂ storage potential to establish a foundation for future research and applications. By providing a detailed assessment of CO₂ storage potential, this study will help advance the utilization of oil and gas reservoirs in the KGF for CO₂ sequestration, an area that has previously been dominated by hydrocarbon exploration. The primary objective of this study is to estimate important petrophysical parameters for the evaluation and strategy of CO₂ storage by utilizing these datasets. The methodology integrates the sequential Gaussian simulation (SGS) to produce numerous iterations of 3D petrophysical property models. This approach accounts for the heterogeneity and inherent uncertainty present in the reservoir formations. It is expected that this methodology will yield a more accurate representation of the geological conditions, thereby enhancing the estimation of CO₂ storage potential in the KGF situated in the faulted area of the MIB.

Geology and stratigraphy are crucial in the context of CCS initiatives, as they significantly influence the feasibility, safety, and long-term success of CO₂ storage projects. A comprehensive understanding of the geological and stratigraphic characteristics of the area enables precise site characterization and sustainable CO₂ storage. The KGF, located in the Sindh Province of Pakistan within the MIB, is situated on the southeastern side of the Jacobabad High at coordinates 27°10′ latitude and 69°19′ longitude (Figure 1). Discovered in 1987 and commencing production in 1995, the field’s development has been driven by its favorable structural and stratigraphic framework, which is essential for both hydrocarbon exploration and potential CO₂ storage [32].

The structural framework of the KGF was shaped by three major tectonic events, which resulted in the current faulted architecture of the reservoir. Of particular importance for CO₂ storage are the normal and strike-slip faults that compartmentalize the reservoir into distinct blocks. The first event, occurring during the Late Cretaceous, involved substantial uplift and erosion oriented on NNE–SSW, although the extent of rotation remains a topic of debate. This event is evident in seismic sections as a truncated base Tertiary unconformity [33]. These faults may act as both barriers and potential migration pathways for CO₂. The second tectonic event, characterized by the Late Cretaceous uplift and erosion, coupled with the Late Paleocene wrenching, contributed to the formation of structural traps and flower-type structures, thereby enhancing the potential for CO₂ containment [34]. The third tectonic event, involving the uplift of the Jacobabad-Khairpur High, has further influenced the stratigraphy, creating combined structural-stratigraphic traps that augment the storage potential in the Lower Goru sands [13]. The Cretaceous strata in the Kadanwari area, particularly the Sembar Formation, consist of shales rich in organic material [35]. These Cretaceous shales in the MIB are recognized as the primary source of gas [36]. The KGF reservoir comprises a sandy Lower Goru Segment, sealed by a shaly Upper Goru Sequence (Figure 2). The sands, deposited along the shoreline, exhibit an NS orientation [37]. The Lower Goru Formation, which serves as the primary reservoir, consists of alternating sand and shale packages. The sandy Lower Goru interval, characterized by high porosity and permeability, is sealed by the shales of the Upper Goru Sequence, providing the necessary trapping mechanism for CO₂ storage. The organic-rich shales of the underlying Sembar Formation serve as the source rock, while the Lower Goru sands, deposited along the shoreline, provide a favorable reservoir for CO₂ injection.

The seismic data used in this study comprise of 10 km² post-stack, time-migrated seismic reflection cube stored in the SEG-Y format. Detailed information about the 3D seismic cube, including in-lines, crosslines, latitude and longitude range, trace, and amplitude range, are presented in online supplementary Table S1 and S2. The well-log data, along with the well headers, were provided in digital logging ASCII standard and notepad file formats, as indicated in online supplementary Table S3 and S4, respectively. Before proceeding with seismic data analysis, the available 3D seismic cube in SEG-Y format and wireline log data were integrated into Petrel modeling software launched by Schlumberger’s in 2021 for building 3D structural and seismic attribute models. The Interactive Petrophysics (IP) software package (2021) was used for petrophysical analysis. The dataset underwent quality control and harmonization in a defined database.

A comprehensive methodology was employed to enhance understanding of reservoir characteristics by combining 3D structural and petrophysical modeling. This technique integrated geological models with structural constraints, seismic attribute implications, and petrophysical properties. The methodology, depicted in Figure 3, involves interpretating seismic and well-log data, followed by reconstructing 3D models. In this study, petrophysical parameters were simulated using the SGS algorithm and static models were developed to incorporate pore volume calculations. Finally, CO₂ storage potential estimation and its probability distribution were conducted for the estimated storage potential of CO₂.

The Cretaceous zone between the G-sand and E-sand in the KGF, MIB, was analyzed using available 3D seismic data. Conventional seismic interpretation involves two main steps, interpretation and modeling. The interpretation step includes generating synthetic seismograms (online supplementary Figure S1), picking horizons, and interpreting faults. The modeling step involves converting picked horizons to 3D two-way time contour surfaces to interpret structural trends and identify prospective leads. Interpreting geologic horizons is a crucial step in the 3D structural modeling workflow, enabling the inclusion and extension of horizons beyond existing points to define the boundaries of the 3D structural model. The interpreted cross-sectional slices of 3D seismic cube provide a visual representation of the displacement and deformations of strata (Figure 4(a)).

Fault interpretation and modeling are crucial for constraining horizon interpretation. Prior to fault interpretation and 3D fault surface modeling, the precision of geologic fault boundaries was reviewed for accuracy using dip magnitude and seismic edge enhancement attributes. The predominance of normal faults in the Cretaceous stratigraphic sequence was interpreted on dip magnitude and edge enhancement cross-sections by aligning reflector discontinuities with a consistent sense and amount of displacement. Figure 4(b) shows a 3D seismic cube depicting the fault surfaces that impact the breakdown of the stratigraphic architecture. The interpreted fault system in the study area follows a typical pattern, trending from northwest to southeast network, which plays a role in determining the distribution of depositional facies, the compartmentalization of reservoirs, and the upward migration of hydrocarbons.

The technique of 3D structural modeling consists of three primary stage’s fault modeling, gridding (pillar), and horizon generation. The initial stage incorporates the analyzed faults into the structural model as a foundation for constructing the 3D structural grid. The next step, pillar gridding, involves creating the structural grid by building essential pillars based on the faults model. The final phase requires creating vertically arranged layers by inserting the interpreted seismic horizons into the gridded model.

Reservoir parameters were determined through a comprehensive analysis of well logs using the (IP software) package (2021) (online supplementary Figure S2). The analysis calculated the volume of shale (V_sh), total porosity (P_hi), permeability, and water saturation (S_w). Initially, facies were identified at well locations by integrating borehole logs, primarily focusing on gamma-ray (GR) logs calibrated into GR zones. The GR values above 80 API were used to identify shale facies, GR values from 50 to 80 API corresponded to shaly sand facies, and GR values below 50 API indicated the sand facies.

The volume of shale (V_sh) is a critical parameter that directly affects a reservoir’s porosity and permeability, significantly influencing the overall reservoir quality. Accurate accounting for V_sh is crucial to avoid inaccuracies in the dependent properties. In this study, the GR approach has been used to calculate the volume of shale, as follows:

where I_GR stands for the GR index; GR_log represents the reading of the GR log; GR_min and GR_max are the lower and upper limits of GR reading in shale, respectively.

Porosity is a crucial property of reservoir rock that measures the amount of available space that fluid occupies. In this study, the lithology of the reservoir is relatively simple due to the thick layers of sand present. Therefore, we asserted that the density log provides a dependable and well-understood tool for calculating total porosity. While acoustic and neutron logs are also common methods for reservoir porosity estimation, we preferred the density log over them for several reasons in the case of the Lower Goru Formation reservoir. The Lower Goru Formation reservoir is primarily composed of sandstones with a uniform mineralogy and texture. As a result, the density log, which measures the bulk density of the formation, provides a reliable estimate of porosity. In contrast, the acoustic log measures the travel time of sound waves through the formation, which can be affected by lithology and fluid saturation in sandstone formations. Furthermore, the acoustic log can be influenced by variations in the texture of the sand grains, leading to inaccurate porosity estimates. Therefore, in this case, we confidently asserted that the density log provides a more reliable tool for calculating total porosity, as follows:

where $∅T$ denotes the total porosity, ρ_b, ρ_ma, and ρ_f represent the bulk density (as measured by a density log), the density of the rock matrix, and the density of the saturation fluid, respectively.

Effective porosity is essential to determine because it significantly affects the permeability estimates. Neutron density porosity data were used to calculate the total porosity value, as follows:

where V_sh is the shale volume, $∅eff$ is the effective porosity, and $∅T$ is the total porosity from the neutron-density log.

The water saturation (S_w) is a measure of the pore volume filled with water. To calculate S_w, various models are available such as Archie’s Equation, which works best for clean sand, the Indonesian Equation, which is suitable for “dirty formation,” and the Dual Water Model, which depends on the lithology. As thinly interbedded shales can reduce the quality of the Lower Goru sand intervals, the Indonesian equation was used to estimate S_w, as follows:

where S_w is the water saturation, R_t, R_sh, and R_w are the true, shale, and water resistivity of the respective formation obtained from resistivity logs, m and a are the cementation and tortuosity factors, typically taken as 2 and 1 for sandstones, respectively.

Permeability quantifies the degree to which pore spaces within a rock are interconnected, thereby influencing the rock’s ability to transmit fluids. The current petrophysical analysis, therefore, relies on calculating permeability from established equations. In the absence of core data, the Wyllie-Rose approach is commonly used, as follows:

where K is permeability in mD, $∅T$ represents total porosity, S_Wi indicates irreducible water saturation, and the constants a = 10,000.0, b = 4.5, and c = 2.

Geostatistical modeling of rock petrophysical properties involves two main steps, that is variogram analysis and Kriging estimation. The first step entails identifying and modeling the spatial structure of the petrophysical data using a variogram. This step captures spatial correlations and directional continuity of the property being modeled. In the second step, the petrophysical properties are estimated using simulation methods, such as SGS. The characteristics of the variogram identified in the first step directly influence the accuracy and reliability of the simulation.

In geostatistical modeling, the experimental variogram is not directly applied for simulation. Instead, it is fitted to a theoretical variogram model, which ensures smooth and continuous predictions. The theoretical model describes the spatial correlation structure of the variable being analyzed. In this study, we employed the spherical variogram model, which is widely adopted due to its simplicity and effectiveness in capturing spatial dependence in reservoir properties. The spherical variogram model is described as follows:

where h is the lag distance, C is the sill (the maximum variance), and a is the range.

SGS is a robust and widely employed geostatistical method in reservoir modeling, valued for its flexibility and computational efficiency in generating spatial data simulations. The process begins with transforming the original dataset into a normal distribution, ensuring that the data conforms to the assumptions underlying SGS. Known data points are then integrated into the simulation model at their respective locations. For each unsampled point, neighboring data, including previously simulated values, are identified. Subsequently, Kriging is employed to estimate the value at each unsampled location, utilizing the spatial correlation of surrounding data points and providing an associated Kriging variance that quantifies the uncertainty of the estimate. Using Kriging, an estimate of the value at each unsampled location is calculated based on the surrounding points, along with the corresponding Kriging variance, as follows:

where $Y∗(μ)$ is the Kriged estimate at location $μ$⁠, λ_β is the Kriging weight, $Y(μβ)$ is the known value at locations $μβ$⁠, and n is the number of known data points used in the estimation.

A random residual R(⁠$μ$⁠) is drawn from a normal distribution with a mean of zero and a variance σ²_sk (⁠$μ$⁠), representing the uncertainty at location $μ$⁠. The Kriging variance σ²_sk is calculated as follows:

where c(0) is the variance of the random function or the nugget effect, Y_α is the Kriging weight, and c(⁠$μ,μα$⁠) represents the covariance between the unsampled location $μ$ and the known data locations $μα$⁠.

After calculating the residual, R(⁠$μ$⁠) is added to the dataset, ensuring that the spatial continuity of the simulated values is maintained through covariance checks. The simulation proceeds by randomly selecting the next location, allowing the process to traverse the grid sequentially. Once the entire model has been simulated, the values are back-transformed from the normal distribution to the original data distribution. By performing multiple realizations using different random seeds, SGS enables the generation of a series of spatial distributions, facilitating the quantification of uncertainty in the model.

The potential for storing CO₂ is determined by analyzing the petrophysical property from the potential sites under investigation. Different methods exist to calculate the CO₂ storage potential of structural trapping, capillary trapping, and solubility trapping [38, 39]. This study estimates the CO₂ storage potential of the reservoir, utilizing the MRST (MATLAB Reservoir Simulation Toolbox) software which was launched in June 2016. The estimation of CO₂ storage potential in a geological reservoir using petrophysical properties can be effectively carried out as follows:

where V_PV is the total pore volume (m³), S_wir is the irreducible water saturation, B is the formation volume factor, $ρCO2$ is the density of CO₂, and E is the CO₂ storage efficiency. The total CO₂ storage potential is then calculated by summing the individual grid cell potentials.

Pore volume (V_PV) is the total volume of the pores within the reservoir rock, as follows:

where V_pv is the total pore volume, Vi is the volume of cell i, ϕi is the porosity of cell i, and N is the total number of cells in the reservoir model. Actually, understanding the uncertainty surrounding reservoir geological modeling is crucial for decision-making. In the case of CO₂ storage potential uncertainty, it is crucial to calculate the pore volume. To assess uncertainty, the Petrel software was used to conduct 100–400 simulations, with the variogram parameter modified to calculate the volume of the pores. The results for P10, P50, and P90 were then utilized to assess the risk associated with the CO₂ storage potential.

The depth conversion process has transformed the two-way travel time of the seismic waves into actual depths. The 2D depth maps (online supplementary Figures S3(a) and S3(b)) reveal a prominent vertical structural characteristic in the center, where depths are at their shallowest. This indicated the presence of a geological ridge or an anticline structure, commonly associated with potential gas reservoirs. Gas tends to accumulate at these elevated structural features. The contour lines show dense spacing in specific regions, signifying abrupt variations in elevation, while in other regions, they exhibit sparse spacing, indicating a more gradual incline. The steep gradients are particularly noticeable around the structural high, typical of anticlinal structures. The structural high and the location of the wells suggest that this is an area of interest for hydrocarbon exploration, as gas reservoirs are commonly associated with such features. The influential data visualizations (online supplementary Figures S3(c) and S3(d)) are typically used to assess the quality and reliability of the geological model, especially in areas where data density is sparse or where there are uncertainties in subsurface interpretation. The influential data maps help identify regions where the model is well constrained by available data and areas where additional data or refinement may be needed.

Examination of the alignment and arrangement of these fault surfaces reveals that the tectonic conditions of the neighboring plate boundary significantly influence the geological setting of the KGF in MIB. This influence is continuous, affecting the stratigraphic structure and hydrocarbon accumulation in the reservoir window. The G- and E-sand horizons of the Lower Goru intervals have been interpreted from the seismic reflection data, providing a visual representation of the displacement and deformations during the Cretaceous period in the KGF, MIB. Based on the study of the Early Cretaceous geologic sequence, it is evident that the identified layers (such as G- and E-sands) exhibit normal fault extension-related horsts, half-graben, and a graben structure, as shown in online supplementary Figure S4(a). Field scale analysis indicated that the core part of the 3D seismic cube has undergone more significant deformation compared to the borders of the cube. Consequently, the number and complexity of faults decrease as one moves from the center toward the sides of the field, as shown in online supplementary Figure S4(b).

The concentration and distribution of hydrocarbons in these normal faults are highest in the center section of the seismic cube, where the faults are closely spaced. The fault plane profiles exhibit the cutoffs of the hanging wall and footwall, which are valuable for understanding seal behavior in relation to hydrocarbon potential, as well as assessing leakage and movement of stored CO₂. Nevertheless, the E-sand interval appears to exhibit greater deformation compared to the G-sand interval.

For an accurate assessment of CO₂ storage estimation of the KGF in the MIB, it is crucial to have a deep understanding of its 3D geological structural and stratigraphic features, as well as the dynamic behavior of the reservoir and the presence of tectonic expansion faults. We conducted an analysis and visual representation of the spatial structure of the Cretaceous sequence zones. This included examining the presence of faults, determining the horizontal and lateral extent of layers, and identifying the separation of sequences into compartments. During the computation of the model, deliberate faults are introduced into the 2D cross-sections to limit the horizons at each fault offset. The detailed analysis of the reservoir model in the depth domain (Figure 5), along with the 3D geological model and 2D cross-sections, uncovers that the geological structural complexity is a result of various tectonic phases of deformation, such as poly-phase deformation spanning from the Early Cretaceous to recent times.

The model exhibits a remarkable level of complexity in its structural mechanics, demonstrating both extensional and transformative movements. The intricate mechanics of the structure are driven by a series of normal faults that slope in a distinctive NW–SE orientation. The normal fault patterns display clear characteristics of brittle deformation, with Cretaceous sequences being displaced relative to each other by a specific amount. The cross-sections exhibited a prominent control over structural domains, which primarily included horst and graben from left to right. Furthermore, the thickness of the reservoir intervals exhibits a noticeable increase in the NE direction. The thickness from the E-sand level to the G-sand level decreases towards the center, suggesting the presence of syn-sedimentary tectonics. These variations in thickness and structural features are important considerations for evaluating the CO₂ storage potential of the reservoir. The arrangement and characteristics of the modeled faults in the Cretaceous sequences of the G-sand and E-sand intervals enable us to identify potential migration pathways for CO₂. The geometric pattern of these faults could serve as crucial conduits for the vertical movement and trapping of injected CO₂. The presence of hydrocarbon accumulations in the Lower Goru Formation intervals (such as G- and E-sands) indicates the effectiveness of the reservoir in trapping buoyant fluids. This is a positive indication that the reservoir may also be suitable for the trapping and storage of injected CO₂ through structural, stratigraphic, and residual (capillary and dissolution) trapping mechanisms.

Petrophysical properties are crucial for assessing the CO₂ storage potential of identified storage units within a specific area, as they provide essential information regarding the overall quality of these units. The primary petrophysical parameters identified in this study, as shown in online supplementary Figure S5, include porosity (ɸ), volume of shale (V_sh), water saturation (S_w), and permeability. These findings indicate that specific intervals within the G- and E-sand formations exhibit favorable petrophysical characteristics for CO₂ storage, characterized by adequate porosity, permeability, and relatively lower water saturation.

In property modeling, well-log data are distributed across grid cells, thereby enhancing the scale of the well logs by populating specific properties within the model’s grid cells. SGS was employed to distribute the petrophysical characteristics of the reservoir, for example, shale volume, effective porosity, permeability, and water saturation, among the wells and the surrounding areas. The petrophysical parameters were upscaled and simulated using the petrophysical modeling technique in the Petrel software package.

The variogram analysis presented in online supplementary Table S5 and illustrated in Figure 6 provides a detailed understanding of the spatial continuity and anisotropic behavior of key reservoir properties, employing a spherical model for permeability, porosity, water saturation, and shale volume in three directions. These parameters are critical for accurately capturing the heterogeneity of subsurface formations, particularly in reservoir modeling applications. The nugget effect, which is less than the horizontal distance from the sampling point, is minimized or adjusted to zero. At this stage, the data are processed during the variography stage, with variograms plotted in vertical and horizontal directions using a spherical model for effective 3D reservoir properties within the reservoir area. Variography is the method used to study and describe the spatial variations of parameters, based on a technical variogram that analyzes data variability in relation to distance. Consequently, reservoir properties exhibit greater similarity over shorter distances compared to longer distances. Online supplementary Table S5 displays the results of the variogram for permeability, porosity, water saturation, and shale volume. Following the variogram operation, reservoir properties can be effectively modeled.

Once the log data has been upscaled to the 3D geo-cellular grid, the properties of each cell intersecting the trajectory of a well can be inferred between the wells in the grid. Therefore, each cell contains a numerical representation of the selected property. The parameters of porosity, permeability, shale volume (V_sh), and water saturation (S_w) were assigned to the 3D grid using SGS algorithm specifically tailored for the areas between the G- and E-sand layers within the Lower Goru Formation intervals.

The V_sh model, as shown in Figure 7, depicts the distribution of shale and the corresponding shale volume throughout the mapped reservoir. Variability in shale volume is a crucial factor in assessing the reservoir’s potential for hydrocarbon or CO₂ storage. Formations with lower shale content typically demonstrate higher porosity and permeability, making them more suitable for fluid storage or hydrocarbon extraction. Conversely, formations with higher shale content usually exhibit reduced porosity and permeability, which can impede fluid movement and limit storage potential.

The S_W model, as shown in Figure 8, displays multiple variations characterized by both increasing and decreasing trends. These observed differences can be attributed to the diverse nature of the rock formations in the study area. Water saturation measurement is a critical for CO₂ storage, serving as an indirect indicator of the available pore space for CO₂ injection. Lower water saturation values are more favorable for CO₂ storage, as they suggest higher saturation levels of hydrocarbons or gases, thus allowing greater pore space for CO₂ injection.

The porosity model, as shown in Figure 9, illustrates the lateral heterogeneity in porosity resulting from face changes. The 3D upscaled porosity models reveal a favorable range of porosity values, with porosity increasing toward the central region, where a notable faulted zone of the Lower Goru Formation is located. These models were utilized to assess the characteristics of the reservoir in the study area. The presence of elevated porosity values within the faulted zone is advantageous for CO₂ storage, providing ample pore space for the injection and confinement of CO₂.

Modeling of permeability data of the G- and E-sands employed the same algorithm and processes as those used for porosity modeling, including verifying whether the regionalized variable is stationary and normally distributed. The SGS method was utilized to simulate the upscaled permeability values across the zone of interest, specifically between the G- and E-sands. The results indicated a satisfactory range of permeability in millidarcies (md) for all three wells (Figure 10). Permeability is a crucial factor in CO₂ storage as it governs the injection and movement of CO₂ within the reservoir. Higher permeability values are preferable for effective CO₂ input and dispersion across the reservoir.

Reservoir forecasting models are subject to various uncertainties, with quantifying this uncertainty being one of the most challenging aspects of modeling [40]. To address this issue, we incorporated geostatistical parameters to account for petrophysical property uncertainties, resulting in extensive geological uncertainty operations. We conducted 200 operations to estimate the pore volume using Petrel software by varying the variogram parameters (online supplementary Figure S6). The Gaussian sampling approach was employed to compute the percentile of the pore volume.

Assessing the storage potential of CO₂ in oil and gas reservoirs is easier and more straightforward than in coal bed methane as well as in deep saline aquifers. Moreover, hydrocarbon reservoirs exhibit a discrete nature rather than a continuous one, which enables CO₂ storage in these reservoirs. The potential for storing CO₂ was calculated, which takes into account the summation volumes of the grid pores, the density of CO₂, the formation volume factor, and the storage efficiency coefficient.

The 3D geological model, along with well logs, provides the essential petrophysical properties, such as porosity, permeability, and net pay thickness. The spatial distribution of the CO₂ storage potential within the reservoir was visualized using the MRST software (Figure 11). The visualization reveals that the highest storage potentials are concentrated in the central and northeastern regions of the KGF, where the porosity and net pay thickness are relatively greater. Lower storage potentials are observed in the southwestern and southeastern parts of the KGF, which exhibit lower porosity and net pay.

The assessment of the possible CO₂ storage potential was performed by analyzing the P10, P50, and P90 statistics of petrophysical properties, derived from 200 realizations. Table 1 displays the calculated pore volume and CO₂ storage potential based on the P10, P50, and P90 statistics. The minimum estimated storage potential, P10, is 10.13 million tons, while the maximum estimated storage potential, P90, is 101.54 million tons. Despite the challenges associated with injecting CO₂ into this nontraditional gas reservoir, this study indicates its potential as a suitable site for CO₂ sequestration. Gas reservoirs offer several benefits for CGS, including enhanced oil or gas recovery (EOR/EGR) and a reduction in the rate of CO₂ emissions.

Furthermore, the CGS community is unsure of the specific capabilities of these formations. Therefore, the potential of a reservoir was assessed by evaluating its potential for CO₂ storage. CO₂-EOR/EGR in tight oil/gas reservoirs offers the benefit of providing an extra storage location within the CGS network. Due to constraints in data availability, our understanding of CO₂ fluxes in hydrocarbons reservoirs is limited. However, the findings of this study offer an initial indication of the potential for CGS in unconventional oil- and gas-bearing formations.

The implementation of this comprehensive strategy will be expected to tackle the difficulties posed by a lack of geological information and the complicated structural characteristics in the faulted reservoir areas. Furthermore, the objective of this study is to develop and verify a systematic process for assessing the CO₂ storage potential and its associated uncertainties by the utilization of integrated 3D structural and petrophysical modeling. This approach can also be utilized for evaluating other reservoirs in Pakistan. The findings of this study fill the gap in existing research on CO₂ storage in Pakistan and lay the groundwork for future initiatives in CGS in the country, specifically in the MIB region.

The findings of this study, which estimated the CO₂ storage potential of the faulted reservoir in the Lower Goru Formation of the KGF, MIB, can be contextualized alongside similar basin studies conducted globally. Although no studies in the MIB have directly assessed CO₂ storage potential, previous researches in this basin have primarily focused on reservoir characterization and hydrocarbon exploration of the Lower Goru Formation [29, 41, 42], employing methodologies that differ significantly from those required for comprehensive CO₂ storage assessments. However, earlier study suggested that Pakistan possesses significant CO₂ storage potential, with estimates of approximately 1.6 Gt CO₂ in its gas fields upon depletion, facilitated by a favorable alignment between CO₂ sources and potential storage sites [43]. In comparative studies, a study conducted in the Eastern Potwar Basin of Pakistan utilized a reservoir simulator to model CO₂ flooding under various scenarios, for example natural depletion, CO₂ injection, and water injection, which indicated that the CO₂ injection site could store over 9 billion cubic feet of CO₂, significantly improving in gas recovery [44]. Globally, research on the Shahejie Formation in China also employed 3D petrophysical modeling and probabilistic approaches to assess CO₂ storage potential, identifying a range of 15.6 to 207.9 × 10⁹ tones [45]. Similarly, a study in the Masila Basin, Yemen, integrated static and dynamic modeling to reduce geological uncertainty, revealing storage potentials ranging from 4.54 to 81.98 million tones for static models, while dynamic simulations estimated potentials between 4.95 and 17.92 million tons [46]. In the Bredasdorp Basin of South Africa, static models estimated storage potentials between 0.71 and 1.62 million tons, with sensitivity analyses indicating that pressure depletion significantly enhances storage potential [2]. When compared to these studies, the storage potential estimates for the KGF (ranging from 10.13 to 101.54 million tons) are consistent, particularly when accounting for reservoir complexities such as faults and compartments. However, while the KGF study emphasizes static modeling and geological uncertainties, the dynamic assessments conducted in the Eastern Potwar Basin, Masila, and Bredasdorp studies highlight the necessity for future investigations in the KGF to incorporate dynamic simulations, thereby improving the accuracy of CO₂ storage potential estimates. These comparisons underscore the importance of utilizing probabilistic and dynamic approaches across various basins for more effective CO₂ storage assessments.

In this study, an integrated 3D structural and petrophysical modeling approach was employed to assess the CO₂ storage potential in the Lower Goru sand intervals, specifically between the G- and E-sand zones of the reservoir in the KGF, Pakistan. The results indicate the KGF’s significant potential for CO₂ storage, which could significantly contribute to mitigating Pakistan’s carbon footprint and supporting CGS initiatives. The key findings and conclusions are as follows:

1. The interpretation of seismic data and the reconstruction of a 3D structural model reveal that the KGF’s reservoir system is characterized by a network of faults. These faults, oriented in a NW–SE orientation, significantly influence the distribution of various sedimentary layers, the compartmentalization of the reservoir, and potential CO₂ migration pathways.

2. Favorable reservoir properties were identified through petrophysical analysis of well-log data, indicating that the target intervals may contain recoverable hydrocarbons and are suitable for CGS. These subsurface conditions provide viable storage sites for injected CO₂, thereby aiding in the reduction of greenhouse gas emissions.

3. SGS was employed to generate 3D property models, accurately capturing the spatial distribution of key petrophysical properties, for example, porosity, permeability, shale volume, and water saturation. Uncertainty analysis was conducted by incorporating uncertainties in petrophysical property through variogram parameterization. This is crucial for estimating petrophysical property distributions, which are essential for CO₂ storage potential calculations.

4. The CO₂ storage potential was estimated using a customized formula that accounts site-specific characteristics. The results indicate a storage potential ranging from 10.13 million tons (P10) to 101.54 million tons (P90), with an average potential of 53.58 million tons (P50).

This study represents a pioneering effort to assess the CO₂ storage potential in Pakistan’s reservoirs, particularly in the faulted region of the MIB. The findings provide a foundation for future CGS initiatives in the MIB region and contribute to development of a comprehensive workflow for evaluating CO₂ storage potential in complex geological settings. Future research should focus on addressing the challenges associated with CO₂ injection in faulted reservoirs, such as potential leakage risks and fault reactivation. Additionally, variations in combined pore volume and irreducible water saturation can significantly impact CO₂ storage potential. Future studies should incorporate a range of irreducible water saturation values across multiple realizations to perform a detailed sensitivity analysis, which could improve the understanding of how these parameters interact and contribute to uncertainty in storage estimates, thereby enhancing the robustness of CO₂ storage potential assessments. Dynamic modeling and simulation studies would also be beneficial to investigate the long-term behavior of injected CO₂ and optimize injection strategies.

The dataset used in this study is not publicly available due to a data privacy agreement with the CNOOC Research Institute Co. Ltd. However, it can be obtained from the corresponding author upon reasonable request.

The authors declare no conflicts of interest.

Syed Yasir Ali Shah: Methodology, Software, and Writing—original draft preparation. Jiangfeng Du: Investigation and Data curation. Sayed Muhammad Iqbal: Software and Data curation. Linze Du: Validation and Writing—review and editing. Umair Khan: Writing—original draft preparation. Baoyi Zhang: Conceptualization, Validation, writing—review and editing, and Funding acquisition. Jingqiang Tan: Conceptualization and Writing—review and editing.

This study was supported by grants from the National Natural Science Foundation of China (grant no. 42072326), the Hunan Provincial Natural Resource Science and Technology Planning Program (grant no. 20230123XX), and the Scientific Research Project of Geological Bureau of Hunan Province (grant no. HNGSTP202301).

Special thanks are extended to the MapGIS Laboratory Co-Constructed by the National Engineering Research Center for Geographic Information System of China and Central South University for providing the MapGIS® software (Wuhan Zondy Cyber-Tech Co. Ltd., Wuhan, China).

Figure S1 shows well-to-seismic tie analysis for KD-6 well, integration of well log data and seismic response. Figure S2 depicts petrophysical interpretations of wells KD-06, KD-10, and KD-11. Figure S3 presents two-dimensional depth maps of (a) G and (b) E sands and influential data of (c) G- and (d) E-sand horizons. Figure S4 shows interpreted fault and horizon surfaces along the 3D seismic cube. Figure S5 depicts petrophysical log interpretation and estimation reservoir properties for wells KD-06, KD-10, and KD-11, and Figure S6 shows pore volume uncertainty assessment according to 200 realizations. Table S1 depicts start and end line numbers of the 3D seismic cube. Table S2 demonstrates three-dimensional seismic cube lines in the study area. Table S3 presents utilized wells in the study area. Table S4 shows utilized geophysical logs and their applications in this study, and Table S5 depicts the variogram parameters of permeability, porosity, water saturation, and shale volume in three3 directions.