Stochastic modelling of spatial variability of petrophysical properties in parts of the Niger Delta Basin, southern Nigeria

Three-dimensional models of petrophysical properties were constructed using stochastic methods to reduce ambiguities associated with estimates for which data is limited to well locations alone. The aim of this study is to define accurate and efficient petrophysical property models that best characterize reservoirs in the Niger Delta Basin at well locations and predicting their spatial continuities elsewhere within the field. Seismic data and well log data were employed in this study. Petrophysical properties estimated for both reservoirs range between 0.15 and 0.35 for porosity, 0.27 and 0.30 for water saturation, and 0.10 and 0.25 for shale volume. Variogram modelling and calculations were performed to guide the distribution of petrophysical properties outside wells, hence, extending their spatial variability in all directions. Transformation of pillar grids of reservoir properties using sequential Gaussian simulation with collocated cokriging algorithm yielded equiprobable petrophysical models. Uncertainties in petrophysical property predictions were performed and visualized based on three realizations generated for each property. The results obtained show reliable approximations of the geological continuity of petrophysical property estimates over the entire geospace.


Introduction
Evaluation of petrophysical properties is a vital component in reservoir characterization in the oil and gas industries. The conventional petrophysical properties needed for hydrocarbon exploration, characterization, field development, evaluation and reservoir management are lithology, porosity, permeability, net-to-gross and fluid saturation (Emujakporue 2017). Detailed petrophysical property characterization is critical during the oil and gas exploration, production, management and surveillance stages. Parameters such as porosity, permeability and volume of sand and shale that are of interest to geoscientists and engineers are products of complex chemical and physical processes. It is significant to understand their spatial orientations and scales as well as the uncertainties that characterize these variables in space from one well to another in a particular field, for optimal hydrocarbon production. Permeability which is a dynamic property of the reservoir is difficult to characterize due to the irregular nature of pore structures occasioned by depositional and diagenetic alterations. The spatial distribution of rock permeability in heterogeneous reservoir systems is problematic since its determination has no direct solution (Fegh et al. 2013). Generally, two direct and reliable ways of determining permeability are the laboratory measurements and well test methods (Noah and Shazly 2014), which are expensive, limited to either the well bore or certain sections of the well bore and cannot be applied to all the wells drilled in the field. Alternatively, estimation of permeability from porosity-permeability cross-plot generated through linear correlation of porosity and permeability data is being used and has gained wide patronage in the oil and gas industries (Jennings Jr and Lucia 2001;Rezaee et al. 2006). However, this method can provide adequate result in sandstone reservoirs with little or no complexities and may not be relied upon when considering reservoirs with higher degree of heterogeneities due to episodes of diagenetic alterations (Edigbue et al. 2015;Ebong et al. 2019).
The application of indirect method, i.e. geostatistical methods in complex reservoir situations can provide information such as rock types, fluid contents and other physical processes within the formation. It involves the use of seismic sections, wireline logs and measurement while drilling and yields reliable results in a less expensive and faster way and can be performed over a large area (Dowell et al. 2006). Results generated from direct methods which are localized are used to estimate petrophysical properties elsewhere beyond wells via stochastic and fuzzy logic processes (Esmaeilzadeh et al. 2013). Geostatistical assisted modelling which is one of the many methods often applied to reservoir characterization (Bueno et al. 2011) has the ability to integrate several groups of information in the generation of suitable reservoir models that fit any given subsurface geologic condition of interest (Caers and Zhang 2002;Liu et al. 2004). Geostatistical principles also ensure that the geologic realities of the reservoirs are not lost in the course of the model building (Chen et al. 2011). The geologic model of the training image includes lithofacies assemblages from seismic sections, well log and production data as constraints and the petrophysical model which consist of the parameters of each facie (Caers 2002;Gonzalez and Reeves 2007). Akin to the deterministic approach, hard data points are usually conserved where it exists and have been interpreted and soft data where they are useful (Zhou et al. 2014;Wilson et al. 2011). In contrast to the deterministic approach, geostatistics provides several plausible outcomes (Zarei et al. 2011). When performed over a large area of the reservoir, results from such spatial investigations permit better characterization of the reservoir and effective recoverable reserves needed for adequate reservoir management can be estimated (Soleimani et al. 2017). In this study, the stochastic technique which is capable of generating more reliable numerical petrophysical models was used, due to the difficulty of predicting the spatial distribution of these properties deterministically (Ma et al. 2008).
This study which involves the estimation of petrophysical properties is aimed at estimating reservoir properties at well locations and predicting the spatial variability of these properties elsewhere within the field.

