1. Introduction

We thank the Department of Geosciences at the Universidad de los Andes, Colombia, and the operating company for authorization to use the data. We also thank the editorial group of the Boletín Geológico, and the anonymous reviewers who helped improve the manuscript with their valuable comments.

The Llanos Orientales is one of the most important oil basins in Colombia. The production of hydrocarbons and economic profitability of the development of fields in this basin have contributed to the country’s energy supply. For this reason, it is essential for the oil industry to determine the volume of hydrocarbons in a reservoir. This is done through geologic modeling of the reservoir, with which geological and petrophysical aspects are determined and the volume of hydrocarbons in place is calculated.

During the drilling of a well, information is acquired to determine the geological and petrophysical properties of the subsoil by direct and indirect methods. Direct methods consist of core sampling and drill cuttings, corresponding to the extracted rock. These samples are studied to determine the lithology of the drilled formation and the cores are subjected to laboratory tests for direct measurement of the petrophysical properties of the formation, such as porosity and permeability.

Indirect methods consist of quantitative measurements of formation properties. These include electric well logs of gamma rays, density, and resistivity, using a sensor and receiver. These records provide indirect values of fluid saturation in the rock and serve as indicators of the lithology.

The main objective of the present work was to calculate the original oil in place (OOIP) from the geological modeling of a reservoir in two fields, with interest in unit C7 of the Carbonera Formation in the Llanos Orientales basin. The specific objectives were as follows: to conduct stratigraphic correlation of the major units of interest and a structural framework from the regional datum; to study the electric well logs for petrophysical evaluation; to analyze the results of that evaluation to characterize the operational units of the reservoir used in the modeling and perform the volumetric calculation of the OOIP in the area of interest.

To do the above, electric well logs and structural surfaces were used to generate a 3D model of the reservoir. This permitted the calculation of the volume of hydrocarbons via the Petrel software. For the work, data from seven wells were used, which were provided by the operating company. Five wells correspond to Field 1 and two wells to Field 2.

2. Geological setting

The study area is in the northeastern Llanos Orientales basin of eastern Colombia, in the department of Arauca, near the municipality of Arauquita (Figure 1). Llanos Orientales is classified as an asymmetric foreland basin, defined as an elongated region formed between the orogenic belt of the Eastern Cordillera and the Guiana shield (Duarte et al., 2017).

The source rocks, deposited during the Cretaceous, are of marine origin and were buried during multiple orogeny events, which favored the formation of structural and stratigraphic traps (Campos and Mann, 2015). Source rocks of the Lower Cretaceous produced hydrocarbons over time, supporting substantial accumulations in Llanos Orientales basin (Mora, 2015).

The structural geology of the study fields corresponds to an anticline elongated in the northeast direction with closure in three directions, limited to the east by the La Yuca regional fault, which presents a strike-slip. The two fields of study are limited by a normal and local fault that is perpendicular to the main La Yuca fault.

Carbonera Formation. The Carbonera Formation dates from the Miocene at the top and early Eocene at the bottom. In the study area, at the northeast of the basin, this formation has an average thickness of 7053 feet, and is composed of sandstone and shale intercalations that permit division of the formation into eight units. Among these, the even numbers correspond to clay lithologies and the odd ones mainly to sands (Table 1). The formation was deposited in continental-type environments to the east and southeast of the Llanos Orientales basin, and in transitional environments to its west and northwest (Domínguez, 2014).

Figure 1.
Location of study area

a) Red polygon demarcates Llanos Orientales basin and green dot the location of the study area; b) Location of Llanos Orientales basin within Colombia.

Table 1.
Study area column

Lower Carbonera. During the early Oligocene, regional sandstone deposition was caused by an increase in the supply of sediments, which facilitated the deposit of Unit C7 of the Carbonera Formation (De la Parra et al., 2015). This was integrated in only one unit called the lower Carbonera, composed of fine-grained facies, caused by a prograding wedge of braided fluvial deposits (Reyes-Harker et al., 2015) and facies of medium grain.

In the study area, the C7 unit has a thickness of 100–300 feet and the companies divide it into the operational units M1, M2 and M3, which consist of massive sandstone bodies separated by a regional shale called guafita, with intercalations of sandstones-siltstones. The sandstone deposits correspond to a fluvial-deltaic environment that accumulates hydrocarbons in the operational production unit M1.

The marker of the area of interest is the guafita shale, which corresponds to the maximum flooding surface that is homogeneously distributed in the northeastern basin. This marker was used as a regional datum because it was identified in all study wells.

