We investigated the basic geological and petrophysical properties of the multimodal pore systems in the Arab D limestone facies in Ghawar field, Saudi Arabia. The study used more than 500 mercury injection capillary pressure (MICP) data, which were type-curve matched using Thomeer Hyperbolas. The new MICP sample data were drawn from 10 cored wells that transect the Ghawar field from north to south and from a previous fieldwide study with 125 MICP samples. These 500 samples have a very rich statistical foundation in that they were selected using only random decimation within each of the facies from more than 3,500 core plugs all with assigned facies. In addition to MICP data and facies, a former, smaller sample set had both facies and Dunham texture codes. A new view of these pore systems emerged, that is built upon the intrinsic, fundamental and separate maximum pore-throat diameter modal elements named porositons. Porositons are stable, recurring and intrinsic modes in the maximum pore-throat diameter of the carbonate pore systems. Analytical results derived from the MICP data showed that the pore systems of the Arab D limestones can be classified based on porositons. The benefits of this new classification are demonstrated by considering in detail the relationships to geological facies, well-log responses, permeability modeling and simple nuclear magnetic resonance (NMR) well-log response. By analogy to the decoding of the Egyptian hieroglyphics using the Rosetta Stone, the use of porositons enables strong connections to be made between the geological facies, petrophysical and reservoir-flow properties of these complex carbonate rocks. The relationships between the new pore systems categories and the facies were thoroughly tested using north-south field trends.


Understanding pore systems is essential for modeling the properties and performance of hydrocarbon reservoirs. Pore systems provide the primary control on hydrocarbon distribution during reservoir charging. They control the interaction of the rock with fluids through wettability modification, and fundamentally control the hydrocarbon storage and recovery through the properties of porosity, permeability, relative permeability and microscopic-displacement efficiency. Pore systems can be examined in detail with core material and perhaps linked to well-log responses, but to fill the inter-well volumes of reservoir models with pore system properties, they must also be linked to predictive geological parameters. Thompson et al., (1987) in a portentous review of the pore geometrical problems of sedimentary rocks, noted “Prediction of rock properties, such as the transport properties of fluids in the pore space, and the elastic properties of the grain space, requires a set of statistics that embody the relevant physics”; and “More generally, the statistical description of pore geometry awaits definition of relevant statistics. This approach could ultimately tie the geology of rock formation to their reservoir properties, a tie with important consequences for oil exploration and production.”

Multiple and complex pore systems are commonly encountered in carbonate reservoirs. Studies of the Arab D limestones in Ghawar field, Saudi Arabia (Figure 1), have demonstrated the presence and approximate volumes of multiple-porosity types as two-dimensional petrographic information (Cantrell and Hagerty, 1999, 2003; Hagerty and Cantrell, 1990 unpublished report). In this paper, we extract three-dimensional pore geometrical statistics that embody the relevant physics using results from Thomeer Hyperbola (Thomeer, 1960) analysis of mercury injection capillary pressure data (MICP; see Appendix 1 for a review of Thomeer analysis and comparison to similar techniques, e.g. Leverett J Function), conventional core analysis data (porosity and permeability) and qualitative analysis of nuclear magnetic resonance (NMR).

The aspects of the pore system we extract are:

  • (1) Pore-subsystem volume, which has a counterpart in the size-differentiated, point-count data;

  • (2) Pore-subsystem geometrical factor, which has a counterpart in the sorting coefficient or size-distribution width of the size-differentiated point-count data;

  • (3) Maximum diameter of the pore-throat that controls the pore subsystem. This is only determined by Thomeer analysis of MICP data, and only has a weakly defined petrographic counterpart. This diameter value, however, plays a fundamental role in a new pore-system classification scheme, allows the pore-subsystems to be related to height above free-water level in the reservoir (effectiveness) and controls the permeability.

The results that emerge from our analysis intersect several petroleum reservoir disciplines and, more explicitly, allow us to relate their subsurface languages. Our approach is similar to the translation-triangle applied to the decoding of Egyptian hieroglyphs using the Rosetta Stone. This paper focuses on the subsurface languages involving Arab D static reservoir properties: depositional facies, well-log responses and the pore systems. Other translations involve the dynamic properties of permeability, relative permeability and microscopic-displacement efficiency and speak to a petrophysical and reservoir engineering audience. These have been published in the engineering literature (Clerke, 2007).

The paper begins with some general information about the Ghawar field and previous geological studies of Arab D facies and carbonate microporosity. Our analysis is focused on limestones and does not include dolomitic facies, nor do we consider the important role of fracture porosity and permeability. We then introduce the statistics of the maximum pore-throat diameter to characterize the multimodal Arab D limestone facies. We discuss density-neutron measurements, porosity and permeability as related to multimodal limestone pore systems. Using only basics of NMR signal analysis, we demonstrate an alignment of the behavior of maximum pore-throat diameter and NMR-detected pore-body diameters. We explore relationships between facies and porositons, and the calculation of permeability. Finally, we conclude that these new investigations create many new paths for improvements in the evaluation of complex carbonate reservoirs.



The super-giant Ghawar oil field of Saudi Arabia is a N-trending anticline that is 230 km long and approximately 30 km wide (Figure 1). The main reservoir is the carbonate part of the D Member of the Upper Jurassic Arab Formation (Powers, 1968). The Arab-D Member is the oldest of four carbonate – evaporite cycles (from bottom to top, Arab-D, Arab-C, Arab-B, and Arab-A with the overlying Hith Anhydrite), each stratigraphically comprising a lower carbonate unit and an upper evaporite unit of which the evaporite intervals dominated by anhydrite. The D Member carbonate consists of several scales of shallowing-upward cycles dominated by burrowed mudstones and wackestones in its lower part and transitioning to skeletal packstones and grainstones and, ultimately, ooid grainstones in the upper part.

