Undoing the paradigm of lithofacies control on carbonate reservoir properties through advanced rock typing Open Access

Carbonate reservoirs, marine and continental, account for more than 60% of the world’s petroleum resources. They display highly heterogenous properties so that the knowledge base and forecasting capability lags behind that of the more predictable siliciclastic reservoirs (e.g., Burchette, 2012). This is due principally to the metastable mineralogy of carbonate rocks; calcite, aragonite, dolomite, and other carbonate minerals are constantly subject to dissolution, replacement, and precipitation from the time of sedimentation, through shallow burial to deep burial and then uplift. Therefore, carbonate pore systems are complex, multiscale, and multimodal (matrix, vugs, and fractures) and rarely resemble those associated with their original lithofacies. Over the last few decades, teams of petrophysicists, geologists, and reservoir engineers have toiled to unravel such relationships through iterations of carbonate rock-typing methods with variable success (see, for example, Gomes et al., 2008; Rebelle et al., 2009; Zhao et al., 2013; Aliakbardoust and Rahimpour-Bonab, 2013; Skalinski and Kenter, 2013; Rebelle and Lalanne, 2014; Kadkhodaie and Kadkhodaie, 2018; Turkey et al., 2018; Mohammadian et al., 2022; Krivoshchekov et al., 2023).

The main challenge in rock typing is determining which geological attributes are driving reservoir parameters and, equally important, their associated spatial and juxtaposition rules for distribution in the “white” interwell space for static or direct dynamic reservoir characterization. Traditionally, lithofacies attributes such as depositional environment, mineralogy, and texture have been utilized as primary controls because they resonate with geologists’ knowledge and experience. In this approach, petrophysical parameters were mapped to depositional attributes, which provided spatial rules (e.g., Al Awadi et al., 2017; Jeong et al., 2017).

Many mature (brown field) reservoirs have extensive mercury injection capillary porosimetry (MICP) data, which permit the use of pore throat radius distributions (PTRDs) (e.g., Thomeer, 1983). Often, PTRD models are interpreted to contain a measure of diagenetic influence that is nested in the primary controlling depositional attributes or rock types (e.g., Grötsch et al., 1998; Neo et al., 1998; Basioni et al., 2008; Lehmann et al., 2008; Rebelle et al., 2009; Lawrence et al., 2010; Warrlich et al., 2011; Salahuddin et al., 2016; Strecker and Ibrahim, 2017; de Ribet et al., 2018). Clearly, this approach has unresolved limitations in quantitatively defining the role of diagenetic modification to PTRDs and rock types, and hence, predictability (e.g., Francesconi et al., 2009; Skalinski et al., 2009).

Hybrid models of depositional and diagenetic drivers also exhibit mixed results for competing predictability in the log domain and spatial juxtaposition since such drivers typically crosscut and have very different spatial preferences (e.g., Skalinski et al., 2009; Finlay et al., 2014; Saneifar et al., 2015; El Din et al., 2018; Kumar et al., 2019) unless resolved by assigning separate regions (controlled by stratigraphy and/or geometric domains) in the reservoir model (e.g., Saneifar et al., 2015; Skalinski and Kenter, 2015).

To address the above challenges in a more systematic and comprehensive way, Skalinski and Kenter (2015) proposed a workflow that first determines the reservoir type (RT; following Ahr, 2006) to identify geological control (e.g., diagenetic versus texture) and then builds rock types using a combination of pore typing in the core-to-log domain, validated by dynamic data. This approach was successfully applied by Saneifar et al. (2015). These authors jointly investigated a total of four different carbonate reservoirs and observed the generally limited influence of primary lithofacies on the resulting petrophysical groupings and their spatial distribution. Instead, they proposed late burial diagenetic processes controlling most of the petrophysical variation and suggested ways to capture the associated spatial trends and juxtaposition rules.

Recent work by Kenter et al. (2022) on core- and reservoir-scale data sets seems to confirm this. A reverse engineering approach on more than 380 reservoir zones of conventional marine carbonate systems across stratigraphy and from different geographic locations was used, and they observed that the relationship between petrophysics and depositional parameters is weak at best, whereas diagenetic modification in the underlying studies was noted by individual studies but largely ignored for its complexity (Figure 1). A review of a large (>250,000 samples) core-scale data set also failed to identify any statistically relevant global control other than for gross texture (mud supported versus grain supported), mineralogy (limestone, dolomite, anhydrite, and clay), and regional and stratigraphic-specific trends for porosity and permeability distributions (Figure 2). These observations fall mostly in line with those by Ehrenberg (2022), who singled out burial depth and temperature as dominant controls. Texture and mineralogy were interpreted as a secondary control on porosity while attributing significant local to regional deviations to porosity preservation due to hydrocarbon migration and/or early onset of overpressure.

This paper builds on those elements by Saneifar et al. (2015) and Skalinski and Kenter (2015), shared in the public domain, such as (but not limited to) Ahr’s (2006) reservoir-type concepts and multiresolution graph-based clustering (MRGC) electrofacies (e.g., Serra and Sulpice, 1975; Serra and Abbott, 1980; Ye and Rabiller, 2000, 2001; Davis, 2018) and has two main objectives. The first objective is to introduce the advanced rock-typing (ART) approach, and the second is to apply ART to a series of case studies. To cover an as broad as possible set of depositional settings and mineral variation—and therefore propose a universal rock-typing solution—the analyzed case studies include one Late Jurassic shallow water trimineral and marine Arab C-D reservoir, two Lower Cretaceous Barremian-age monomineral and marine Kharaib-2 reservoirs (all in the United Arab Emirates and marine carbonate depositional systems), and one Lower Cretaceous Aptian-age trimineral lacustrine reservoir (Angola) (Figure 3A).