3. Method

We used information from seven wells, five in Field 1 and two in Field 2, that have directional files, electric well logs of gamma ray (GR), resistivity (RT), and density (RHOB). Likewise, we took as reference the regional datum and its structural surface, the faults in the two fields of study, the permeability equation as a function of effective porosity found from data of tests on rock core plugs, and the J functions determined from reservoir capillary pressure curves. The cores were not available for facies analysis. In the following, we describe and explain the steps in the study.

3.1. Stratigraphic model

We conducted a stratigraphic correlation based on the regional datum of the study area. The datum was located at each well and the well tops of the operational units that belonged to the C7 unit of the Carbonera Formation were defined as: M1, M2_1, M2_2, and M3. These were identified based on the response of the GR log.

In the analysis of electrofacies, we identified the facies of the depositional environment of the studied intervals according to patterns of electric well logs established and accepted in the geological literature.

To determine the facies, we analyzed forms of the GR log, which indicated a relative variation of grain size related to the depositional environment. The large values of the log indicate the presence of clays and small values the presence of sands. The main forms of the log responses are divided into cylinder, funnel, and bell. Each form can be associated with more than one depositional environment.

3.2. Structural model

The structural surfaces of operational units M1, M2_1, M2_2, and M3, were generated from the structural surface to the top of the guafita regional datum and isopachous maps of each unit that were obtained from the correlation between wells. We created fitted surfaces to the well tops of the units by adding or subtracting the average thickness of each surface. This operation implies an addition of thicknesses if the unit is above the guafita surface, or a subtraction if it is below.

3.3. Fault model

We used fault segments supplied by the operator company as a result of seismic interpretation, which correspond to deep fault lines of the reservoir. By interpolating these segments, we obtained the plane of the fault in depth. The segments become fault pillars, a set of lines with upper, middle, and lower nodes that define the shape and location of the fault surfaces in the model.

3.4. Petrophysical model

The petrophysical properties used for modeling the reservoir were effective porosity, permeability, and water saturation. For the porosity model of units M1, M2_1, M2_2 and M3, we used effective porosity (ϕe) furnished by the operating company, which is calculated from the total porosity (ϕT) based on the density log and corrected by the clay volume (Vsh). This model is related to the reservoir facies. The permeability model (k) is based on an equation that relates porosity:

Here, m is the slope and c a constant that indicates the cutoff point with the permeability, when porosity has a value of zero (Satter et al., 2008). This equation is established from core samples in the laboratory. Then, based on the permeability equation, the modified Lorenz method was applied (Craig, 1972) to determine the various rock types in the reservoir.

Rock typing consists of grouping rocks that have similar petrophysical properties of the different flow units comprising the reservoir (Guo et al., 2005). For this, we used the modified Lorenz method, which normalizes and graphs values of the effective porosity and permeability. Each slope on the graph indicates a type of rock.

Thereafter, we grouped the rock types of the formation using the flow zone indicator (FZI) method, which uses the fraction of effective porosity ϕz and the rock quality index RQI, based on the relationship between permeability and porosity (Lis-Śledziona, 2019).

The values from Equation (4) are grouped to classify and graph the effective porosity and permeability ranges for each rock type. The trend equation is a mathematical expression of permeability as a function of porosity. A permeability function was established for each rock type via a permeability vs. porosity graph, and using a software calculator, the functions were applied to generate the permeability model.

The water saturation model starts from the rock types and water-oil contact, related to the capillary pressure curves and functions J(Sw) provided by the company. These functions were determined from the capillary pressure curves of the cores, which were averaged to obtain a representative curve for the area of interest, called the J function. This is a dimensionless function that correlates the capillary pressure to the reservoir pore structure in terms of porosity and permeability (Donnez, 2007). The functions are specific to each rock type and are implemented in the software to generate the water saturation model of each field.

Here, J(Sw) is the J function, Pc the capillary pressure, k the permeability, ϕe the effective porosity, ϑ the contact angle, and σ the interfacial tension.

3.5. Construction of 3D geologic model

Using the previous models, we constructed the geologic model with standardized and normalized data. This model includes the construction of a cell grid with dimensions X, Y and Z, limited by the structural surfaces at the top of each unit (M1, M2_1, M2_2, and M3) and by the fault surfaces.