Many articles have been written on the super-giant Ghawar field, its geology and the performance of the Arab D Reservoir (e.g. Powers, 1968; Mitchell et al., 1988; Al-Husseini, 1997; Stenger et al., 2003; Figure 2 after Lindsay et al., 2006), but few have integrated its pore systems with the geology and reservoir performance. Moreover, although thousands of well logs, production surveys, cores and extensive production histories are available as data to characterize the Arab D Reservoir performance, the fundamental understanding of porosity and permeability and prediction of permeability remains a challenge.

With reference to the Ghawar field, Cantrell and Hagerty (1999, 2003) documented observations on the common presence of microporosity in the Arab Formation. They concluded: “Microporosity occurs throughout the Arab Formation of Saudi Arabia, and affects the log response, fluid-flow properties and ultimate recovery of hydrocarbons in these reservoirs.” Their qualitative microporosity observations from petrographic data, were supported by displays of the pore-throat histograms from a sample subset on which MICP data were available (Figure 3, from Cantrell and Hagerty, 1999).

Early work on the presence of microporosity in carbonates by Herling (1968), contrasted the difference between carbonate and clastic sedimentary rocks. He noted, “The good correlation between specific surfaces calculated from grain-size distribution and measured with BET method (adsorption isotherm method of Brunauer et al., 1938) for quartz sands and powdered quartz is no longer valid for calcareous sands and silts.” Pittman (1971) discussed one type of microporosity in carbonates as resulting from boring and perforating actions performed by blue-green algae into carbonate grains. A similar observation was made by Bathurst (1966) for skeletal sand grains of the Bimini lagoon.

Cantrell and Hagerty (2003) found four types of microporosity in the Arab D using an operational definition of microporosity as pores approximately 10 microns in diameter or less. Their four types of microporosity are: (1) microporous grains, (2) microporous matrix, (3) microporous fibrous to bladed cements, and (4) microporous equant cements, with microporous grains being the most volumetrically significant microporosity type. These authors stated, “Scanning Electron Microscope examination of pore casts and fractured rock surfaces reveals that a variety of skeletal and non-skeletal grain types are microporous. The microporosity–forming process transforms different grain types into grains that are similar with respect to their internal fabrics. Microporosity thus consists of a network of highly interconnected, uniform-sized straight tubular to laminar pore-throats that intersect with less elongate, more equant pores.”

Mitchell et al. (1988) characterized the Arab D carbonates in terms of five limestone facies plus Dolomite: Skeletal-Oolitic, Cladocoropsis, Stromatoporoid–Red Algae–Coral, Bivalve-Coated Grain-Intraclast and Micrite. These facies represent a classification that utilizes fauna, assemblages of Dunham (1962) textures, major and minor grain types, sedimentary structures, pore types and diagenetic modification. A key result that emerges is that these facies subdivide the MICP data set by pore-system properties more clearly than other available facies descriptors. This establishes an important and previously nonexistent link between the depositional geology-facies and the pore systems, and hence the static and dynamic properties of the reservoir.

Clerke (2004) modified the limestone facies of Mitchell et al. (1988) and arrived at six facies (Figure 2, Appendix 2). In particular, he divided the skeletal-oolitic facies into those above and below an important sequence boundary identified in core descriptions by C.R. Handford (T. Keith, personal communication, 2001). His modification also separates-out the burrowed-micritic facies. In the present study, the Clerke (2004) classification is used unless otherwise indicated.


Thompson et al. (1987) discussed the many studies of pore networks “composed of pipes of widely varying sizes, which are distributed randomly along the links of the network”; and noted, “there are no experimental data to contradict the assumption of random distribution of pores.” The data we have collected and present is likely the first and the most comprehensive data to show a deeper and non-random structure in the pore network parameters of the Arab D limestone.

Two data sets were used in the present study. The first consisted of 125 mercury injection capillary pressure (MICP) samples characterized by the textures of Dunham (1962) and the facies of Mitchell et al. (1988). This data set was compiled by Hagerty and Cantrell (1990, unpublished report) and quantitatively analyzed (Clerke, 2003, 2004; Clerke and Martin, 2004). It is dominated by packstone (Dunham, 1962), but otherwise contains a fairly uniform selection of other Dunham (1962) textures. The collection is also fairly uniform by facies except for undersampling of the Cladocoropsis facies.

Thomeer Hyperbolas (Thomeer, 1960) were fitted to this first MICP data set using computerized spreadsheets (Appendix 1; Clerke and Martin, 2004). This analysis first yields the number of Thomeer Hyperbolas – used here as equivalent to ‘pore system’ – required to fit the MICP data from each sample. We call this integer the pore system modality. Thirty-five percent of the samples required a single hyperbola (Figure 4a), 62% required two Thomeer Hyperbolas (Figure 4b) and 3% required three hyperbolas. To limit trivial occurrences of multiple pore-systems in one sample, we added to the Thomeer MICP-fitting process, the requirement that a volume of at least one unit of porosity be present for a significant second or a third pore system.

Table 1 shows the occurrence of the pore-system modality (number of Thomeer Hyperbolas required for each sample) for the Dunham (1962) textures. Key results are that mudstone is dominantly monomodal, while grainstones and mudlean packstones are dominantly bimodal; packstone and wackestone showed no definitive pore system modality. Table 2 shows the occurrence of the pore-system modality for the Mitchell et al. (1988) facies. Stromatoporoid-Red Algae-Coral and Cladocoropsis are almost always bimodal. Note also the presence of 60–40 or 70–30 splits for Skeletal-Oolitic, Micrite, Dolomite and Bivalve-Coated Grain-Intraclast. Table 3 shows the pore-system modality when the Skeletal-Oolitic facies are distinguished into above-and-below the sequence boundary (Figure 2; Clerke, 2004) and the value of this split will be evident later when microporosity types are defined.