Stochastic reservoir modelling in brief
Petrophysical properties such as porosity, permeability and fluid saturation are the requirements for building 3D grid during flow simulation. These properties are limited to well locations and can pose a great deal of uncertainty in grid block property assignment. Building alternative numerical models or images of these petrophysical properties, taking into consideration the unknown aspects of the spatial distribution is generally referred to as stochastic reservoir modelling (Deutsch 1992). The traditional geostatistical approach to reservoir property modelling is through the sequential simulation of facies and/or petrophysical properties (Arpat 2005). The sequential simulations usually performed for stochastic reservoir modelling are the sequential Gaussian simulation which is parametric and the sequential indicator simulation, i.e. non-Gaussian and nonparametric approach (Xu et al. 1992). These two are performed using computer codes (Deutsch and Journel 1998). The sequential simulation path is such that it visits each node of the model and simulated values are drawn from conditional distribution of the values at the node given the neighbouring subsurface data and previously simulated values (Arpat 2005). Simulation techniques generate multiple realizations of unknown spatial distribution of observations, without losing the originality of the observed data. The differences between realizations provide quantitative measure of uncertainty (Vasquez 2014).
Since sequential simulation implies the reproduction of desired spatial properties through sequential use of conditional distributions, we can represent the primary variable of interest distributed over a field A as, z 1 (u) and the coordinate vector, u ∈ A . The primary variable in this work is the well log data from which porosity, permeability, fluid saturation, u, were derived on a three-dimensional scale. The integrals of z 1 (u) are associated with total pore volume. Randomly distributed values z 1 (u ), = 1, … , n 1 are available at well locations u ∈ A , and the seismic data z 2 (u � ), = 1, … , n 2 of related nature are available at much larger number n 2 locations, u � ∈ A . The purpose of this coupled approach is to produce several highly resolved maps of the distribution z 1 (u) over the field, A, making the most of both hard {z 1 (u ), = 1, … , n 1 } and soft {z 2 (u � ), = 1, … , n 2 } data. The term "hard data" is used to emphasize the fact that the modelling method should exactly reproduce the point data obtained from wells at its location. On the other hand, "soft data" which is the seismic section is used as constraint to guide the process and may not be reproduced exactly.
Data conditioning is significant especially when seismic data is utilized in estimating petrophysical properties. Hence, both hard and soft data must be conditioned before usage. Hard data conditioning may not only mean an exact reproduction of the original data but also requires the generation of adequate continuity around the defined area as dictated by the geological continuity model (Arpat 2005). For instance, if the geology of the reservoir shows evidence of thin, continuous horizontal layer of shale, a well data which indicates the shale value at a particular location should be honoured by generating such shale layer at neighbouring data location. Variogram-based simulation algorithms always exactly honour hard data values via the use of kriging, which is an exact interpolation method (Lima 2005). Soft data which refers to data sets obtained from indirect measurements, e.g. seismic sections, represents the "filtered" view of the subsurface heterogeneity (Arpat 2005). Since the filtered view is not the exact physical response, it is usually approximated by a mathematical model, which often times is referred to as "forward model". Soft data sets are usually integrated using variogram-based methods (Deutsch and Journel 1998).