We constructed the model grid using the fault surfaces as the framework. This grid was divided into three skeletons, upper, middle and lower, related to the nodes of the fault columns. The grid is composed of cells with a rectangular prism shape that the software defined in the I, J and K directions, corresponding to the directions X and Y (I and J) with depth K, equivalent to Z, with sizes defined on each axis. The grid has 337 cells in direction I, 334 in direction J, and 178 layers in direction Z. This results in a total of 20 035 324 cells, corresponding to a geocellular model of high resolution.

Horizons and layers were created by incorporating structural surfaces to the top of each unit (M1, M2_1, M2_2 and M3) in the 3D grid, and vertically subdividing each unit limited by surfaces by a determined number of layers. The number of vertical layers of each unit was determined from the minimum thickness to preserve geological bodies of minimal thickness, which corresponds to the vertical resolution of the model.

After determining the facies, we did the upscaling, which consists of taking information from the scale of well logs to that of the geologic model. These scaled cells were the main input to the 3D modeling. Each cell traversed by a well corresponds to only one value of the scaled property of the wells, either discrete or continuous. For facies, these are discrete quantities because there can only be an integer value associated with a facies. Petrophysical properties are continuous variables that take on decimal values within a specific range. The scaled properties in this case are the facies and effective porosity. For the permeability and water saturation, we used equations implemented in the software calculator, which are a function of the previously established properties.

After scaling the facies in each well, we used a geostatistical analysis module, which was applied to the discrete values of the facies in order to spatially populate the cells. This module is based on the deterministic or stochastic modeling algorithms – sequential indicator simulator (SIS) and object modeling. The SIS algorithm generates a pixel model of the facies distribution, following the probability tendencies of these and extrapolation using variograms obtained for each facies of each unit. Object modeling represents the geometry of geological bodies such as channels, bars and fans, in order to reliably represent the distribution of those bodies. This information is acquired via an analysis of seismic interpretation (seismic attributes), data of outcrops and/or known analogue fields, which facilitate acquisition of the directions and dimensions of modeled bodies.

After modeling the facies, we modeled the petrophysical properties. For the porosity model, we used the scaled log and conducted a geostatistical analysis, generating variograms for the porosity of each facies of each unit. The permeability was modeled using an equation related to the reservoir porosity, from which was determined the rock types. Each type had its range of permeability and porosity. Then we found a function for each rock type and extrapolated to the model by applying this equation in the software calculator. The J functions of unit M1 were used to obtain the water saturation model, which is based on the identified rock types and permeability model.

3.6. Contacts and volumetric computations

For each field, we determined the water-oil contact through the resistivity log and tests of well production. We used the volume calculation module for the OOIP and selected the oil-water contact, water saturation model, porosity model, facies model, and oil volumetric factor. From the above models, the software calculated the OOIP of each field.

4. Results and discussion

4.1. Facies model

The study interval corresponds to the operational units M1, M2 and M3, among which M1 is the unit of interest. A stratigraphic correlation was carried out at the level of these units, taking as a starting point the regional datum of the study area. Figures 2 and 4 show this correlation, wherein the tops M1, M2_1, M2_2 and M3 are identified based on the response of the GR and density logs. Unit M2 corresponds to thin intercalations of sands and shales. Unit M2_1 consists of clean sands, according to smaller GR values. Unit M2_2 had larger values in this log, owing to the clay content in the sands. Units M1 and M3 are relatively homogeneous and correspond to packets of sandstone with intercalations of shales. We are only interested in the top of unit M3, so it was assigned a lower limit corresponding to the base of the model.

Based on regional studies and those of the cores of neighboring fields of the operating company and of similar fields, it was established that the studied reservoir was deposited in

Figure 2.

Figure 2.
Stratigraphic correlation of the seven wells

Field 1: P1, P2, P3, P4, P5; Field 2: P6, P7

a fluvial-deltaic environment. Four main facies were determined, channel, deltaic front, overflow fan, and floodplain. Figure 3 shows the facies model.

The channel facies correspond to massive sand bodies of units M1 and M3, and smaller sand packets for the channels of M2_1. The signature of clean sand channels in the GR log were identified with cylinder shape, indicating small GR values and uniform lithology. These facies are the most important of the study because they have greater flow and storage capacity, owing to their petrophysical properties of porosity and permeability.

Floodplains are flat areas adjacent to the round of a river, or long-range areas composed of fine sediment lithologies. The floodplain facies are identified by large values in the GR log, which correspond to areas of fine-grained lithology such as shales and are local seals of the reservoir. These facies are present in operational units M1, M2_1 and M2_2, and to a lesser extent in M3.