The second data set (Rosetta Stone Data) used 10 cored and described wells with full conventional well-log suites in a NS-transect across the Ghawar field (Figure 1). All of the core plugs (approximately 3,500) from these 10 wells were cataloged by the facies of Mitchell et al. (1988). Then seven to nine samples from each facies per well were selected at random to give 90 samples per facies, thus resulting in a uniform sample density for each facies. MICP experiments and analysis was conducted on these samples and the resulting data set was combined with the 125 samples of Hagerty and Cantrell (1990, unpublished report) as listed in Appendix 3.

The Thomeer Hyperbola analysis was then applied to the combined MICP data (Appendices 1 and 3) using computerized spreadsheets (Clerke and Martin, 2004). Thomeer parameter, Pd, when plotted on a logarithmic axis [Log(Pd) domain] showed the most significant variation, and it was considered the controlling parameter for classification. We observed four distinct and separable modes (Figure 5) and fitted each with a Gaussian distribution, termed a Porositon.

Pore-size ranges and microporosity have been assigned many names and size scales (Choquette and Pray, 1970; Pittman, 1971; Cantrell and Hagerty, 1999). In contrast, Figure 5 shows that maximum pore-throat sizes of the pore systems have a few natural and distinct modes (four) and hence self-define its size-based classification. The Thomeer Hyperbola curve-matching tool is the only one which allows the microporosity volume and its largest controlling pore-throat to be quantified separately from other pore volumes.


The four modes in the distribution of the maximum pore-throat diameters (Figure 5) correspond to one macroporosity (M Porositon) and three microporosity types (Type 1, 2 and 3 Porositons). The position of the four porositons is a property of the porous medium and not determined by ad hoc criteria. The term “microporosity” is used here in the sense of Swanson, (1985). “Micropores in reservoir rocks are defined as pores whose dimensions are significantly smaller than those contributing to the rock’s permeability”. Later sections of this paper and other publications (Clerke, 2007; Buiting and Clerke, in preparation) support this naming convention by demonstrating the lack of contribution of the micropores to the measureable permeability when the M Porositon occurs, i.e. macroporosity.

Results from the fitting of the data to Gaussian distributions are shown in Table 4. The distribution parameters for the four porositons are: mean, width (standard deviation) and the best Log(Pd) cutoff parameter separating the distributions. Notice from Table 4 that the bulk of the porosity on the average is carried in the M Porositon (17.1 pu) followed by the Type 1 Porositon (5.57 pu) and finally by Types 2 and 3 (2.22 pu) Porositons. The computed Thomeer permeability for the mean values is shown (see Appendix 1 for permeability calculation). The mean Thomeer pore geometrical factor, G, ranges from 0.51 to 0.13, with the highest values (widest distribution and poorest sorting) in the M Porositon (0.51 + 0.19). The characterization of carbonate “heterogeneity” is captured by the width of the M Porositon and by the amount of multimodality present in the pore systems.

Core Plug Porosity and Permeability, and Multimodality

In core plugs the pore systems of the Arab D limestone matrix are commonly bimodal and composed of a single instance of macroporosity and some amount and type of microporosity from one of three possible types (Clerke, 2004, 2007; Clerke and Mueller, 2006; Buiting, 2007; Buiting and Clerke, in preparation). We examined the trend of total plug porosity in the limestone samples as compared to the pore system modality. Figure 6 shows the modality of the maximum pore-throat diameter data with a red curve representing the rolling-window average of the data versus porosity. The data indicate that bimodality is common when the porosity exceeds 7 pu. Trimodal pore systems start to occur, along with bimodals, in the porosity range from 12 to 24 pu, and above 24 pu the systems are mostly bimodal. In Figure 7, the maximum pore-throat diameter modality is compared to core plug permeability data using the rolling average (red line). Two distinct ranges occur: (1) from 0.01 to 10 mD, which is both monomodal and bimodal in nearly equal proportions; and (2) above 10 mD, bimodality dominates with occasional trimodality.

Density-Neutron-Derived Porosity and Multimodality

In Figure 8, the plug maximum pore-throat diameter modality versus porosity data of Figure 6 is superimposed on the conventional density-neutron well-log crossplot using multi-well Arab D log data. The combination of density-neutron logs in a known limestone matrix might show some subtle evidence of maximum pore-throat diameter modality. The deviation of the data trend from the limestone line, as porosity increases, may be related to the increasing presence of dual-porosity systems through unknown effects such as flushing phenomena and/or flushed-zone, residual-oil saturation. Porosity, determined from these two porosity tools (density-neutron), appears to only weakly indicate the modality of the pore system. Hence, we observe that these “conventional” petrophysical data sets at most give minimal indication of multimodality. The characterization of multimodal pore systems in reservoir limestones will require new methods for full evaluation by well logs. Re-interpretation of existing (density-neutron) well-log suites even with core data and powerful reprocessing (e.g. neural nets and self-organizing maps) likely cannot begin to address this aspect of the evaluation.

Porositon Combinations

The porositon combination naming scheme for multimodal pore systems follows the MICP access order, i.e. it follows the order of increasing pressure or decreasing largest controlling pore-throat, i.e. M, 1, 2, 3. Hence a MICP trimodal pore systems can be described by the increasing Pd series: M–1–2, but not 2–M–1.

The major porositon combinations still exhibit order when viewed on the traditional porosity-permeability crossplot (Figure 9) shown here with an overlay of constant Pd using equation A3. Hard trends and clear clusters are not present but there are still definite patterns. The M–1 bimodals with their various proportions of macro- and microporosity occupy a large area in the upper right corner of the plot. Moving to the left and reducing porosity but retaining high permeability, we find the M–1–2 and M–2 pore systems with the additional presence of micrite associated micropores (2, 3). Trending down to the left from the M–1 cluster the Type 1 monomodals are intersected, then Type 1–2 and 1–3 bimodals and then at very low porosities and permeabilities, Type 2 and 3 monomodals associated with micrite.