Estimation algorithm
Kriging is a spatial interpolation technique commonly used in estimating varying petrophysical properties of reservoirs based on variograms (Matheron 1962). It is an estimation technique used to populate reservoir layers with continuous varying petrophysical values by means of interpolation in an unbiased pattern with minimum variance (Lima 2005). Kriging is a linear, least square estimation paradigm on which all algorithms are based. Unlike other spatial interpolation techniques, kriging provides a window for assessing the uncertainty associated with its prediction. The kriging estimator can be used to define non-stationary variable and multiple covariates (i.e. collocated kriging). Kriging technique is a generalized regression algorithm in which an unknown value is estimated from a linear combination of the primary data z 1 .
where z * 1 is the kriging surface which is a regression hyperplane fitting a ( n 1 + 1)-dimensional calibration scattergram of values z 1 (u) versus z 1 (u + h ) at distances h = u − u elsewhere, = 1, … , n 1 . In practice, this calibration is limited to two-point statistics which relates any two values z 1 (u) and z 1 (u + h) as the covariance. The weights denoted by , which determine the regression, are represented by a system of normal equations. The more robust collocated kriging (cokriging) estimator is an extension of the kriging estimator that includes other sets of data different from z 1 . In this case, if n 2 seismic data z 2 (u � ) is coupled with n 1 well data z 1 u , the cokriging estimator for the primary variable at any unknown point u can be expressed as Once again, the calibration is limited to two-point statistics which relates any two values z(u) and z(u + h) , such that z can be of type z 1 or z 2 and, n 1 and n 2 weights (1) and (2) , respectively, are given by a set of normal (cokriging) equations. The kriging estimator only differs from the cokriging estimator in a practical sense, in that instead of the one covariance function, we now have four inference and consistent modelling covariance functions (Onyekwelu 2013). The covariance functions are as follows: . It is usually difficult to handle, but with modern workstations increased dimensions of the cokriging system can be handled.
Matrix instability is one of the implementation problems associated with full collocated kriging approach when a dense secondary data is involved (Deutsch 1992). The closeness and large autocorrelation of contiguous seismic data (i.e. z 2 ), as opposed to the large separation distances and poor autocorrelation between well data (i.e. z 1 ), creates unstable cokriging matrix, i.e. close to singularity (Onyekwelu 2013). Collocated secondary data z 2 (u) with z 1 (u) value to be estimated tends to obscure the influence of farther secondary data. The solution to this problem is to retain only the collocated secondary data at each location u to be estimated, hence, making n 2 = 1 and the estimate z * 1 (u) as where m 1 = E{Z 1 (u)} and m 2 = E{Z 2 (u)} are two stationary means, and At this point, no small covariance value associated with the highly redundant secondary data exists, but it still requires the inference of C 12 (h) . The standardize form of the collocated cokriging under the Markov model (Xu et al. 1992) can be written as where the two stationary variances are 2 1 = Var{Z 1 (u)} and 2 2 = Var{Z 2 (u)}. All regression algorithms, including kriging, full cokriging, kriging with an external drift model or collocated cokriging with a Markov-type cross-covariance model, are low-pass filters which tend to yield an over-smoothed image of the actual spatial variability of the primary data z 1 . Such characteristics are suited for static volumetric calculations such as mapping of reservoir top, but may under-represent extreme values (conditional bias) when applied to dynamic modelling, e.g. spatial distribution of permeability modelling (Onyekwelu 2013). It may also remove certain patterns of spatial connectivity arising from barriers and flow paths for such primary data as permeability. To overcome these shortfalls, stochastic simulation (i.e. full-pass) algorithms that reproduce full spectrum (i.e. covariance) of spatial variability can be applied. The sequential Gaussian simulation with collocated cokriging algorithm provides a more reliable approach since it involves sequential steps in drawing alternatives, equiprobables and joint realizations of the random variable component from random function model. The various realizations represent possible images of the spatial distribution of regionalized variable over the domain. Each realization represents the properties that have been employed in developing the model. Details of this method can be found in Cekirge et al. (1981), Xu et al. (1992), Deutsch (1992).