Because this paper deals with a challenging topic that crosses several disciplines and fields of expertise, several topics are out of scope: (1) the description of lithofacies and methods of lithofacies groupings or lumping, (2) the geostatistical techniques for static and dynamic simulations, and (3) detailed analysis of the diagenetic processes responsible for the petrophysical behavior of the case studies. Coding methods for lithofacies descriptions are critical, need to be systematic and consistent, and are part of the Discussion section. The actual descriptions used in the case studies were taken at face value (mostly from established geoscientists and companies), and the paper simply cannot include a detailed synthesis for reasons of space and focus. Stochastic simulation techniques for static and dynamic models have been published elsewhere (e.g., Deutsch, 2002; Caers, 2005; Pyrcz and Deutsch, 2014) and stop short of the generation of petrophysical probability cubes, which are more critical to the ART objectives.

The ART method is a pragmatic and data-driven workflow consisting of five steps, starting at the finest observational scale and progressing to larger scales (Figure 3B). Step 1 tests the influence of lithofacies, diagenesis, and fractures in the core-to-log domain scale (comparable to that described by Skalinski and Kenter, 2015) and sets the vertical scale of the RT designation. Step 2 uses electrofacies to generate petrophysical rock categories (PRCs) that integrate statistically relevant geological data from step 1. Step 3 integrates pore types (PTs) from special core analysis in laboratory (SCAL) and/or borehole image (BHI) logs, and nuclear magnetic resonance (NMR) logs, where available, determine the final PRCs (FPRCs) and set the horizontal scale of the RT designation. Step 4 generates an FPRC probability cube through kriging interpolation or via a reliable porosity cube from acoustic impedance and deploys sequential indicator simulation (SIS) to assess spatial trends and juxtaposition. Finally, step 5 identifies the underlying geological drivers, depositional and/or diagenetic, that explain the spatial distributions and relationships documented in step 3 and that can be used to distribute reservoir properties in static or dynamic reservoir models (not in the scope of this paper).

The Arab C-D Formation case study is used to illustrate the workflow (Figures 4–8). It should be noted that the ART workflow is a collaborative process between carbonate geologists and petrophysicists, with the latter executing most of the analytical core-to-log domain procedures in specialist software applications. The authors have attempted to provide language, rationale, and procedures that both serve competencies and facilitate communication.

The step 1 test, guided by supervised learning (i.e., category of machine learning to handle large volumes of data), investigates the predictability of the lithofacies associations (FAs) within the existing cored intervals using a given well log suite (i.e., gamma ray, density, neutron porosity, P-wave slowness, photoelectric, deep resistivity). The FA designations are typically provided by geologists and are assumed to be unbiased (see Skalinski and Kenter, 2015) and reproducible by others (see example in Figure 4A). Note that the origin and definition of FAs are not within the scope of this paper.

Step 1 predictive analytics are performed using the k-nearest neighbor algorithm (k-NN; Fix and Hodges, 1951; Cover and Hart, 1967), which classifies the different lithofacies per closeness, based on their well log signatures into a multidimensional space (i.e., number of features within the well log suite) and within a specific region defined by a selected number of neighbors (i.e., k number). Because low k and high k can cause overfitting and underfitting, respectively—hence, model instability—it is paramount that the appropriate value be determined (Hastie et al., 2009). A k-fold cross-validation procedure (e.g., Yuan and Yang, 2019) is used to estimate the performance of the prediction model for a given range of k-number values (from 1 to 50, for example). The k-number giving the best prediction score during this evaluation is selected as the optimal one, which is then used to build the predictive model. Once data processing and optimization are completed, the predictive model is built on the training data set and then tested on the validation set.

The resulting scores provide an estimate of the relative influence of lithofacies or better associations (FAs) on the log domain. The targeted prediction score, which has been considered to be statistically acceptable, is set to 80%. Outcomes can be displayed on the vertical axis of the ternary diagram in Figure 3C or visualized as contingency tables, in stacked bar diagram format (respectively, Figure 4B, C). Figure 3C graphically captures the relative contribution to flow of three typical carbonate pore systems in a ternary diagram (following Ahr, 2006): (1) depositional matrix, (2) diagenetically modified matrix, and (3) fractures or high-permeability pathways. The corners of the ternary diagram represent RT end members whose identification greatly simplifies understanding petrophysical behavior and geological influence. Skalinski and Kenter (2015) labeled the “center” area as hybrid or RT4 and argued that RTs could be defined from various data types such as but not limited to core-to-log analysis, petrography, SCAL data (MICP), NMR and BHI logs, or any combination of these. Furthermore, they proposed that RTs are associated with unique spatial trends and patterns. Figure 3C shows FA prediction results—the relative vertical positions—for several published marine reservoirs (Saneifar et al., 2015; Skalinski and Kenter, 2015) as well as for the four case studies analyzed in this paper. The relative horizontal position is determined in step 3, which includes determination of the pore system. Mean prediction scores equal to or higher than 80% suggest that FAs (and therefore their pore system) are strongly linked to the physical properties recorded by well logs (log domain) and should provide the basis for further analytics (blue zone for RT1 in Figure 3C). Mean prediction scores between 40% and 80% imply that both FAs and secondary modifications (i.e., either diagenetically modified matrix and/or fractures/high k pathways) control the log domain (see vertical 40%–80% interval in Figure 3C). Some heterogeneity may be explained by FAs, whereas others may be fully controlled by diagenetic modification (matrix and/or fractures). Note that only one case study, Wafra Eocene 1, is in this zone (Figure 3C). Finally, scores below 40% indicate the absence of influence by lithofacies and suggest that the flow behavior is fully controlled by diagenetic modification of the pore system (see Figure 3C). Step 3 identifies the relative horizontal position (diagenetic matrix versus fracture/high-permeability pathways) in the Figure 3C ternary diagram.