The results from Table 4 and the porositon combinations were used to build average porositon combination capillary pressure curves with uncertainties (Clerke, 2004). The measured capillary pressure curve data have been disassembled using the superposition capability of the Thomeer hyperbolas (see Appendix 1) followed by the porositon classification. The average capillary pressure curve for each porositon combination can similarly be reassembled using the averages and standard deviations of each porositon’s Thomeer parameters and a forward superposition of the Thomeer hyperbolas. Figures 10 and 11 demonstrate the forward-modeled mercury capillary pressure curves for the major porositon combinations.

Facies and Porositon Combinations

In this section we examine the Arab D limestone pore systems and facies in terms of combinations of the four basic porositons. Two common combinations of macroporosity and the microporosity types are shown in pore-throat diameter histogram view (Figure 12). In the 125 sample data set of Hagerty and Cantrell (1990, unpublished report), we observed that the sequence boundary represents a distinct facies break that is associated with the microporosity Types 1 and 2 (Table 5). The Skeletal Oolitic facies above-and-below the sequence boundary have completely different microporosity types (Figure 3). Above the sequence boundary, Skeletal Oolitic, Cladocoropsis and Stromatoporoid-Red Algae-Coral facies share the common occurrence of the Type 1 Porositon. The strong correlation of the microporosity types against the facies (Table 5) is in stark contrast with the weak correlation with Dunham (1962) textures (Figure 13). It is seen that grainstones can either contain Type 1 or 2 microporosity as can Mud-lean Packstone and Packstone. Wackestones contain only Type 2 microporosity. Thus facies descriptors (Clerke, 2004) contain the most petrophysical pore system information.

Figure 14 shows the nine significant porositon combinations with the occupancy frequency coded by the facies. These values are obtained after using a noise threshold cutoff of 5.2%, e.g. M, M–1, M–1–2, M–2, 1, 1–2, 1–3, 2, 3. The project’s experimental design included multiple independent facies assignment cross-checks for all core plugs and thereby generated an internal facies error threshold of 5.2%.

Three of the facies tend to occur in the upper part of the reservoir interval (Skeletal Oolitic above-the-sequence boundary, Cladocoropsis and Stromatoporoid–Red Algae-Coral) and for these, only four porositon combinations are required with M–1 being dominant (Figure 14). Note that the monomodal, M, occurs only in the Skeletal Oolitic above-the-sequence boundary, and that the presence of Type 2 microporosity (Micrite associated) in the M–1–2 and M–2 is associated only with Stromatoporoid–Red Algae-Coral. An image of an M–1–2 Stromatoporoid–Red Algae-Coral trimodal core plug (Ahr et al., 2005) is shown (Figure 15) along with its MICP data and Thomeer fit, and the corresponding pore-throat presentation of the data and its fit.

For the facies occurring below-the-sequence boundary, a much more complicated situation occurs. The Skeletal Oolitic and Bivalve-Coated Grains-Intraclasts facies require five and six, respectively, of the nine porositon combinations (Figure 14). Similar in terms of their M–1 content, Bivalve-Coated Grains-Intraclasts contains 1–3 and 2, which do not occur in the skeletal oolitic facies, but which contains 1–2 not found in Bivalve-Coated Grains-Intraclasts. The Micrite facies is clearly dominated by Type 2 and Type 3 Porositons; 45% and 40% respectively, with a 15% presence of Type 1 Porositon (Figure 14). The intraclast content of the Bivalve-Coated Grains-Intraclasts facies presumably accounts for this presence of micritic pore systems.

Ghawar Field Trends by Facies and Porositon

The Rosetta Stone study was designed to give both field level and well-level statistical results. At the well level, the sample set allows a study of the north to south variation in the maximum pore-throat diameters by facies and by variation within the range of maximum pore-throat diameters defined in each porositon.

In Figure 16, the mean of the logarithm of the entry pressure Log(Pd) for each of the 10 Ghawar wells is shown for each of the facies. The data show a steady decrease in the Log(Pd) for Micrite microporosity from north to south, or a coarsening of the maximum pore-throat diameter, though all of these pore-throats are small. In contrast, the Log(Pd) steadily increases, fining from north to south for the other five facies, which are controlled by the M Porositon. Essentially, the largest pore-throats get smaller to the south. The picture is considerably different when viewed from the porositon perspective (Figure 17). The NS porositon trends are essentially flat except for the M Porositon, where Log(Pd) again shows an increase to the south. Hence the Type 1 to 3 Porositons appear to be uniform over the whole field.

When superimposing the Micrite facies trend over the porositon trends (Figure 18), we can reconcile these observations. Micrite is composed of Porositons Types 1 to 3, but dominated by Types 2 and 3. Hence the NS-trend for Micrite can only be explained by its enrichment with Type 2 and/or Type 1 compared to Type 3 porositons to the South. This indicates that the Micrite facies, being a varying mixture of the three microporositons cannot be used as a field-wide uniform petrophysical calibrator for well-log normalization. However, the Porositons 1, 2 or 3 can be used for fieldwide calibration and normalization.


Thomeer (1983) found an empirical equation for air permeability based on the three Thomeer parameters (Appendix 1, equation A3). We used his equation to compute permeability for all the samples using only parameters from the first pore system, hence neglecting the smaller pore-throat modes. In a sample comprised of multiple Thomeer Hyperbolas (pore systems) written in decreasing order of maximum pore-throat diameter, e.g. M–1–2, the first pore system is M. We compared the calculated and measured values for the monomodal (require one Thomeer Hyperbola) samples and the bimodal (require two Thomeer Hyperbolas) samples (Figure 19). The good match implies that the permeability of the Arab D limestones is dominated by the properties of the first, commonly M porositon. Detailed comparisons between monomodal and bimodal sample data showed no detectable permeability contribution that could be attributed to the neglected smaller modes. This observation starts to explain the poor results that would arise from a conventional total porosity-permeability approach (Delfiner, 2007) to the data as shown in Figure 9 where at 25 pu, the permeability ranges from 10 to 3,000 mD.