An overflow fan is generated when a river breaks its barrier and pulses of sand are released, creating fan-shaped sedimentary bodies. These bodies are recognized in the logs by a sawn-off funnel shape with a coarsening-upward sequence. This is because they have fine-grained deposits at the base and, as the river

Figure 3.
Plan view of layer 1 facies modeling

Figure 4.

Figure 4.
Stratigraphic correlation with the relationship of the facies defined from seven wells, based on the response of well logs

Field 1: P1, P2, P3, P4, P5; Field 2: P6, P7

releases sands, they increase toward the top. These facies are present in the upper part of unit M1 and are in a fluvial environment.

The deltaic front facies mark a transition zone between a continental and marine environment. The deltaic front corresponds to the shallow part of a delta receiving sediment from river channels and, for this reason, in the log analysis there was a coarsening-upward response in the shape of a non-sawed-off funnel. It was also identified that they are heterogeneous bodies that are interspersed with channels due to the channel-front interaction of waves and sediment redeposition. These intercalations are evident in unit M2_1.

4.2. Structural model

The structural model of the fields corresponds to an anticline elongated in the northeast direction with closure in three directions. It was obtained from the structural surface in depth at the top of the guafita marker and the faults affecting it. The top surfaces of each unit were created using isopachous thicknesses. Figure 5 shows a 3D view of the structural surface at the top of the producer unit M1.

4.3. Petrophysical model

The modeled porosity is discriminated by the facies of each unit. Figure 6 shows a plan view of this variable, where it is seen that the model adequately represented the porosities of the facies. The channel facies has high porosity values of 22%–27%, the delta front facies 18%–22%, and the facies floodplain < 5%.

To find the model permeability (K), we used the data on reservoir permeability furnished by the operating company.

We determined the rock types of the reservoir with the modified Lorenz method. Figure 7 shows three different slopes, corresponding to rock types categorized in ascending order of porosity: Rock 0, related to floodplain facies; Rock 1, related to channel facies; Rock 2, related to deltaic front facies.

We obtained an FZI range for the Rock 2 type at 0–1.65, for Rock 1, 1.66–11, and for Rock 0, 12–8312. For the permeability model, equations were established for each rock type using a graph of permeability vs. effective porosity.

Figure 5
Structural surface to the top of producer unit M1.

Figure 6.
Plan view of porosity modeling, layer 1

Figure 7.
Modified Lorenz graph

A permeability model discriminated by rock type and effective porosity was constructed. Figure 8 shows a plan view of the modeled permeability, where it is seen that the facies model and rock types are coherent. The high permeabilities (> 100 mD) correspond to rock types 1 and 2, which are related to the channel and deltaic front facies.

The following table shows average porosity and permeability values for each facies present in each unit.

The water saturation model has as input the permeability model by rock type, executing the J functions for operational unit M1. Based on information of pressure, logs and production tests, we determined the water-oil contact at Field 1 at a true vertical depth subsea depth of −7848 feet. At Field 2, this was −7777 feet.

Table 2.
Average values of porosity and permeability models for each facies in each unit

J functions were run for each field in the software calculator. Figures 9and 10 depict the water saturation model for unit M1 of the two fields studied. The area of interest corresponds to saturations from 20% to 30% in the structurally higher part of the reservoir, corresponding to channel facies.

4.4. Calculating OOIP volume

The software included the previous models to do the volume calculation of the OOIP of unit M1 of each field. An oil volume factor of 1.06 was used for both calculations. For Field 1, we obtained a value of six million barrels of OOIP, and 9 million barrels for Field 2. Figures 11and 12 show OOIP maps, indicating the area of oil accumulation.

The values obtained were used to calculate the current recovery factor, based on August 2019 production data. For this we divided the accumulated production of each field by the calculated OOIP. Field 1 has a current recovery factor of 1.5% and Field 2 has 22%.

5. Conclusions

From the stratigraphic correlation of the wells, we identified four units of the lower Carbonera Formation, M1, M2_1, M2_1, and M3. The unit of interest was M1, because in it lie clean sands that store hydrocarbons.

The deposit environment of the reservoir is fluvial-deltaic and the most important facies is that of the channel in operational unit M1. There, clean sands predominate, which are identified as the main pool of the reservoir. The floodplain facies are hallmarks of the reservoir that are present mainly in units M2_1 and M2_2. Facies of the deltaic front reservoir and overflow fan are lower quality reservoirs in comparison with the channel facies.