Those FAs determined to be predictable with an 80% (or better) score now labeled FA categories (FACs) will be used to define PRCs further in the process, and each will have a spatial context (proportion and/or trend and juxtaposition). Those FAs falling below this score are lumped together to reach the statistically acceptable threshold of 80% and labeled FACs as well. Such lumped groups are purely statistical and have no geological or spatial context. This iterative approach completes the determination of FACs in the cored intervals, which are statistically robust and predictable by logs (Figure 5A, first two columns). These new categories capture the lithofacies, including the secondary overprints (i.e., diagenesis and fractures) on the depositional pore system. Once determined in cored sections, the FACs are subsequently propagated into the uncored sections in all of the involved wells, using the previously built k-NN predictive model.

Electrofacies analysis, unlike step 1, uses unsupervised learning within the same well log suite. The resulting electrofacies association categories (EFACs) are statistically stable and highly predictable. However, they represent the physical properties of rock only and lack any geologic meaning. The unsupervised learning technique used is the Gaussian mixture model algorithm (e.g., Hadavand and Deutsch, 2020). The principle of Gaussian mixture model is first to consider the data set as a mixture of a finite number of Gaussian distributions, and then tune the key hyperparameters (i.e., mean and covariance) to optimize those distributions. As for all unsupervised algorithms, the finite number of distributions is required as an input (i.e., the number of clusters) and obviously unknown a priori. To determine the optimum number of clusters, the elbow method (see Bholowalia and Kumar, 2014) is used. Once determined, EFACs are nested into the previously defined FACs to integrate EFACs with, if statistically persistent, geologic behavior contained in FAC (Figure 5A, B). The resulting PRCs from FACs and EFACs are statistically stable and are highly predictable by well logs. To link the PRCs to permeability values, a supervised k-NN model is built for permeability prediction. The k-NN model is trained by the well log suite already used in the previous step, supervised by permeability values from conventional core analysis (CCA/routine core analysis). This computes a predicted permeability log for all of the wells involved.

The integration of PTs is a key step for capturing the petrophysical heterogeneities that influence matrix-scale and larger-scale permeability and flow behavior in most carbonate reservoirs. Such heterogeneities are the sum of primary, depositional fabric–related pore systems modified by secondary processes (e.g., compaction, diagenesis, fracturing). Therefore, determination and prediction of PTs is essential for reliable carbonate rock typing, including the identification of spatial trends. As displayed in the ternary diagram (Figure 3C), flow behavior can be matrix driven (depositional or diagenetic) and/or influenced by a secondary large-scale pore system, including vugs and fractures. One or another (or both) pore systems will be investigated, depending on the data scenario (i.e., data available). Two main data types for pore typing are discussed here: MICP data (i.e., core plug data), to capture the core-scale matrix pore system, and specific well log data such as NMR and BHI logs, to investigate the secondary larger-scale pore network. Each case is illustrated in the following paragraphs, and results are embedded in the PRCs. One common objective of MICP measurements is to calculate capillary pressure (i.e., Pc, ability of a formation to retain water) and estimate fluid saturation. Another important application of MICP consists of deriving PTRDs that are tied to flow units and matrix permeability (e.g., Purcell, 1949; Swanson, 1977; Brown, 2015). To capture the full range of pore throats, several requirements are important: (1) MICP samples need to be representative of the reservoir heterogeneity space and only high-pressure MICP data (pressure up to 60,000 psi or ∼414 MPa, PTR down to 0.002 µm) are accepted; (2) data are rejected if the pore volume is lower than 0.5 cm³, but data from samples with pore volumes ≥1 cm³ can be used with confidence; and (3) depth matching between well logs and the selected MICP core plugs needs to be precise and reliable. When those conditions are met, the selected MICP data are used to compute pore size distribution functions, and MRGC (Ye and Rabiller, 2000) is used to determine PT categories. Like the EFAC step, the elbow method is used to tune to the optimal number of clusters. The resulting PT clusters represent distinct PTRDs and relate to the spectrum of flow properties (Brown, 2015). In the core domain, PTs are populated using a k-NN model that is supervised by the PTs (previously defined from MICP data clustering) and trained by permeability and porosity from CCA (e.g., Saneifar et al., 2015). In the log domain, PTs are propagated using a second k-NN model supervised by the populated PTs in the core domain and trained by the predicted-permeability log and porosity log. Figures 6 and 7 show an example of the progression of the analysis with the Pc curves and pore size distribution models following PT definitions from MRGC (Figure 6A), PTs prediction in the core domain (Figure 6B), nesting of the PTs into PRCs following prediction in the log domain displayed as bar diagrams and graphics (Figure 7A, B), and review of resulting FPRC (see details in the next section) proportions in the 250 wells used for this study (Figure 7C). In addition to providing an estimation of porosity and fluid saturation, NMR wire-line log data are sensitive to pore size (not pore throat size). Resulting T_1,2 distributions are estimators of pore size, and their distributions can be captured by applying an unsupervised learning (MRGC) on the T_1,2 well log. However, NMR should be treated with caution regarding the small radial depth of investigation of the tool (R = 0.5–1.5 in. or ∼12.7 to 38.1 mm), and therefore limited sample volume, truncation of relaxation times underestimating larger pores, and reading drilling mud properties where encountering caliper deviation caused by large pores (i.e., caves).