Location and geology
The petroliferous Niger Delta Basin (NDB) lies between Longitudes 4.0° and 8.8°E of the Greenwich Meridian and Latitudes 3.0° and 6.5°N of the Equator (Fig. 1). It is regarded as one of the world's largest regressive deltas and situated along the nose end of the northeast/southwest Benue Trough-an African Cratonic quasi-linear Cretaceous sedimentary depression (Doust and Omatsola 1990;Reyment 1965;Cratchley and Jones 1965;Mascle 1976). The NDB is bounded by the Cameroon Volcanic Line to the east and the eastern-most West African transform-fault passive margin, Dahomey Basin to the west. It covers ~ 300,000 km 2 in area of land extending towards the Atlantic Ocean (Kulke 1995). Wu and Bally (2000) classified the NDB as a classical shale tectonic province due to the presence of overpressured shales and shale diapiric structures associated with the area. The overall sediment volume in the Niger Delta is ~ 500,000 km 3 (Hospers 1965) with ~ 10 km sedimentary thickness around the depocenters (Kaplan et al. 1994).
The NDB formed along a failed arm of the triple junction system (aulacogen) was initially developed in the late Jurassic following the breakup of the Gondwana into the African and South American plates (Burke 1972;Whiteman 1982). The southwestern coast of Nigeria and Cameroon which harbours two of the arms formed the West African passive continental margin, while the third failed arm developed to form the Benue Trough (Lehner and De Ruiter 1977). During Cretaceous to Tertiary time, synrift sediments were accumulated within the basin with the Albian age sediments being the oldest dated sediments. Several episodes of transgression and regression led to the deposition of marine and marginal marine sediments and carbonates (Doust and Omatsola 1990 phase. Renewed subsidence resulting from the separation of the continents paved the way for the sea to transgress into the Benue Trough. Progradation of the Niger Delta clastic wedge continued during the Mid-Cretaceous into the depocenter located on top of the deformed continental margin at the spot where the triple junction was situated. During the Late Cretaceous, sediment progradation was interrupted by episodes of transgression (Whiteman 1982).
The fault system in the Niger Delta is predominantly normal faults resulting from movements of ductile, over-pressured, deep-seated marine shale (Fig. 2). These processes over time have led to the deformation of the Niger Delta clastic wedge to a large extent (Doust and Omatsola 1990). Several of these faults which were due to slope instability along the continental margin were formed during delta progradation. These affected sediment dispersal, due to the syndepositional episodes during the basin evolution. At depths, the faults flatten onto master detachment plane near the summit of the over-pressured marine shales around the base of the Niger Delta sedimentary succession. Pockets of complex structures in isolated areas are indicators of the density and style of faulting. Flank and crestal folds (i.e. simple structures) occur along individual faults. Hanging-wall rollover anticlines were built-up, due to listric-fault geometry and differential loading of the clastic wedge above over-pressured shales. Several complex structures, truncated by group of faults with varying amounts of throw, include collapsedcrest structures with domal shape and strongly contrasting fault dips at depth. The NDB is dominated by offlap cycles of fluvio-marine sedimentary fill in a stepwise pattern. This reflects sedimentary progradation towards the south, from the edge of the continent in the direction of the oceanic basement. The Tertiary Niger Delta stratigraphic evolution and underlying Cretaceous strata are well documented in Short and Stauble (1967). The three major subsurface lithostratigraphic units of the Niger Delta are the Akata, Agbada and Benin Formations (Fig. 3). The Tertiary Niger Delta stratigraphic evolution and underlying Cretaceous strata are well documented in Short and Stauble (1967).

Materials and methods
The data set includes processed soft copies of inlines and crosslines 3D seismic sections, ten composite wireline well log data (Fig. 4) and velocity check-shot data. The entire seven suites (i.e. caliper, gamma-ray, spontaneous potential, resistivity, density, neutron and sonic logs) were used. The sonic and density logs were used to generate the synthetic seismogram which was used for well-to-seismic tie and check-shot editing. The spontaneous potential, resistivity and gamma-ray logs were used for lithologic and horizon delineation and correlation. The Schlumberger Petrel 2014.2 software and Techlog were used for data modelling. Petrophysical properties such as volume of shale, total and effective porosities, permeability and water saturation were first estimated empirically. Shale volume was estimated from the gamma-ray index by delineating the clean sand line from gamma-ray logs. For linear responses, shale volume is equal to the gamma-ray index, which gives an over estimated value of shale volume in a particular formation (Asquith and Gibson 1982). The nonlinear empirical response that provides a more reliable result (i.e. shale volume value is lower than that derived from the latter) was used. It takes into consideration the geographical area and/ or age of formation (Asquith and Krygowski 2004). Larionov (1969) correlation for Tertiary unconsolidated sands used in estimating shale volume ( V sh ) is represented as where IGR is the gamma-ray index, GR log is GR value of formation measured from log, GR min is GR minimum value of clean sand, and GR max is the GR maximum value of shale at zone of interest. Total porosity was estimated from density and sonic logs (Eq. 11).
where density of rock matrix, ma = 2.65(g∕cm 3 ) and average density of formation fluid, f = 1.0(g∕cm 3 ) . The effective porosity was estimated by introducing the volume of shale percent into Eq. 11 (Atlas 1979). Thus, The Indonesia model (Eq. 13) which can be applied nearly everywhere was used to estimate water saturation across the reservoirs (Poupon and Leveaux 1971). Fig. 3 Stratigraphic column of the Niger Delta Basin. (Redrawn from Doust and Omatsola 1990) where R w is water resistivity at formation temperature (Ω-m), R t is resistivity reading on the deep log (Ω-m), e effective porosity (fractional), m is the cementation exponent, n is the saturation exponent, R sh is the resistivity of shale (Ω-m), and V sh is the volume of shale (fractional). Permeability was estimated using Eq. 14 as where S wirr is the irreducible water saturation (Tixier 1949).
Although permeability relates to porosity, it does not totally depend on porosity, but on the volume and size of the interconnecting pores or capillaries (Nelson 1994). The measure of the potential of productive portion of the entire reservoir as percentage or decimal fraction is referred to as net-togross (NTG). The NTG ratio (which is between 0 and 1 or expressed as percentage) is the proportion of the gross rock volume that can store or produce hydrocarbon. It is usually defined by a cut-off on a permeability or porosity plot. The risk in this technique lies in the changes in recoverability during the life of the field and recovery technique employed (Worthington and Cosentino 2003). Hence, NTG can change during field life, which can be misleading (Egbele et al. 2005). NTG was estimated using Eq. 15 where the net thickness = gross thickness − shale thickness. Well correlation was performed using gamma-ray and resistivity logs, and lateral continuity of the reservoir sands was delineated (Fig. 5). Equations 9 to 15 were used to estimate petrophysical properties (Fig. 6). The seismic wavelet extraction was performed using the statistical wavelet extraction technique. Well-to-seismic tie which shows how well the synthetic trace derived from the convolution of seismic wavelet and reflectivity series from well logs correlate with seismic data was performed using synthetic seismograms of wells which had complete sonic and density logs. Figure 7 shows the modelling workflow. Quality control checks for unit consistencies and missing intervals were performed on the check-shot data, sonic and density logs prior to usage. Good ties with minimal time shift of ~ 15 ms were achieved when the synthetic seismograms were matched with the well logs apiece. Manual horizon picking which allows the interpreter to determine the directional behaviour of the tracking process was employed (Ebong et al. 2019). Both the inline and crossline of the seismic volume were mapped at intervals of 10 m, and two horizons were identified. Interpreted faults across the 3D seismic vintage were normal and listric growth faults that are characteristic of the Niger Delta Province. The convergent interpolation algorithm was utilized in the generation of time maps (Fig. 8) which were later converted to depth maps using velocity modelling processes and velocity data. Input fault sticks arising from seismic interpretation were used in building the structural models of the reservoirs. Lateral grid cell dimensions and reservoir boundaries (i.e. pillar gridding) were defined. Horizon models were built using depth surfaces with seismic and well tops as control. Each reservoir model was zoned and layered using an average thickness of 2 ft to capture the heterogeneity observed on the log curves during upscaling. The faults were modelled and pillar gridded to generate 3D grid. The modelled faults were quality checked with interpreted fault sticks to ensure consistency in geometry. A lateral grid dimension of 50 m by 50 m was used. Horizon models were built utilizing well tops and depth converted surfaces. The resultant horizons were quality checked by superimposing the contours of the input seismic interpreted surfaces of the horizon models to ensure consistency in geometry. The reservoir models were zoned using isochores generated from well tops. The zones created were further subdivided into layers. Input logs were scaled-up within the property modelling domain using the appropriate averaging algorithm. In the data analysis domain, data transformation of parameters like porosity was performed. Variogram parameters were measured in the vertical, major and minor ranges (Fig. 9). The facies logs were scaled-up into the grid cells using the "Most of" algorithm. Porosity, net-to-gross (NTG) and water saturation were upscaled using the arithmetic averaging method at all reservoir levels. Permeability data was upscaled into the 3D grid using the geometric averaging method. All the input property logs (facies, porosity, NTG, permeability and water saturation) were upscaled such that each resultant upscaled log retained the heterogeneity captured by the input logs. Facies models were built from the upscaled facies using sequential indicator simulation (SIS). Porosity was built by conditioning the facies model using sequential Gaussian simulation (SGS) algorithm. This was done because of the innate relationship between porosity and facies. Permeability was modelled and cokriged with porosity model and distributed using the SGS algorithm. By this, the permeability model will follow the porosity distribution as porosity will provide closer correlation than facie distribution (Lucia 2007). Water saturation was modelled using SGS algorithm and collocated cokriged with porosity and conditioned to facies. The hydrocarbon saturation, which is usually the difference between unity and the fraction of water saturation, was estimated. Equiprobable realizations of some petrophysical properties are shown in Figs. 10,11,12,13,14 and 15. Due to non-availability of core information, the root-mean-square error (RMSE) technique was considered as statistical validation tool to compare observed data derived from well log analysis and those predicted by the geostatistical tool (Al-Mudhafar 2017).
where X obs is observed values and X model is modelled values at time/place i.

Results and discussion
The interpretation of well logs and time surfaces revealed two reservoirs, i.e. K200 and K300 with depth to top of about 8050-8530 ft and 8190-8800 ft, respectively ( Fig. 5 and Table 1). The wells were used to delineate the various lithologic units. The lithofacies were predominantly the sand/shale sequence typical of the Niger Delta Basin. Studies show that depositional facies pattern controls the distribution of reservoir properties (i.e. lithofacies-dependent properties such as porosity and permeability), and the lithofacies defined were used to constrain the porosity distribution across the field. Although porosity varies with lithofacies, porosity statistics with lithofacies as constraint will exhibit less variability (Ma et al. 2008). Studies also show that porosity can provide more accurate correlation of permeability than the facies (Nelson 1994;Abbaszadeh et al. 2003;Lien et al. 2006;Lucia 2007). Permeability which is a measure of interconnectivity of pores can be deduced more precisely using porosity as constraint.
Hence, porosity provides more closer estimate than depositional facies based presumption, as the latter changes with time. Additionally, permeability and porosity in clastic rocks (i.e. clean sands) exhibit a direct relationship (Rahmouni et al. 2014), so anisotropy and heterogeneities in petrophysical properties are facies dependent, hence, the use of porosity as constraint in permeability estimation for both reservoirs. The variogram model parameter defined along the major, minor and vertical directions was employed in the distribution of the properties of the reservoirs based on the voxel-based method (Fig. 9). The total porosity maps generated across the two reservoirs show distinctive character based the deformational episodes and depth of burial of sediments that characterize the basin (Doust and Omatsola 1990;Owoyemi 2004). The ranges of porosity of interest that are dominant across the reservoirs lie between 0.20 and 0.35 (~ 25% on average) ( Table 2). This range of porosity with colour codes between green and yellow tends to dominate over 80% of the both reservoir units and are classified to be the high pay zones (Figs. 11 and 14). This range of porosity corroborates with reports of Okwoli et al. (2015) and Nwankwo et al. (2015) within the Niger Delta Basin. The streak of blue-pink which depicts lower porosities ≤ 0.10 shows the net pay cut-off. In clastic sediments, this range of porosity (i.e. ≤ 0.10) could be due to cementation and/ or overburden stress since porosity at depths of ~ 10,000 ft will be about this range (Larsen and Chilingarian 1982;Byrnes 1997). Within the Niger Delta Basin, it could be due to the over-pressured shales from the Akata Formation (Osinowo and Oladunjoye 2015). However, Weber and Daukoru (1975) reported evidences of overpressure at shallower depths than the Akata Formation due to rapid sediment loading of the sandy Agbada and Benin Formations on the mobile shales of Akata Formation. Conversely, some of these areas with lower porosities could be sandy-shale or silty-shale portions which may be gas prone areas. The shale volume estimated ranged between ~ 0.10 and ~ 0.25 across the field (Figs. 10 and 13). Permeability values were generally ≥ 100.0 mD (Figs. 12 and 15, Table 2). The porosity and permeability values observed were influenced by the depositional sequences and the differential responses resulting from the various diagenetic regimes within the basin. Water saturation observed across the reservoirs was ≤ 0.30 (Table 2).

Conclusion
A comprehensive petrophysical property evaluation is vital in order to optimize production and reservoir development. Geostatistical modelling affords the opportunity of evaluating petrophysical properties and extrapolating to predict those properties away from the well bore with the assistance of variogram models. The spatial and vertical variability of the petrophysical properties were estimated using seismic and well logs. The geostatistical reservoir simulation method involving the SGS algorithm was used to distribute the properties across each reservoir. The local varying mean estimator was used for porosity and other parameters except permeability which employed the collocated cokriging estimator in its distribution. These yielded high-resolution images of the petrophysical properties in the 3D grid. The validity of the estimator was tested using the RMSE techniques, and the results of the observed and estimated total porosity fitted well with minimum error of the range 0.001-0.64 (Fig. 16). Simulations were performed for each of the properties of the two reservoirs identified with three realizations each to ensure stability of the variance models. 3D maps were generated from which good and poor reservoir qualities zones were identified. K200 shows values of total porosity (0.25-0.35), permeability (> 100 mD), water