Having isolated the M Porositon as controlling the measurable permeability, we investigated the explicit dependence of the permeability on each of the three M Thomeer parameters (Clerke, 2007). The major control on permeability (Figure 20) is the Thomeer parameter, Pd,f, the diameter of the largest pore-throat in the first (ordered by decreasing maximum pore-throat diameter) pore system, which is found by the equation between capillary pressure and capillary diameter:

Pc = 0.58 × [ (σcosθ)/d]

Using values for the Mercury-Air experiment interfacial tension and contact angle (σcosθ), the mercury air pressure in psi; the diameter of the maximum pore-throat (capillary diameter) in microns is:

d throat,max(microns)=214/Pd,f

The first displacement pressure is directly related to the maximum pore-throat diameter, which is the maximum pore-throat diameter of the first porositon when ranked by decreasing maximum pore-throat diameter.

The correlation R2 between permeability and Pd,f is 0.65. Of the three Thomeer parameters, the most important for permeability is clearly Pd,f, or the maximum pore-throat diameter or the largest pore-throat diameter of the first (usually the Macro) pore system. The next best correlation is found between permeability and the total sample porosity, i.e. R2 = 0.55, but other Thomeer parameters could also be used with only slight decrease of the correlation. We choose to continue with the total sample porosity because of the ease with which it is determined from well logs.

Using these results, we propose a new, simple two-term permeability model, which has potential for well-site implementation using properly processed well-log data (Figure 21 and 22):

Log (Measured Permeability) =[a] + [b] × [Log (maximum pore-throat diameter)] + [c] × [Porosity (%)] = -1.544, b = 1.206, c = 0.0727

The measured-versus-predicted plot for the proposed two-term model has a correlation coefficient R2 of 89%. The measured-versus-predicted plot for the proposed two-term model shows excellent results over seven orders of magnitude in permeability (Figure 22). This model has many similarities to the approach of Lucia (1995) who also proposed a pore space-permeability classification based on two variables: the particle size and the interparticle porosity for non-vuggy rocks. This approach has been further developed to include the behavior of the imbibition oil relative permeability in multimodal M–1 pore systems (Clerke, 2007).


After correcting NMR log data for reservoir and borehole fluid and fluid-surface interaction effects, petrophysicsts can derive some information about the pore system or its attributes (e.g. permeability or saturation). Here we deploy a different strategy. We use the knowledge gained from studying the pore systems and pore-throats from our MICP studies of the Arab D limestone and investigate whether they are present in the NMR log data. Although we have not yet performed a rigorous correction for every fluid effect, our initial investigation focused on the conformance of the two data sets at every point of contact. This is supported by information that the Arab D Reservoir apparently has only minimal difference between the NMR properties of the two reservoir fluids at reservoir conditions (R. Akkurt, personal communication, 2007).

The pore system modality detected within the MICP data records the distinct occurrence of multiple-pore systems, each with a range of possibilities for the maximum-pore-throat diameter (see Figure 5b). The use of the words “pore bodies” and “pore-throats” tend to give rise to the concept that the pores themselves are constructed from discrete objects such as connecting tubes and spheres. This is often a useful approximation. In reality, a pore-throat is a critical pore body section whose radius and position limits access to the remainder of the pore body or bodies. The position and size of these pore-throats is governed by many factors such as the arrangement of the grains, grain sizes, sorting, angularity, configuration and cementation.

If pore-throats show multimodality in the MICP data, it is an attractive leap to conjecture that there may also be pore-body multimodality, termed “porobodons”. That is, the modes in the core plug MICP maximum pore-throat diameter spectrum might have counterparts in the NMR pore-body spectrum. The pore-throat and pore-body data are shown in Figure 23a for an M–1 pore system core plug where laboratory NMR measurements have been performed at 100% water saturation at reservoir temperature after which MICP data were acquired from the core plug.

At reservoir conditions and with reservoir fluids, the pore-body modality from an NMR well log in the Arab-D is displayed versus total porosity in Figure 23b. A pore-body modality indicator is easily computed by recognizing that it is equivalent to half the number of inflection points in the NMR signal. Alternatively, a threshold crossing-counter technique can also be applied. The NMR data from even one well is many times more abundant than the more than 500 MICP samples and it also represents a completely different scale of investigation. Yet the two plots (Figures 6 and 23b) are extremely similar, which tends to support the attractive ‘porositon-porobodon’ conjecture.

To further examine the “porositon-porobodon” conjecture, we anticipated its form by computing the histogram of the largest controlling pore-throats (Figure 5a) weighted with the average porosity of the histogram bin (Figure 24). Figure 24 shows a first estimate of what an NMR spectrum from the Arab D limestones might resemble if the porositons and porobodons do have the conjectured correspondence.

The porosity-weighted, pore-throat histogram is a basic estimate because the NMR device responds to the amount of porosity governed by a characteristic hydrogen spin-decay time rather than a pore-throat diameter bin. In Figure 24, we see that the M Porositon (to the left of 4.6 microns pore-throats) is heavily weighted by its abundant pore volume. The common Type 1 microporosity (from 4.6 microns to 0.35 microns pore-throats) is weighted by a smaller bin porosity; porosity Types 2 and 3 (less than 0.35 micron pore-throats) are now even lower. An NMR pore-body spectrum might look similar to this if the correspondence to porositons is correct. The x-axis values in terms of milliseconds of NMR hydrogen spin-decay times can only be estimated. The very large pore-throats (and associated pore bodies) will have long, many-second decay times. The Type 2 and Type 3 microporosity will have very fast NMR decay times (few milliseconds). For comparison, an actual NMR pore body spectrum for the Arab D is shown in Figure 25. There is indeed a distinct similarity between what we have estimated (Figure 24) and what is obtained (Figure 25). At each of these three points of contact between the two data sets; single core plug, modality versus porosity, porositon modeled spectral shape, the data are very similar.

