Petrophysical interpretation in shaly sand formation of a gas field in Tanzania

An onshore gas field (hereafter called the R field—real name not revealed) is in the southeast coast of Tanzania which includes a Tertiary aged shaly sand formation (sand–shale sequences). The formation was penetrated by an exploration well R–X wherein no core was acquired, and there is no layer-wise published data of the petrophysical properties of the R field in the existing literature, which are essential to reserves estimation and production forecast. In this paper, the layer-wise interpretation of petrophysical properties was undertaken by using wireline logs to obtain parameters to build a reservoir simulation model. The properties extracted include shale volume, total and effective porosities, sand fractions and sand porosity, and water saturation. Shale volume was computed using Clavier equation from gamma ray. Density method was used to calculate total and effective porosities. Thomas–Stieber method was used to determine sand porosity and sand fraction, and water saturation was computed using Poupon–Leveaux model. The statistics of the parameters extracted are presented, where shale volume obtained that varies with zones is between 6 and 54% volume fraction, with both shale laminations and dispersed shale were identified. Total porosity obtained is in a range from 12 to 22%. Sand porosity varies between 15 and 25%, and sand fraction varies between 33 and 93% height fraction. Average water saturation obtained is between 32 and 49% volume fraction.


Introduction
An onshore gas field (hereafter called the R field-real name not revealed) is in the southeast coast of Tanzania. An exploration well R-X was drilled vertically using water-based mud (WBM) to a total measured depth (TD) at 3489 m in the Cretaceous (Aptian-Albian) aged sands formation ( Fig. 1) of the R field. The commercial, sweet, dry gas was detected in the well R-X in the Miocene-Oligocene aged formations as confirmed by resistivity logs and drill stem tests. This is the target formation for petrophysical properties interpretations; however, minor oil shows were detected during drilling further down at the bottom of the well in the Eocene, Paleocene, and Cretaceous aged formations.
The target formation in the stratigraphy of the area ( Fig. 1) can be categorized as shaly sand (silts and clays with thin turbidite sands) of deltaic marine deposits inside the canyon slope settings (TPDC 2011; RPS Energy 2018) as shown in Fig. 1b. In the regional stratigraphy (Fig. 1a), the interval additionally includes reefal and chalky limestones. No core was acquired in the target depth interval. The section below the target interval (Eocene, Paleocene, and Cretaceous formations) consists of thin sands interbedded with shales and extrusive igneous rocks.
Petrophysical properties (e.g., porosity, permeability, water saturation, and net reservoir thickness) are fundamental to reserves estimation (or static property modeling) and production forecast which is equally important for R field. These properties can be derived from well logs, formation flow (pressure) tests, or cores through formation evaluation, or combined cores with wireline logs and/or formation pressure test(s) analysis, which is the general practice (Hill 2017;Ellis and Singer 2007;Tiab and Donaldson 2016). Formation evaluation using wireline logs is an appropriate means to obtain petrophysical properties for the R field with well R-X in the target interval where no core was acquired.
The primary well-logs' measurements (properties) used for porosity determination are formation density and nuclear magnetic resonance (NMR). Other commonly used measurements are resistivity, neutron, sonic (acoustic), and dual combinations (cross-plots) of density, neutron, and sonic logs (Hill 2017;Tiab and Donaldson 2016;Ellis and Singer 2007). Properties accounting for and included in models to compute water saturation are electrical conductivity (inverse of resistivity), porosity, formation water salinity, lithology-dependant fitting parameters (cementation or porosity exponent m and saturation exponent n), and temperature (Thomas 2018;Santamarina et al. 2019;Hill 2017;Archie 1942). These properties can vary vertically (with depth) and laterally depending on the lithology-type and texture with the degree of sorting, compaction, cementation (e.g., type, quantity, and distribution of shale or clay minerals), and dissolution (Thomas 2018;Lin et al. 2017;Gong et al. 2016;Zheng et al. 2015;Magara 1980;Smith 1971;Athy 1930). Porosity (or compaction) of sedimentary rock decreases with increasing depth (overburden pressure) whereby for shale it decreases rapidly at shallow depth and slowly at deep depth. The decrease in porosity for the sandstones at shallow and deep depth is fairly constant, wherein cementation can significantly affect porosity change (Hakiki and Shidqi 2018;Magara 1980;Smith 1971).
Permeability is derived mainly using NMR, formation testers, and cores; however, reasonable estimates can be obtained using porosity-permeability relationships. The net reservoir thickness (or net to gross ratio) is essentially derived from well logs with cutoff values of shale volume, porosity (and/or permeability), and water saturation.
In the existing literature, there is no layer-wise published data of the petrophysical properties in the R field. Therefore, the layer-wise interpretation of petrophysical properties was undertaken to obtain inputs to a rock property model to build a reservoir simulation model. The properties can also be used to estimate saturation-height relationship(s) (capillary pressure) and relative permeability. Measurements (variables) available for the study are gamma ray, bulk density, compensated neutron, and resistivity (induction, laterolog, and microresistivity). These data were environmental corrected. The properties estimated include shale volume, total and effective porosities, sand fractions and porosity of sand, and water saturation. The shale volume with its distribution and influence on the porosity change was elucidated.
Established methods were used in the evaluation of the properties considered. The equation of Clavier et al. (1971) was used to compute shale volume from gamma ray. Density and neutron-density methods were used to calculate total and effective porosities. Thomas and Stieber (1975) method was used to determine shale distributions, sand fractions, and sand porosity. This method was selected to account for the thin-bed effects on the log measurements in the thin turbidite sands. Water saturation in the sand was computed using the model of Poupon and Leveaux (1971). The model includes shale volume and shale resistivity to account for shale effect to the formation conductivity which is relevant for this interpretation. The details of established methods are presented in "Appendix A".