The main objectives of BHI techniques are the determination of PTs (typically not captured by NMR) and PT sizes, in addition to the interpretation of structural and sedimentary features. Especially in carbonates, topology of fractures (density and aperture) and vugs (density and size) are paramount to flow prediction. As for MICP and NMR analyses, BHI PTs are defined by applying unsupervised learning (MRGC). The pore typing results complete the RT identification and position in the ternary diagram (Figure 3C).

The FPRCs are defined by nesting the PTs into the previously defined PRCs that result from FACs and EFACs (see Figure 7B). The PTs could be from MICP (matrix pore system) or well logs (NMR or BHI) or the combination of all, depending on the data scenario. An example of a set of FPRCs is shown in Figure 8. Pie charts for mineralogy and ranges in porosity, permeability, and impedance for the remaining highly predictable FPRCs are shown in Figure 8 and form the basis for a relative ranking of the reservoir quality (Figure 7B). A final assessment of the relationship between FPRCs and FAs can be investigated by mapping the latter to the FPRCs (Figure 7A).

The previous step revealed the relationships, if any, between FPRCs and depositional parameters manifested by the FAs in the core-to-log domain. This step distributes FPRCs in the earth (reservoir) model via a reliable acoustic impedance-derived porosity cube or application of a kriging interpolation, which has minimal human bias. The resulting product is entirely data driven and allows the investigation of FPRC spatial trends, juxtaposition rules and extract processes, and drivers underpinning the spatial distribution, and, finally, the application of those relationships in the realization of one or more static and dynamic reservoir models. To feed the interpolation, upscaled FPRCs are used for vertical wells only because they have sufficient spatial distribution across the area of interest (reservoir domain). To guide the kriging mean, variograms for each FPRC are calculated to estimate their spatial continuity (i, j), and the vertical proportion curve (VPC) is used as a vertical trend (k). The result is a probability cube for each FPRC that shows primary apparent and consistent spatial trends and juxtaposition (Figure 10B, C). Figure 10D is the result of a stochastic simulation using SIS constrained by the FPRC probability cubes. Using multiple seeds did not change the results significantly, which therefore suggests that the resulting cube/zones are statistically stable. However, in this case study (Arab C-D in the United Arab Emirates), it must be noted that other than for the anhydrite FA-1, none of the FPRCs has a meaningful relationship with core-based FAs (Figure 9B) and associated depositional models that provide spatial context published by, among others, Strohmenger et al. (2004, 2006) and Van Buchem et al. (2010).

To review and identify processes that explain the observed spatial trends, a comprehensive review of the FPRC probability cube is required. Figure 10A shows the VPC of the FPRCs for the Arab C-D case study, and the trends very much resemble those of the lithofacies that are driven by depositional sequences resulting from changes in accommodation space. Although the waning and waxing of proportions may indeed be comparable, the FPRCs have, as explained earlier, no link with FAs. Interval 1 and interval 2 in Figure 10A (FPRC vertical proportions) represent very different parts of the reservoir. Corresponding FPRC lateral trends for these selected intervals are displayed as probability maps in the deterministic interpolation (Figure 10B, C). Here, horizontal intersections with the gas–oil and oil–water transitions clearly show that FPRC areas follow zones truncated by fluid contacts and therefore have a preferential occurrence associated with these pore fluids. In fact, the primary control on FPRC distribution appears to be fluid type rather than anything depositional, as documented by the analysis in previous steps. The observation that the best reservoir quality is mainly located near the crest in the hydrocarbon leg and decreases toward the flanks is believed to be a primary diagenetic modification trend associated with the free water level (FWL). This was reported earlier by, for example, Hollis et al. (2017) and Eherenberg et al. (2024) as gas and oil progressively arresting cementation by calcite during migration into the structure, whereas in the water leg, such cementation and associated compaction continued reducing porosity and generating stylolites. Earth models of the FPRCs (Figure 10D) and corresponding porosity, permeability, and water saturation were populated with primary trends and juxtaposition rules that overwhelmingly reflect (late) burial diagenetic modification. Despite the general absence of or, at best, weak, primary depositional trends, the models are very robust because the FPRCs are better than 80% predictable in the 250 wells available and can easily be updated with new well data. Saturations are predicted by FPRCs, and effective porosity can be estimated from those values and directly populated, skipping conventional saturation height modeling. In addition, if desired, the FPRCs can be upscaled to a dynamic grid and supplemented by relevant properties such as relative permeability, fluid properties, and others.

In addition to published work on this topic (Figures 1A, 2), this paper contributes another four case studies to test the geological link with the core-to-log domain and, more important, to investigate the observed trends and spatial juxtaposition of the FPRCs. As explained earlier, one of those case studies, the Middle East Arab C-D, was used for illustrating the workflow in the previous sections (Figures 4–8). The following is a brief review of the four case studies and integration with the four previously published case studies (Figure 3C, modified after Ahr, 2006 and Skalinski and Kenter, 2015) with a rock-typing approach that used similar published elements and concepts (Saneifar et al., 2015; Skalinski and Kenter, 2015). Generic information on all of the reservoirs is summarized in Table 1.