The evidence is encouraging that the maximum pore-throat modality seen in core plug MICP is also present in the much larger rock volumes investigated by the NMR well logs. Even after upscaling, the porositons may have equivalent porobodons. This is a fruitful path for future investigations.


Mercury injection capillary pressure (MICP) data from the Arab D limestones in Ghawar field (Appendix 3) were analyzed using Thomeer Hyperbolas (Appendix 1). The analysis identified four discrete modes of maximum pore-throat diameter values termed porositons; the largest of these is macroporosity (M Porositon) followed by three types of microporosity (Porositons 1, 2 and 3). Modality is a new and important indicator for carbonate reservoir characterization because it demonstrates that the presumably continuous maximum pore-throat diameter axis can be well-represented using only four discrete modes.

The pore systems of the Arab D limestone (Appendix 2) have been described as combinations of these porositons, with 70% of the samples showing bimodal or trimodal behavior. Three of the six facies have nearly identical porositon combinations, leaving four “porositon combination effective” facies. The best reservoir quality shows a common M–1 bimodal pore system, with the M Porositon carrying the measurable permeability. The M–1 matrix pore system acts as a dual porosity – single measureable permeability system.

The number of modes in each of the facies and the type of microporosity are new and important geological and petrophysical attributes that are not accessible with conventional well logs but accessible through MICP analysis. A fruitful approach to NMR well logs processing is to extract pore-body modality (porobodons) that are related to the porositons. The conjectured porositon-porobodon correspondence suggests that a very similar behavior in terms of mode number, position and occupancy is present in the NMR signal, and early evidence is given to support this. This paper suggests that NMR data might be used to support facies identification.

Additionally, the dominant presence of multimodal pore systems in the best reservoir intervals of these carbonates requires that the static and dynamic property models be fully generalized for multimodal pore systems. We have demonstrated one generalization for permeability; showing that only attributes of the M (Macro) Porositon are significant. The generalization to imbibition oil relative permeability in multimodal M–1 pore systems has been published elsewhere (Clerke, 2007).

Using the Rosetta Stone analogy, this study was designed and executed as a three language decoding exercise, which enabled a strong connection to be made between data in the domains: geological facies, pore systems and reservoir-flow properties. The overall workflow and analysis processes as demonstrated for the Arab D limestones can most likely be applied to other multimodal carbonate pore systems and reservoirs when the appropriate data are collected.


Thomeer (1960) developed a method for the analysis of mercury injection capillary pressure (MICP) data, which was used primarily within Shell (Thomeer, 1983; Smith, 1992; Hawkins et al., 1993). He observed that the data from the MICP experiment, for simple rock types, could be represented by a hyperbola when plotted on Log-Log graph paper. The data are: (1) volume of mercury injected, and (2) applied pressure between the wetting (air) and non-wetting phases (mercury) as mercury intrudes into the pore space fraction. For the rock sample, the bulk volume is known, so that the volume of mercury injected can be re-expressed as a fraction of the total sample bulk volume, Bv. To fit the hyperbola to the data, the value of two asymptotes, B, Pd, are required. The Thomeer Hyperbola is shown in Figure A1, and can be expressed as:

Log ( Bv/B)×Log(Pc/Pd)=K, where K is  a hyperbola shape factor

The asymptotes are: B is the percent bulk volume occupied by mercury at infinite applied pressure and Pd, the displacement pressure required to first intrude mercury into the largest pore-throat. Thomeer chose to express constant K = Log [exp (- G)] such that Equation A1 becomes:

Bv/B=exp [G/Log(Pc/Pd)]

where G is the Pore Geometrical Factor and determines the shape of the hyperbola (Figures A1 and A2). In practice, Bv and Pc data from the MICP experiment are fit by equation (A2) to determine Pd, B and G, for individual samples.

Thomeer (1983), using a weighted regression on data from 279 rock samples, found a relationship between the three parameters and air permeability (Ka):


This equation reproduces the measured permeability to within a multiplicative uncertainty of 1.82 (Figure A3 and Table A1). With this result, Thomeer demonstrated why permeability-porosity crossplots fail, especially in carbonates, because important details about the pore network are missing in the relationship.

Spreadsheet Implementation

Unlike the plastic overlays that were originally used to implement this technique, today the Thomeer analysis method can be implemented in an ExcelTM spreadsheet, available from many of the MICP data providers (Clerke and Martin, 2004). Hyperbolas are fit using the Solver function to minimize the error signal derived from the difference of the actual versus predicted capillary pressure curve. Both the data and its derivative can be used to fit the Thomeer Hyperbola, analogous with type curve matching in Pressure Transient Analysis. Superposition of Thomeer Hyperbolas for multiple pore systems is carried out by parsing the pore system in the pressure domain. All of this is implemented in a highly interactive way. The Thomeer permeability is continuously generated for comparison to the measured permeability. The result for the 125 MICP Arab D Ghawar Hagerty and Cantrell (1990, unpublished report) samples is shown in Figure A4.

A basic result from the Thomeer analysis of the MICP data is a record of the number of Thomeer Hyperbolas (used here as equivalent to ‘pore system’) required to match the data. We call this integer the “pore system modality”. The minimum would be one and the maximum encountered in our study was three. To limit trivial occurrences of pore system multimodality, we add to the Thomeer MICP fitting process the requirement that a volume of at least one unit of porosity be present for a significant second or even third pore system.

Closure Correction