Figure 8.
Plan view of layer 1 permeability modeling

Figure 9.
Layer 1 of water saturation modeling for unit M1 of Field 1

Figure 10.
Layer 1 of water saturation modeling for unit M1 of Field 2

Figure 11.
Field 1 OOIP Map

Figure 12.
Field 2 OOIP Map

Petrophysical properties such as effective porosity indicate that unit M1 is the most prone to accumulate hydrocarbons, as it has the greatest porosities in the channel facies (22%–27%). We integrated the stratigraphic, structural, petrophysical models through 3D geological modeling, from which the OOIP was calculated in each field.

Field 1 has 1.5% recovery factor. This indicates that the sands of unit M1 can produce oil with recovery factors up to that of Field 2, which is 22%.

Acknowledgments

We thank the Department of Geosciences at the Universidad de los Andes, Colombia, and the operating company for authorization to use the data. We also thank the editorial group of the Boletín Geológico, and the anonymous reviewers who helped improve the manuscript with their valuable comments.

References

Campos, H., & Mann, P. (2015). Tectonostratigraphic Evolution of the Northern Llanos Foreland Basin of Colombia and Implications for Its Hydrocarbon Potential. In Petroleum geology and potential of the Colombian Caribbean margin (pp. 517-546). American Association of Petroleum Geologists. https://doi.org/10.1306/16531948M1083651

Craig, F. F. (1972). The Reservoir Engineering Aspects of Waterflooding. SPE Monograph 3. Society of Petroleum Engineers.

De la Parra, F., Mora, A., Rueda, M., & Quintero, I. (2015). Temporal and spatial distribution of tectonic events as deduced from reworked palynomorphs in the eastern Northern Andes. AAPG Bulletin, 99(8), 1455-1472. https://doi.org/10.1306/02241511153

Domínguez, A. (2014). Análisis Sismoestratigráfico de la Formación Carbonera Miembro C7 (Municipios de Orocué y San Luis de Palenque, Cuenca Llanos Orientales), Colombia [Undergraduate thesis]. Universidad Nacional de Colombia, Bogotá, Colombia.

Donnez, P. (2007). Reservoir rock properties. In Essentials of reservoir engineering (pp. 19-20). Editions Technip.

Duarte, E., Bayona, G., Jaramillo, C., Parra, M., Romero, I., & Mora, J. (2017). Identificación de los máximos eventos de inundación marina miocenos y su uso en la correlación y análisis de la cuenca de antepaís de los Llanos Orientales, Colombia. Boletín de Geología, 39(1), 19-40. http://dx.doi.org/10.18273/revbol.v39n1-2017001

Guo, G., Díaz, M. A., Paz, F., Smalley, J., & Waninger, E. A. (2005). Rock Typing as an Effective Tool for Permeability and Water-Saturation Modeling: A Case Study in a Clastic Reservoir in the Oriente Basin. SPE Annual Technical Conference and Exhibition. https://doi.org/10.2118/97033-MS

Lis-Śledziona, A. (2019). Petrophysical rock typing and permeability prediction in tight sandstone reservoir. Acta Geophysica, 67, 1895-1911. https://doi.org/10.1007/s11600-019-00348-5

Mora, A. (2015). Petroleum systems of the Eastern Cordillera, foothill basins, and associated Llanos basin: Impacts on the prediction of large scale foreland and foothill petroleum accumulations. AAPG Bulletin, 99(8), 1401-1406. https://doi.org/10.1306/bltnintro032615

Reyes-Harker, A., Ruiz-Valdivieso, C. F., Mora, A., Ramírez-Arias, J. C., Rodríguez, G., De la Parra, F., Caballero, V., Parra, M., Moreno, N., Horton, B., Saylor, J., Silva, A., Valencia, V., Stockli, D., & Blanco, V. (2015). Cenozoic paleogeography of the Andean foreland and retroarc hinterland of Colombia. AAPG Bulletin, 99(8), 1407-1453. https://doi.org/10.1306/06181411110

Rider, M. (2000). The geological interpretation of well logs. Rider-French Consulting Ltd.

Satter, A., Iqbal, G. M., & Buchwatter, J. L. (2008). Practical enhanced reservoir engineering. PennWell Corporation.

Vayssaire, A., Abdallah, H., Hermoza, W., & Negri, E. G. F. (2014). Regional Study and Petroleum System Modeling of the Eastern Llanos Basin. AAPG Search and Discovery.