This case study was used as documentation for the generic ART workflow steps in the previous section, and the following will be limited to the lessons learned in terms of factors controlling the multiscale petrophysical behavior. The offshore reservoir is Upper Jurassic (Kimmeridgian) in age and deposited as a prograding marine shallow-water ramp or platform with a suite of lithofacies (and groupings) based on texture, sedimentary structures, and dominant grain types (Grötsch et al., 2003; Strohmenger et al., 2004, 2006; Van Buchem et al., 2010) (Figure 4A). Although regarded as a standard for comparable depositional systems in the Middle East subsurface, no continuous outcrop analogues are available to confirm the interpreted FAs and their spatial trends and juxtaposition. Table 2 summarizes the main results and observations and confirms the weak link (average prediction score is ∼35% and only a single FA has a prediction score better than 80%, anhydrite) between geological attributes and core-to-log domains (Figure 4B). Attempts to predict diagenetic “facies” similarly showed a weak link and failed to reach any predictable threshold. Seven better than 80% predictable FPRCs capture the core-to-log heterogeneity (Figures 8, 9). These FPRCs reveal spatial trends and juxtaposition (Figure 10) that suggest a first-order fluid-type control, followed, in relevance and not in stratigraphic time, by an early diagenetic (dolomite associated) pore system modification, confirming the conclusions published by Morad et al. (2012) and Hollis et al. (2017). Utilizing these spatial rules (fluids and dolomitization) resulted in a reservoir model and associated property distribution that is robust, unbiased, and easy to update.

The upper Thamama B, or Kharaib-2, is of Lower Cretaceous (Barremian) age and is believed to be deposited in a mostly aggrading ramp of platform setting or on an epeiric platform (Ehrenberg et al., 2024). The FAs are, once again, based on earlier classifications by Strohmenger et al. (2004, 2006) and were recently updated by Tendil et al. (2022). It is not within the scope of this paper to discuss the explicit nature of those classifications but rather test whether they hold up (are predictable) in the petrophysical domain. Two different fields with different data scenarios, one with a limited number of cored wells (35 only) available and one with 84 cored wells out of a total of nearly 480 wells, were used for testing and yielded very similar results. The results from the data-dense case study are discussed here. The prediction of lithofacies resulted in an average score of 47%, which confirmed a weak link between lithofacies and core-to-log domains (Figure 11A). Further lumping or grouping reached better than 80% scores for two groups that represent dense argillaceous and reservoir-prone lithofacies (Figure 11B). Subsequent electrofacies generation produced six EFACs, one of which equals the dense argillaceous FAC1 and the other five contained within the reservoir-prone FAC2. Since the reservoir is monomineral—all limestone—the EFACs can be renamed PRCs. Figure 11B shows the progression from lithofacies to FAC to PRCs, and the proportion of the PRCs in the reservoir is displayed in Figure 11C. The MRGC analysis of PTs from MICP resulted in two highly predictable (more than 80%) and distinct PTs (Figure 12A) that were propagated to core and log domains (Figure 12B) and nested in PRCs (Figure 12C). A bar diagram display of FPRCs allows the reduction to eight dominant categories (Figure 12D). Figure 13A and B summarizes, respectively, porosity and permeability ranges for the eight FPRCs and allows a relative ranking of increasing reservoir quality (Figure 13C). Figure 13D displays the relationship between highly predictive FPRCs and original FAs. Because the FAs have a weak link with the core-to-log domain, spatial influence and possible geologic control can be assessed only by kriging of the predicted FPRCs in all of the available wells. Like the Arab C-D in case study 1, the VPC of FPRCs in Figure 14A shows a stratigraphic waning and waxing of reservoir-quality–based FPRCs that resemble or suggest an intrinsic relative sea-level signal. The FPRC7, associated with FAC2 or dense argillaceous lithofacies, is dominating the top and base parts of the Kharaib-2 interval. Both FPRC1 and FPRC10 appear to represent intermediate baffles with a sedimentary origin; however, cross sections of proportions of kriged (in 480 wells) FPRCs in a spatial context suggest a very different interpretation (Figure 14B). Here, pastel color rendering of FPRCs highlights a predominance of FPRC8, FPRC5, and FPRC9 (respectively, from “good” to “best” reservoir quality) in the center and northeastern crest regions, with FPRC5 near the bottom. In an intermediate vertical position and dominant is FPRC8, and FPRC9 is near the top (Figure 14B, sections 1 and 3). These FPRCs pinch out near or above the established FWL. Toward the southwestern part of the structure, a transverse slice of reservoir (Figure 14B, section 2) shows dominance by moderate reservoir quality (absence of pastel green color, FPRC8) and is roughly bound by transverse faults with moderate offsets. These first-order spatial observations of FPRCs suggest a predominant influence of fluid type controlling the distribution of reservoir quality; the oil phase is interpreted as having preserved reservoir from cementation, whereas the petrophysical properties are degraded below the FWL (Figure 14B, sections 1 and 3). More complex diagenetic processes probably have taken place in the fault-bound segment (imaged by Figure 14B, section 2), where the reservoir quality preservation may have been counterbalanced by fluid circulation–related cementation along certain fault segments. Although comparable observations were made earlier by Morad et al. (2012) and Hollis et al. (2017), and more recently by Ehrenberg et al. (2024), the ART method presented in this paper maps reservoir quality (represented by the FPRCs) in all available wells, quantifies any link with FAs, and reveals the spatial patterns and juxtaposition between FPRCs at field scale. In addition, the reported FPRCs will assist in the identification and mapping of the widely published tilting of paleo-FWLs in the Middle East (e.g., Heydari-Farsani et al., 2020). Moreover, the spatial mapping of the FPRCs may also provide a complementary understanding of the structural control on diagenetic overprint.