Methodology and workflow
Methodologies used were implemented in Techlog software, and a summary of a workflow is represented in Fig. 2. Details of the workflow are as follows:

Zonation
Eleven sand (shaly sand) zones were identified between 1493.22 and 2150.06 m depth from the rotary Kelly  bushing (RKB) which include the reservoir interval. The lithological indicators used to identify these zones were gamma ray responses and separation of density and neutron log curve. Relative low gamma ray responses and negative separation of a plot of density and neutron curves in a compatible scale can be an indicator to the presence of sandstone or shaly sand formations (Ellis and Singer 2007). The thickness of the zones ranges from 1.8 to 16.9 m. The depth interval and thickness of each zone are presented in Table 1, and sample zones are shown in Fig. 3. There is a thick shale layer (19-147 m) between zones as

Formation temperature
Formation temperature was determined using the Horner plot ( Fig. 4) with maximum well recorded temperature during logs run. Real formation temperatures from the Horner plot were used to calculate temperature gradient (0.0305 °C/m) presented in Fig. 5. This gradient was used to calculate a continuous profile of bottom-hole temperature (BHT) for calculation of temperature-dependent properties.

Shale volume computation
Shale volume was computed using Clavier equation (Eq. 1) and neutron-density (Eq. 2) methods. Properties of 100% shale were taken at a depth interval between 1413 and 1485 m RKB. This is an interval just above the top zone selected and has an even shaped plot of gamma ray log. Properties of 100% sandstone matrix were taken at 1994.24 m RKB depth. These properties of 100% sandstone, shale, and water are represented in Table 2.

Fig. 4
Horner plot of bottom-hole temperature; where markers (black circles, red triangles, blue diamonds, and magenta crosses) are maximum well recorded temperature at depths shown, and solid lines are logarithmic least-squares fit to the data. At least three maximum thermometer records of the logging tool are preferable to extrapolate formation temperature using the Horner plot