In the process of obtaining a core-plug or rock-chip sample, an artificial boundary is placed into the pore geometry – the cylindrical or outer-plug surface. The plug saw cannot cleave grains of the sample and so leaves a rough outer surface where grains may have been removed. If the sample contains very large pores, these can be dissected by the sample/plug surface in ways not representative of the formation. These surface irregularities and pore-body dissections create an initial storage of mercury that must be filled before the normal pore-system geometry is encountered by the mercury. This volume artifact requires a closure correction, and is one of the most troublesome issues in the analysis of MICP data, particularly since a small relative volume change in the large pore-throats strongly affects the permeability calculation. This effect and its correction are illustrated for sample 30 (Hagerty and Cantrell 1990, unpublished report), which also exhibits a dual-pore system at about 1,000 psi (Figure A5).

An important feature of the interactive Thomeer spreadsheet is the ability to pick and revise the closure correction to iteratively solve for the Thomeer parameters while comparing the computed and measured permeability. The spreadsheet also incorporates the ability to select the closure correction and visualize the corrected data and the hyperbola in the derivative or Pore-Throat Size Distribution space. The selection of a closure correction volume (red bar) is illustrated for a bimodal sample on the Pore-Throat Diameter (derivative) plot in Figure A6. Data is in black, Thomeer hyperbola in red.

Application of Thomeer Analysis

The use of the Thomeer method for the Arab D Reservoir at Ghawar is focused on three issues: (1) facies-dependent saturation-height modeling; (2) free-water level determination; and (3) facies-dependent recovery estimation. The general process of defining petrophysical rock types based on Thomeer analysis of capillary pressure behavior is illustrated in Figure A7. Fitting of the MICP data is done for all samples and the resulting parameters are statistically analyzed, PRT’s are defined and then forward modeled.

The result of the process shown in Figure A7 is displayed in Figure A8, which shows five clusters (PRT’s) in the Thomeer parameter space (not the Arab D) and the forward model of each cluster’s average Thomeer parameters. This application is a specific branch of the broader range of applications of MICP and capillary pressure data; namely: (1) facies-dependent saturation-height modeling; (2) conductive minerals; (3) fresh or varying formation water salinity; (4) very low porosity; (5) complex pore systems – fractured, vuggy, microporous; thin sands; (6) determination of depth to water contacts; (7) evaluation of downdip potential; and (8) predictions in the absence of or poor resistivity logs.

Properties of Thomeer Hyperbolas

A single pore system can be represented by one Thomeer Hyperbola and is characterized by three parameters: Pd, B and G, each of which has measurable errors to characterize uncertainty. This is in contrast to other methods (e.g. Leverett J Function, see next page) that are based on MICP data, but introduce additional uncertainties when scaling the data to measured permeability and porosity. The Thomeer parameters are intuitive: largest pore-throats, sorting of pore-throats, and total amount of porosity. They describe the pore system in a similar manner to grain system: i.e. largest grains to largest pore-throats, sorting of grains (e.g. Trask sorting), and sorting of pore-throats (G).

Petrophysical Rock Type (PRT) can be defined as an object in the Thomeer parameter space (Pd, B and G), i.e. a similar pore system. The parameters for PRTs can be operated on statistically (mean and standard deviation) to create average capillary pressure curves and their error bounds. Estimates of recoverable hydrocarbon can be guided by the knowledge of a complex composite pore system and hence by PRT.

Air permeability can be computed and predicted from the pore network parameters, Pd, B and G, to within a multiplicative uncertainty of 1.8, and this can be compared to a measured permeability. A Thomeer forward-modeled capillary pressure curve can be generated from insight into the generating parameters that can come from a variety of sources of rock data, from cores to cuttings and the Shell Rock Catalog (Thomas et al., 1995).

Thomeer Hyperbolas parameterize the process of filling the pore space with mercury (nonwetting phase) while displacing air (wetting phase), a drainage process. Suitable adjustments can be made to convert this behavior to other fluid systems in drainage (oil-water, gas-water).

Thomeer Hyperbolas can be combined or superposed to quantify complex pore systems. Superposition can be used to quantify porosity types and entry pressures, i.e. microporosity. The Thomeer function parameterization of MICP data has been shown to upscale (Buiting, 2007).

Square Root K/Ø Methods: Leverett J Function and Flow Zone Indicator (FZI) Technique

Whereas the Thomeer method has been used for many years in Shell, in other companies the Leverett J Function method (Leverett, 1941) and the related FZI technique (Amafuele et al., 1993) have been more extensively used. Here we discuss and compare these methods. The Leverett J Function characterization of a capillary pressure curve is defined by:

J(Sw) = (0.217 Pc/σ)×Κ/ø)

where Sw is water saturation, K is permeability, ø is porosity, Pc is the capillary pressure and σ the interfacial tension between the wetting and non-wetting fluids. The rock’s pore system is characterized by K/ø, and its square root is termed the pore system speed as it is the ratio of the flow parameter to the storage parameter. The Leverett method seeks to reduce the mathematical complexity of the complete behavior of many capillary functions by finding a scaling relationship among them. It is important to note that the use of this parameter to characterize the whole pore system in a complex pore system (mixture of macro, meso and micro pores) is incomplete, because the flow behavior of the pore system is strongly dominated by the large pores (low capillary pressure).

The FZI and Leverett J Function techniques are related in that they both use the √(K/ø) parameter, but in the former the height dependence is removed and replaced by a fixed value. The Reservoir Quality Index (RQI):

RQI = 0.0314 ×Κ/ø), which is J(Sw) when (0.217 ×Pc)/sigma=0.0.314

or when Pc = 4.5 psi (again low capillary pressure). The FZI method invokes the use of a constant RQI versus NPI (Normalized Porosity Index) slope.

This ratio, FZI, measures to what extent it is true that when we increase porosity by replacing matrix by pores, we do this by adding pore systems of the same speed √(K/ø), and thereby keep FZI constant. Here we make contact with the capillary bundle model of Purcell, because when pore systems of the same speed are added, it is equivalent to adding one or more tubes of the same type to the bundle.