The Aptian-age Pre-Salt along the South American conjugate margin is typically considered to be a strikingly different depositional system of shallow water lacustrine makeup. Despite many publications, no consistent depositional model has been adopted (e.g., Lima and De Ros, 2019; Herlinger et al., 2020), and the link with reservoir quality remained unresolved for a long time. However, numerous recent publications report and emphasize diagenetic modification as a key parameter shaping the pore system and associated matrix and excess permeability behavior (e.g., Fernández-Ibáñez et al., 2022; Mendes et al., 2022; Quediman et al., 2022; de Lima et al., 2023; Mimoun and Fernández-Ibáñez, 2023; Wennberg et al., 2023). This case study builds on previously published work by Cazier et al. (2014) and Saller et al. (2016) and focuses on the application of ART to identify and quantify the process(es) controlling reservoir quality. Figure 15A shows the depositional model following the concepts developed by Saller et al. (2016) for a set of nearby exploration wells. Because this field is in an early development stage, with only three wells at the time of this project, the core, rotary sidewall core, and cuttings-based descriptions represent less than 10% of the log domain and have very narrow ranges in log properties (Figure 15B), not capturing the full width of petrophysical properties. This not only makes it very challenging to contrast and compare the range of lithofacies and their proportionality (Figure 15C) to the model by Saller et al. (2016) but it also renders any lithofacies prediction meaningless. As a result, the first step of estimating the depositional influence on the core-to-log domain is abandoned for lack of appropriate data, and no statement can be made as to the nature of the relationship. The following step resulted in the generation of nine exclusive electrofacies that consequently can be regarded as PRCs (Figure 16A). Analyses from borehole imagery, lost circulation, and wide caliper observations (“bad hole”) substitute the lithofacies as continuous data and result in the definition of three PTs. These represent conductive fracture density, pore size and shape, and bad holes, which can be predicted by logs (Figure 16B–D) nested in the PRCs (Figure 17A) and lumped together to represent three highly predictable FPRCs (Figure 17B). Figure 17C shows the attribute distributions across these FPRCs. No direct evidence can be found that supports a consistent and statistical relevant relationship with the depositional framework. Instead, the strong ties to PTs, bad hole depositional and mineralogy suggest that diagenetic modification has a predominant role in porosity development and flow behavior of such reservoirs. Preliminary population of FPRCs in the acoustic impedance porosity cube strongly suggests a primary depositional distinction between “shallow” and “deep” crosscut by lineaments of vertical bodies that may well represent fault corridors, none of which are clearly visible in the three-dimensional (3-D) seismic cube.

Evidence has been presented that sheds light on the factors or processes that control carbonate reservoir quality in the subsurface. One is at core scale, one is at core-to-reservoir scale, and one is centered around the rock-typing approach presented in this paper. The global core scale data are numerous (>270,000 data points) and presented in both burial depth versus porosity and porosity-permeability crossplots (Figure 1). In a burial context, the primary control is depth (and likely temperature, as proposed by Ehrenberg [2022]). When filtered (normalized) for burial depth, no single depositional attribute or mineralogy appears to have statistically significant control in the porosity-permeability space. In part, this is due to the absence of systematic data on effective porosity (Rabiller and Gerges, 2020), which would likely improve existing trends by considering the porosity contributing to flow using MICP analysis. In the absence of such data, however, the only attribute that seems to have a mild influence is texture. Mud-dominated textures appear to have lower permeabilities than grain-supported textures at similar porosity values. It should be noted that distinct trends can be present for certain sedimentary rock types in individual reservoirs, although such observations are rather rare.

The integrated core-to-reservoir scale data of more than 380 reservoir zones with static and dynamic models suggest a similar absence of obvious depositional partitioning in flow behavior (Figure 2). Scatterplots of porosity versus permeability for “inner” and “outer” shallow-water depositional domains or “slope” domains show comparable ranges and relationships (Figure 2). When further split by lithofacies, for example, lagoonal deposits show petrophysical property ranges and relationships similar to those of carbonate sand shoals (e.g., dominated by ooids or rudists), and although the rock-typing methods vary greatly among the data sets, one would expect that their number would overwhelm such bias.