Total and effective porosity
Fluid densities for porosity calculation were determined using pressure gradients (∇p's) of 0.0991 bar/m (0.438 psi/ft) and 0.0131 bar/m (0.058 psi/ft) for water and gas, respectively, obtained from RPS Energy (2018) report. The in situ ∇p's were determined by using pressures which were recorded during repeat formation tests (RFTs). The ∇p's were converted to densities using the relationship ∇p = f g , where ρ f is the fluid (water or gas) density and g is the standard gravity (9.80665 m/ s 2 ). The corresponding values of densities obtained are represented in Table 2. Effective fluid density in flushed (invaded) gas zones for porosity calculation was estimated using invasion correction Eq. 3 with assumed flushed zone water saturation (S xo ) of 0.5-0.8. The invasion was considered to be proportional to the mudcake thickness (an analogy of mud filtrate) from the caliper log. Bulk total porosity and effective porosity of each zone were computed using density (Eqs. 4a-4c) and neutron-density porosity equations (Eqs. 5a-5d). Density total and effective porosity were used in the next calculations.

Evaluation of sand porosity and sand fraction
Total porosity and shale volume were used as inputs in Thomas and Stieber (1975) method to compute sand porosity and sand fraction (VSS) as shown in Fig. 6. The method assumes that within the interval investigated there are only two rock types, a high porosity clean sand and a low porosity pure shale, and in situ porosity is a result of mixing the two rocks. Also, pure shale sections above and below the sand and shale mixed in the sand are mineralogically the same. These assumptions were adopted in this interpretation. Shale volume from Clavier equation was used in this evaluation. Clean sand porosity of 26% and pure shale porosity of 16.4% were taken from depth intervals mentioned in shale volume computation Section to develop endpoints of Thomas-Stieber plot. Other parameters (clean sand and shale gamma ray, and shale resistivity) which are required to develop Sensitivity analysis on porosity of sand and sand fraction was performed by changing shale endpoint in two ways. First scenario involved increasing porosity of shale in the base case by 10% (Fig. 7a). Second scenario, the gamma ray in shale was reduced by 4% (Fig. 7b) from the values used in the base case (Fig. 6). The results are represented in Fig. 8a and b for the sensitivity on sand porosity and sand fraction, respectively.

Computation of water saturation
Water saturation was computed using Archie (1942) equation (Eq. 6) and Poupon and Leveaux (1971) equation (Eq. 7). Density total porosity was used in the Archie (1942) equation and density effective porosity in the Poupon and Leveaux (1971) equation. The formation water resistivity for water saturation calculation was determined by solving Archie (1942) equation in water zones (S w = 1) with assumed tortuosity factor (a) equal 1.0, saturation exponent (n), and cementation exponent (m) equal 2. Note that, Picket plot or ratio method for R w calculation (Eq. 8) is not applicable in the depth interval interpreted because microresistivity measurement is missing. The R w value of 0.1778 Ωm was obtained at 93.06 °C. Salinity of formation water, 12,612 ppm, was estimated using the Schlumberger (2013) chart Gen-6, and then, a profile of R w was generated by using the Arps (1953) equation given on the chart. The results of shale volume, sand fraction, porosities, and water saturation computed are displayed in Fig. 9 for sample zones.

Results and discussion
Bulk shale volumes computed are represented in Table 3 in terms of P90, P50, P10, and means. The probability values reported in this interpretation refer that a P90 is a property estimate at the low end of the range (descending cumulative density function). The higher mean value of the shale content of approximate 0.5 m 3 /m 3 was obtained in zones H, J, and K, and lower mean value of approximate 0.1 m 3 / m 3 was obtained in zones A and B. The deviation of probabilistic values of shale content obtained was relatively high (0.46, 0.41, 0.39, and 0.42 m 3 /m 3 ) in zones C, F, G, and K, and relatively low (approximate 0.1) in zones A and I.
The mean value of total porosity was 21.8% in zone E which is relatively higher compared to a mean total porosity of 12.4% in zone D (Table 3). The difference between P10 and P90 of porosities is about 4% in zone I which is relatively low compared to 21% in zone A. The high difference between P10 and P90 in zone A is because of low porosities (about 4-10%) at the bottom of this zone. These low porosities were expected because of relatively higher log readings of bulk density (about 2.4-2.6 g/cm 3 ). Both neutron (less than 10%) and gamma ray (less than 56 GAPI) log readings at the lower part of zone A are relatively low which suggest that the formation is probably cemented, or the shale distributed in dispersed form. The same reason and possibility are provided for the cause of low mean total porosity in zone D.
Sand porosities and sand fractions of the Thomas-Stieber analysis are represented in Table 4. High mean sand porosity of 24.6% was obtained in zone H of which its average sand fraction is approximate equal 0.4 m/m (Table 4). Zone A has high sand fraction (mean value of 0.93 m/m), but sand porosity has a mean value of 22%. A plot of zones mean values of total porosity and gamma ray in the Thomas-Stieber diagram is presented in Fig. 6, where the locus of zones A Essentially, the properties of the sand and shale endpoints are quantified though detailed core analyses in practice. Without detailed core analyses, endpoints can be determined as described in the methodology which are subjected to uncertainties that affect the evaluation of sand fraction and sand porosity if the assumptions of the model are correct. The shale endpoint is generally more uncertain than sand endpoint because of variation of shale matrix density (2.12-3.10 g/cm 3 ) with the clay mineral groups, while the sand (quartz) grain density is typically fixed (2.65 g/cm 3 ) (Craddock et al. 2018;Serra 1984). Therefore, sensitivity analysis on sand fraction and its porosity was performed by changing shale endpoint.
Sensitivities of sand porosity and the sand fraction associated with changes in shale point ( Fig. 7a and b) are represented in Fig. 8a and b, respectively. The 10% increase in shale porosity can decrease the sand porosity of zone C by 0.042 m 3 /m 3 (Fig. 8a), but in other zones, the decrease in sand porosity is less than 0.042 m 3 /m 3 . The change in sand porosity because of the 4% decrease in gamma ray of shale is less than ± 0.033 m 3 /m 3 . Therefore, sand porosity is relatively more sensitive to 10% increase to shale porosity than 4% reduction in the shale gamma ray. Sensitivities of sand fraction to changes in shale point are represented in Fig. 8b. The change in sand fraction because of 10% increase in shale porosity within ± 0.042 for all zones. The change in the sand fraction because of the 4% decrease in gamma ray of shale is between − 0.067 and − 0.007.
The porosity-depth relationship is nonuniform as presented in Fig. 10. For the top three zones (less than 1700 m depth), it is obvious that porosity is decaying with increasing depth and compaction with cementation may have caused porosity loss. At the depth greater than 1700 m, the top trend is interrupted possibly because of the additional diagenetic effect (e.g., dissolution); however, afterward porosities appear to decay. The relatively higher porosity loss in zone D is because of the high degree of cementation. Hitherto, the sand porosity has been evaluated by correcting the shale laminae effect, and the porosity-depth relationship is elaborated in this paper. Further studies are required to investigate the porosity history and factors involved.
The results of water saturation are represented in Table 5 where cutoffs of V shale < 0.50, φ e > 0.08, and S w < 0.60 were used to derive the distribution in gas zones. The cutoffs were adopted from RPS Energy (2018) report, which have been applied in evaluating petrophysical

Conclusions and future perspectives
The petrophysical properties (shale volume, total and effective porosities, sand fractions and sand porosity, and water saturation) were extracted in the layer-wise of R field using wireline logs which were not available at these details in the published literature. These properties can be used in rock property modeling to build a reservoir simulation model, to estimate saturation-height relationship(s) (or capillary pressure) and relative permeability. Hitherto, the layer-wise statistics of the properties can be presented.
The mean value of bulk shale volume varies between 6 and 54% volume fraction across zones. The P90 of shale volume varies between 2 and 40%, and P10 varies between 13 and 74%. Shale laminations and dispersed shale were recognized in turbidite sand layers. The mean of total Fig. 9 Display of results in sample zones. The parameters contained in the tracks from the left of the plot are measured depth (MD) in track 1; sand fraction from Thomas-Stieber method (VSS_TS) and shale volume (VSH_FINAL) in track 3. Sand porosity from Thomas-Stieber method (POR_SS_TS) and total porosity from neutron-density and density methods (PHIT_ ND, PHIT_D), respectively, are presented in track 4. Effective porosity from neutron-density and density methods (PHIE_ ND, PHIE_D), respectively, is presented in track 5; and water saturation from Archie and Indonesian equations (SW_AR, SWE_INDO), respectively, is presented in track 6 porosity varies between 10 and 21%. The P90 of total porosity is between 6 and 18%, while P10 is between 16 and 28%. The mean of effective porosity is between 11 and 21%; P90 is between 9 and 14%; and P10 is between 15 and 27%. The mean of sand porosity is between 13 and 22%; P90 is between 11 and 18%; and P10 is between 17 and 27%. The mean of sand fraction is between 0.35 and 0.94; P90 is between 0.21 and 0.88; and P10 is between 0.50 and 1.0. Average water saturation in gas zones is between 0.32 and 0.49; P90 is between 0.19 and 0.45; and P10 is between 0.41 and 0.60.
A verification is made in the promising use of the Thomas-Stieber method to avoid bypassing thin (less than 3 m) sand layers in the R field and identifying shale distributions.
The sand porosity was determined by using the Thomas-Stieber method with the assumption that shale inclusion is the main porosity destroying factor in Tertiary aged sand-shale sequences. However, other factors (e.g., mineralogical-type, texture, sorting, and compaction) can influence porosity change. Further studies can be undertaken to investigate the porosity history and factors involved.  2013-2019. Financial support of the EnPe program is gratefully acknowledged. The authors would also like to acknowledge the Tanzania Petroleum Development Corporation (TPDC) for the permission to use and publish data of the R field in this project.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Clavier equation for shale volume
Where , V shale is the volume fraction of shale in a zone of interest. GR ma is a gamma ray log reading in 100% matrix rock, GR shale is a gamma ray log reading in 100% shale, and GR log is a gamma ray log reading in a zone of interest (Bassiouni 1994;Clavier et al. 1971).

Neutron-density equation for shale volume
which is a slope of a line joining 100% water point and 100% matrix point. NPHI log and ρ log are neutron porosity and bulk density, respectively, logs reading in a zone of interest. NPHI ma and ρ ma are neutron porosity and bulk density, respectively, logs reading in 100% matrix rock. NPHI f and ρ f are neutron porosity and bulk density, respectively, logs reading in 100% water. NPHI shale and ρ shale are neutron porosity and bulk density, respectively, logs reading in 100% shale (Schlumberger 1991).
(3) f = mf S xo + g 1 − S xo , Fig. 10 Porosity against depth plot, where total porosities (black-circle markers) were derived from density log and orange-triangle markers are the shale laminae-free sand (or shaly sand) porosities derived from the Thomas-Stieber method where S xo is a water saturation in flushed zone, ρ g is a density of gas, and ρ mf is a density of mud filtrate.

Density porosity equation
Algorithm to calculate porosity from neutron-density cross-plot (Schlumberger 1991) If mf < sand , mf ⟨ b and b ⟩0 Where: ρ b is bulk density, ρ lim limestone grain density (typical value is 2.71 g/cm 3 ), ρ mf is mud filtrate or bulk fluid density (typical value is 1 g/cm 3 ), ρ sand is sandstone grain density (typical value is 2.65 g/cm 3 ), ρ dol is dolomite grain density (typical value is 2.9 g/cm 3 ), φ n is neutron porosity log reading, φ d is density porosity.
The condition to decide the lithology line is such that if n ≤ d , limestone-sandstone combination is selected, and if n > d , limestone-dolomite combination is selected. A test is performed for anhydrite such that if n > d , 2.91 ≤ b ≤ 3.5and e ≤ 0.04 . The calculated porosity should be equal to zero (0). Otherwise, the calculation proceeds as follows: firstly, the matrix is converted from limestone to sandstone or dolomite, depending on the lithology line combination decided above. Secondly, apparent hydrogen index (HA) is computed using neutron porosity with porosity of sandstone or dolomite matrix ( VA ): The cross-plot porosity is calculated iteratively using the corresponding line combination equations. The iteration converges when the difference of the latest calculated porosity from previous porosity is less or equal to 0.01.

Archie (1942) equation
where a is a tortuosity factor, m is a cementation exponent, n is a saturation exponent, R w is formation water resistivity, R t is formation resistivity, respectively.

Formation water resistivity (ratio method)
where R mf is a resistivity of mud filtrate and R xo is a resistivity of the invaded zone.