NPI = ø/(1ø), and FZI =RQI/NPI

Square Root K/ø methods have several limitations. They cannot be superposed and do not completely characterize the high-pressure behavior, nor microporosity issues. Hence, recovery issues associated with multimodality are not well-represented. The Leverett J Function is dependent on two laboratory variables, k and ø; their individual uncertainties are compounded. It is not rock-texturally intuitive, does not provide a prediction of the permeability and hence does not provide a quality-control cross-check on the permeability data and measurement. Further, the clustering and inclusion of high-pressure behavior are a subjective visual process not leading to statistical measures. Finally, the approach does not lead to a method for determining the average capillary pressure curves and their uncertainty bands from laboratory data. Ekrann (1999) stated that “The true effective capillary pressure curve is not proportional to the rock Leverett J Function.”


DATA BASE (Excel Spreadsheet available from the author)


We thank the many colleagues who have contributed to this work from the Saudi Aramco Reservoir Characterization Department and the Reservoir Management Department, from Saudi Aramco management to the Saudi Aramco core warehouse. In particular, we acknowledge John Neasham of PoroTechnologies for excellent Mercury Injection Capillary Pressure Measurements and Dave Cantrell and Abdulkader Afifi for reviews and comments. We also thank the Ministry of Petroleum for permission to publish the paper and the careful and thoughtful work of the GeoArabia reviewers and editors and constructive suggestions from Gretchen Gillis. The authors also acknowledge Richard Koepnick for encouragement. The final design of the paper by Arnold Egdane is appreciated.


Edward A. (Ed) Clerke is a Geological Consultant in the Geological Modelling Division/Reservoir Characterization Department of Saudi Aramco. Previously, he held the positions: Head of Petrophysics – Petrophysical Advisor for Pennzoil, Senior Principal Petrophysicist with ARCO in Plano, Texas and Petrophysical Engineering and Research Positions with Shell Oil Co. USA. Ed received his PhD in Physics from the University of Maryland in 1982 and has published articles in Log Analyst, SPE Production Engineering, Physical Review, Physica and the Journal of Physical Chemistry and received the Best Paper Award at GEO 2004 and 2006 for work presented on Arab D Limestone Pore Systems and also received the 2006 Best Paper Award at the SPWLA Carbonate Permeability Topical Conference. Less recently, Ed received the 1993 Best Paper Award from the West Texas Geological Society for work in Permian Basin carbonates. He holds five patents, four for borehole televiewer technology.


Harry W. Mueller III has been a Carbonate Specialist with Saudi Aramco for seven years, most of that time in the EXPEC-ARC Department. Prior to that, he worked in a similar capacity with ExxonMobil for 26 years, much of that time working on problems in Middle-Eastern carbonate geology. He has a PhD in Geology with emphasis on carbonate stratigraphy and sedimentology from the University of Texas at Austin and a BSc in Geology from Birmingham-Southern College.


Eugene Craig Phillips is presently a Petrophysicist for Saudi Aramco located in Dhahran, Saudi Arabia working in Geologic Modeling. Since 1984 he has been primarily involved in core-log integration. He was formerly with Core Laboratories working in Technical Applications Development where he concentrated on core-log integration emphasizing core calibrated log analysis models. In 1994 Eugene joined Western Atlas in the Interpretation, Research and Advanced Development group where he developed NMR log and core applications.


Ramsin Y. Eyvazzadeh is a Petrophysical Specialist in the Petrophysical and Special Studies Unit, Reservoir Description Department in Saudi Aramco. Previously he has held several positions globally as Petrophysicist with Schlumberger and Numar Corporation. Ramsin has a BSc in Petroleum Engineering from the University of Southern California and a MS in Mechanical Engineering from the University of California at Berkeley.


David H. Jones is a Reservoir Engineering Consultant for Quantum Reservoir Impact. He previously spent 13 years with Saudi Aramco’s Reservoir Management Department. He led the multi-disciplinary GIANT project charged with developing and applying leading-edge technologies and methods to the determination of OIIP, ultimate oil recovery and the reservoir surveillance program for the Ghawar field - the largest oil field in the world. David is a recognized expert in reservoir management metrics and reservoir performance. His work in these areas led to the development of in-place systems at Saudi Aramco for tracking of all the companies’ operating reservoirs. Prior to joining Saudi Aramco in 1993, David worked fifteen years at Chevron and Gulf covering a wide range of domestic and international reservoir projects. David received a BSc in Petroleum Engineering from the University of Kansas. He has participated as an invited speaker to corporate and professional forums and is the author of numerous corporate technical documents.


Raghu Ramamoorthy joined Schlumberger Wireline Services in 1982 and served as a Logging Engineer in Egypt and Iran, and later as a Log Analyst in Oman. He received an MSc in Petroleum Engineering from the University of Texas at Austin in 1994. He then spent an extended stint as Research Scientist (Petrophysicist) at Schlumberger-Doll Research in Ridgefield, Conn., where he worked on reservoir characterization and carbonate petrophysics. In 1997 Raghu returned to the field and served as Principal Petrophysicist in Australasia, East Asia and Saudi Arabia. Raghu is currently Petrophysics Advisor in the Schlumberger Carbonate Regional Technology Center based in Abu Dhabi.


Ashok Srivastava received his Bachelor of Technology in 1988 from the Indian Institute of Technology (Mumbai). He started his career in Schlumberger as a Field Engineer in 1988. Currently he is Director of Curriculum for Petrophysics for the Data and Consulting arm of Schlumberger. Prior to this he has been a Principal Petrophysicist and has worked extensively in the Middle East area. He has various publications in formation evaluation and reservoir monitoring. Ashok has two patents, one in carbonate petrophysics and another one on mathematical techniques.