The application of variable rock-typing methods is likely the cause of error in the previous argument because the methodology varies from case to case (hence, making them incomparable). For example, many of the methodologies assume that lithofacies should be one of the major controls influencing the reservoir flow behavior and therefore introduce an initial bias in the analysis. However, the rock-typing methods used for the four published case studies (Saneifar et al., 2015; Skalinski and Kenter, 2015) and those in the present study use similar published elements and concepts and are comparable. The data set of eight case studies with variable stratigraphic age, depositional setting, mineralogy, well density, and data quality (Table 1) unanimously suggests that early diagenesis (dissolution by meteoric fluids or dolomitization by hypersaline brines) has a modest influence, whereas late diagenesis (e.g., dissolution, cementation, recrystallization) caused by fluid circulation resulting from tectonic expulsion, hydrothermal venting, or hydrocarbon migration are the primary processes that shape the final petrophysical properties relevant to reservoir characterization. Although fractures and joints are commonly associated with late burial diagenesis, early synsedimentary fractures in carbonates may remain open or even further corrode during successive late burial events (e.g., Narr and Flodin, 2012). Only one case study (Tengiz Platform) suggests a depositional control for one PRC, whereas a second suggests a hybrid scenario (Eocene Wafra) in which dolomite crystal size controls petrophysical behavior while mimicking depositional partitioning. Finally, a third case (Pre-Salt case study) paints a scenario in which the final PRCs suggest a secondary control by “deep” versus “shallow” environments of deposition, but the primary control is deep burial diagenesis related to high-pressure hydrothermal fluids circulating via fracture swarms. The data presented strongly suggest that subsurface pore systems in carbonate reservoirs are primarily the result of diagenetic modification far away from the depositional parameter ranges from micro- to macroscales.

Traditionally, data collection for reservoir characterization has focused on depositional attributes and PT characterization via core and petrographic descriptions supplemented by meter- to hundreds of meters–scale outcrop analogues for calibration and spatial context. Diagenetic studies focused mostly on the micrometer-to-centimeter scale and resulted typically in detailed paragenetic sequences relating main events to diagenetic models available in the public domain (e.g., Wilson, 1975; Moore and Wade, 2013 and references therein), most of which have inferred from observations on spatial trends in modern examples. The observations from this work suggest that although a vast amount of knowledge has accrued on carbonate diagenesis, paragenetic models, and spatial concepts, there is a general mismatch between such models and what is observed in the subsurface resulting from petrophysical rock-typing applications at reservoir characterization scale and smaller (present study). For example, recent work on the Thamama-B reservoir present highly detailed and corroborated observations on diagenetic processes influencing reservoir quality (e.g., Ehrenberg et al., 2024); however, it remains very challenging, if even possible, to apply such findings in reservoir characterization of related reservoirs. It is proposed that focus and resources be shifted to understanding porous media trends defined by highly predictable rock types in the subsurface in a spatial context first, followed by detailed descriptions and inferences on process as a second step. Reactive transport modeling demonstrated itself to be a valuable tool (e.g., Xiao and Jones, 2006; Whitaker and Xiao, 2010; Xiao et al., 2018) for scoping first-order processes and associated trends and has been deployed successfully in Pre-Salt conjugate margin settings (Jones and Xiao, 2013) and Pre-Caspian Basin (Jones and Xiao, 2006). Outcrop analogues, where they provide multiscale mapping of diagenetic products, may very well offer another opportunity. Recent examples of continuous mapping (using hyperspectral imaging) of the dolomite contribution to Jurassic ramp and platform strata in Morocco and Saudi Arabia allowed the determination of process and spatial patterns of dolomitization at seismic scale (e.g., Dujoncquoy et al., 2022; Gairola et al., 2024).

One possible source of uncertainty in the propagation of depositional attributes may be embedded in traditional sedimentological or depositional models and core description, classification, and grouping of lithofacies and coding practices (Galluccio et al., 2022). Although Galluccio et al. (2022) propose a more systematic classification and coding approach for lithofacies and depositional environments, the anchoring in conventional geomodels from outcrops and modern depositional environments is unchanged. Whereas coding needs to be systematic and consistent, the main source of uncertainty is more likely to originate from the conversion of modern-like depositional landscapes, affected by erosion and burial, to spatial partitioning of lithofacies in the subsurface. For example, modern landscapes of ridge-shaped and migrating ooid shoals near Joulter Cays in the Bahamas (Harris, 2009) do not resemble similar deposits in the Bashkirian, which have (interpreted from core) sheet-like morphologies lacking such ridge shapes and have porosity trends generated by dissolution rather than primary interparticle space (Kenter et al., 2010). The confidence attributed to common and long-lived depositional models for lithofacies (e.g., Strohmenger et al., 2004, 2006), although partially based on nearby outcrop analogues, may need a recalibration with densely cored reservoir data sets. The reality, unfortunately, is that the log domain relationship with geologic parameters from core is generally weak (below statistically acceptable thresholds) and therefore obscures petrophysical distributions in resulting models. Diagenetic modification may explain this lack of correlation; however, the current understanding of such processes primarily addresses geochemical and petrographic characteristics and omits spatial trends and juxtaposition rules, resulting in inadequate models. The absence of reliable outcrop analogues is another factor hindering the development of robust spatial concepts. While lithofacies models have been published in large numbers across the stratigraphic record and mostly from fully accessible outcrop analogues supplemented by information from the modern, their accuracy and reliability for subsurface applications may be questioned due to the limited understanding of paleo environments and mostly qualitative approaches used. Finally, commonly used carbonate rock-typing methods are not tailored to test and handle such uncertainties. By focusing on the statistical predictability and coherency of PRCs first and then determining kriged distributions in reservoir space, spatial trends and juxtaposition rules are extracted that form the basis for assigning appropriate processes that reliably explain reservoir heterogeneity.

The ART workflow presented in this paper allows the investigation of reservoir complexity in an unbiased, statistics-based, data-oriented, and repeatable manner, leading to the identification of FRPCs. Such FRPCs can be distributed spatially within a geologic model, and the resulting trends, shapes, and juxtapositions can be interrogated for the underlying controls on these distributions, which can be unique for a given reservoir. When consistently applied, this road map provides a baseline for a systematic comparison of reservoirs and their property ranges, distributions, and behaviors, which can be used in exploration and production optimization of hydrocarbons, as well as carbon capture and storage.

The technical skill required from the user is not excessively demanding; ART can be run on most commercial petrophysical software applications, such as SLB’s Techlog or Petrel, Paradigm’s Geolog. It also requires machine learning applications in a dedicated programming language, such as Python or C++.

Another strength of the workflow proposed here is its time efficiency; the average time required to complete each of our case studies was approximately 3 months for the petrophysicists and geologists combined, using data that were quality checked and ready before the start.

The principal limitation of ART, in the authors’ view, is not technical but rather related to the natural human resistance to change. Carbonate reservoir characterization is a complex exercise integrating diverse data sets, scales of observation, software platforms, and so forth, carried out by multidisciplinary teams over months, even years. Methods and workflows have been developed over many decades, and innovative approaches demanding a shift or modification of such ways of working often find substantial resistance from individuals, teams, and/or corporations. As an example, many geoscientists are hesitant to accept that lithofacies do not control reservoir quality, whereas many studies show that the dominant guide to the spatial distribution of properties in the subsurface is diagenesis (e.g., the 380 reservoir zones studied by Kenter et al., 2022). Even when diagenetic modification is presented as an important control, the traditional sedimentologic description of such modifications cannot be translated readily into quantitative and spatial information required for reservoir characterization.

This paper presents ART as an approach to generate highly predictable rock types with realistic distributions in 3-D earth models. Advanced rock typing is a pragmatic and innovative methodology for carbonate rock typing and a paradigm shift in how the subsurface characterization community addresses reservoir quality in carbonates for the following reasons.

1. The resulting rock types capture only statistically relevant geological attributes, follow data-driven spatial trends and juxtaposition rules, and can be predicted with an accuracy of at least 80%.

2. It confirms the weak correlation between depositional attributes and the core-to-log domain observed by many previous studies.

3. It provides unbiased assessments of the processes driving first-order changes in multiscale and multimodal porous media.

4. The results of the case studies confirm the dominant role of late burial fluids on reservoir quality.

5. This method is data driven, systematic, and fast to deploy, and it allows meaningful comparisons between different reservoirs or between iterations of a given reservoir model.

6. It identifies unbiased modeling strategies for well placement and reservoir optimization in production settings.

The authors would like to acknowledge the TotalEnergies students and interns who have contributed through various projects: Céline Baral, Fanny Dissard, Lionel Degermann, Anne-Sophie Canivet, Alexandra Loaiza Freire, and Hugo Rain. We express our gratitude to the TotalEnergies carbonate team for discussions and support, and especially to Emmanuel Caroli and Vincent Marlot for their reviews. Philippe Rabiller provided valuable technical guidance and proofread an earlier version of the paper. The TotalEnergies Offshore Angola team, particularly Nikki Morris and Andrew Hoover, was an enthusiastic partner and early adopter of this method, and instrumental in the release of data and approvals for publication. We are grateful to TotalEnergies colleagues Alan McInally, Jean-Yves Frouté, and Francisco Javier Sancho Jaquel for championing the carbonate rock-typing project. We also appreciate the reviewers for their efforts and valuable feedback toward manuscript completion. ADNOC is warmly acknowledged for providing rich datasets, engaging in technical discussions, and for allowing us to publish the results. A warm thank you goes to our dedicated technical and style proofreader, Beatriz Garcia-Fresca. Lastly, and most important, Jeroen Kenter wishes to express his deepest gratitude to his mentor in petrophysics and friend, the late, great Mark Skalinski, to whom this publication is dedicated.

Ousmane Mangane received a Ph.D. (2013) in CO₂ underground storage (University of Montpellier, France). He is a research engineer in “in-situ recovery” for uranium mining. He has extensive experience in petrophysics, well log analysis and in carbonate sedimentology. He focuses on reactive transport in porous and integrated petrophysics for carbonate reservoir characterization. He has developed workflows for carbonate rock typing.

Jeroen Kenter has a Ph.D. in geology (1990) and worked in academia and for Chevron ETC, Statoil ASA, ConocoPhillips and TotalEnergies (carbonate expert). Jeroen’s passion is reservoir characterization, rock typing workflows, outcrop modeling and carbonate prediction through time and space. He mentored eight Ph.D. students, (co)-authored more than 100 publications, holds two patents and retired in 2023.

Tarik Mecheri earned his M.Sc. in engineering (geophysics) from the University of Science and Technology Houari Boumediene in Algiers, Algeria in 1998. He began his career at SLB (formerly Schlumberger) as a geophysicist, focusing on reservoir characterization and geomodelling. Currently, he serves as a technical support consultant for TotalEnergies in the Modeling Simulation Production & Optimization division.

Jean Borgomano earned his Ph.D. in 1987 from Aix-Marseille University, France. He worked for Shell International in exploration–production (1988–2003). He worked as a carbonate expert at Total (2013–2015) and then returned to Aix-Marseille University where, as a professor, he is leading the Carbonate Reservoir Laboratory. His research focuses on the characterization and numerical modeling of carbonate reservoirs.

Christophe Nussbaumer received a Ph.D. in earth sciences from the University of Geneva (Switzerland) in 1997. He joined Total as carbonate sedimentologist and focussed on reservoir quality investigations of various reservoirs in the Middle East, Europe, Asia and Africa. He currently works as carbonate studies coordinator.