3D geo‑cellular modeling for Oligocene reservoirs: a marginal field in offshore Vietnam

This study focuses on constructing a 3D geo-cellular model by using well-log data and other geological information to enable a deep investigation of the reservoir characteristics and estimation of the hydrocarbon potential in the clastic reservoir of the marginal field in offshore Vietnam. In this study, Petrel software was adopted for geostatistical modeling. First, a sequential indicator simulation (SIS) was adopted for facies modeling. Next, sequential Gaussian simulation (SGS) and co-kriging approaches were utilized for petrophysical modeling. Furthermore, the results of the petrophysical models were verified by a quality control process before determining the in-place oil for each reservoir in the field. Multiple geological realizations were generated to reduce the geological uncertainty of the model assessment for the facies and porosity model. The most consistent one would then be the best candidate for further evaluation. The porosity distribution ranged from 9 to 22%. The original oil place of clastic reservoirs in the marginal field was 50.28 MMbbl. Ultimately, this research found that the marginal field could be considered a potential candidate for future oil and gas development in offshore Vietnam.


Introduction
The Cuu Long Basin is a well-known source for exploring and producing hydrocarbons in Vietnam (Thanh et al. 2019). However, oil and gas have been extracted from this basin for a long time. As a result, the large reserve oil fields have been depleted in the Cuu Long Basin (Thanh and Sugai 2021). In this basin, the marginal field is a potential candidate for new resources in Vietnam. However, high cost is a significant issue for us when attempting to develop a strategy to explore and develop the marginal field. To address this, 3D geocellular modeling has been adapted to deeply understand possible reserves (Islam et al. 2021). Furthermore, this modeling technique provides the target reservoir's better manner, such as faults, sub-layers, and the amount of oil reserves in the field (Mike 2009).
Moreover, 3D geo-cellular modeling was used to develop a new method by integrating artificial neural networks (ANNs) and geostatistical methods for static and dynamic CO 2 storage in the Nam Vang field, offshore Vietnam (Thanh et al. 2020). Therefore, this study adopts 3D geo-cellular modeling to explore the oil reserve in one of the marginal fields in offshore Vietnam. The geological model of the Lower Oligocene reservoirs of this field was constructed in this study using the interpretation of 3D * Hung Vo Thanh vothanhhung198090@gmail.com;vothanhhung1990@snu.ac.kr seismic acquisition and data from wells drilled in the field, such as L-1X, L-2X, and L-3X. The present geological model aims to gather all available data to build a comprehensive geostatic representation of the Lower Oligocene reservoirs for a marginal field in offshore Vietnam. For the study, the marginal field was divided into two main zones based on two defined reservoir units (Thanh et al. 2020). They are zone E10 from top E10 to base E10 and zone E20 from top E20 to base E20. Generally, structural maps and faults from acquired 3D seismic data are used to generate structural and stratigraphic frameworks. Due to data availability, the stochastic approach was used to model facies and petrophysical properties (Al-khalifa et al. 2007;Cherpeau et al. 2010). Stochastic algorithms use random seed and input data and provide details of the results with differences (Xin et al. 2016). Moreover, stochastic modeling is more complex than deterministic modeling and takes much longer to run. This indicates more aspects of the input data, specifically its variability (Thornton et al. 2018). This means that particular cells will appear in the results that are not driven by the input data and whose location is purely an artifact of the random seed used (Perevertailo et al. 2015). Therefore, the results will have a distribution that is more typical of the real case, although the specific variations are unlikely to match (Ailleres et al. 2019). This can be useful, particularly when taking the model to the simulation stage, as the variability of a property is likely to be just as important as it is the average value (Kamali et al. 2013). The disadvantage is that some crucial aspects of the model may be random. Therefore, a proper uncertainty analysis is performed with several realizations of the same property (Matias et al. 2015).
Furthermore, it is necessary to model the lithology before distributing the petrophysical properties of the 3D model. Several methods can be used for facies modeling. However, SIS has been recognized as the most valuable and flexible method for lithofacies modeling. In addition, SIS can predict the facies characteristics using well-log information to provide the facies distribution in uncore wells. This advantage is supported by stochastic facies modeling (Al-mudhafar 2017a).
Moreover, SIS has also been used to quantify geological uncertainties in modeling the architecture of complex reservoirs (Seifert and Jensen 1999). In addition, the SIS technique was employed to model the facies in sandrich turbidite systems (Jordan and Goggin 1995). Most recently, SIS has been successfully adapted for lithofacies modeling in the Luhais oilfield in southern Iraq (Almudhafar 2021). In a similar study, SIS was considered in the development of an innovative carbonate facies modeling workflow in one of the UAE fields (Aidarbayev et al. 2020). A previous study demonstrated the advantage of SIS for lithofacies modeling. Therefore, in this study, SIS was used to capture the depositional environment for clastic reservoirs in the marginal field regarding the lithofacies modeling.
The organization of this study is as follows: the field and geological description of the marginal field is outlined in Sect. 2, while Sect. 3 highlights the modeling framework for the paper. Section 4 presents the results and discussion of the modeling process. Finally, Sect. 5 presents the conclusions of the study.

Field description
The Cuu Long Basin constitutes a Paleogene rift basin located off the southeast coast of Vietnam, covering approximately 150,000 km 2 . It is separated from the Nam Con Son basin by the Con Son Swell in the southeast (Lee et al. 2001).
In the current tectonic sketch (Fig. 1), the Cuu Long Basin is situated southeast of the intra-continental plate that developed on the Eurasian continental crust. The Cuu Long Basin is a rift type, a subsided trough in the Paleogene that was founded on the cretaceous regressive continental crust. The Cuu Long Basin was filled with Neogene (N1 2 -Q) passive continental margin sediments. During the Cretaceous (J3-K), the region was occupied in the central part of a magmatic arc of a NE-SW trending from the adjacent continental areafrom the Da Lat zone to Hainan Island in China. The Cuu Long Basin basement consists mainly of igneous magmatic arc (volcanic and intrusive) rocks (Dang et al. 2011b).
The tectonic evolution of the Cuu Long Basin can be divided into several main stages (Thanh and Sugai 2021): 1. Late Cretaceous-Eocene: Pre-rift uplift/initial rifting phase. 2. Late Eocene-Oligocene: Main rifting phase/initial ocean floor spreading phase. This led to the development of the main structural features within the basin, following extensional and transtensional deformations. 3. Early-Middle Miocene: Regional subsidence/renewed rifting. A change marked this from fault-controlled subsidence to thermally controlled, high-rate subsidence. 4. Late Miocene: Partial inversion/regional subsidence.
During this stage, the entire area was dominated by compression in combination with the dextral strike-slip fault system in eastern offshore Vietnam, which probably generated basin uplift/partial inversion. 5. Pliocene-Pleistocene: Regional subsidence/renewed rifting. Diverse tectonic activity, from low to moderate amplitude differential uplift, acted across the basin.
Structurally, the Cuu Long rift basin is an elongated depression that is made up of a series of alternating horst, graben, and half-graben structural features arranged along the strike (the main strike is NE-SW) of the basin (Dang et al. 2011a). The location of the study area is shown in Fig. 1.

Facies and depositional environment
The facies and depositional environment analyses along the wells in the marginal field are based on petrographical analysis, core/cutting description, and log forms. The results of the sequence analyses are shown in Fig. 2. The details of the subsequence are summarized below.

Intra E20-basement subsequence
In this subsequence, the facies association is comprised of fine-grained facies-mudstone of lacustrine/floodplain, minor of very fine to fine sand facies of sheet flood and distal bar, and fine-to medium-grained sand facies of small channels. Regarding the depositional environment, this subsequence is a lacustrine, alluvial plain, and braided channel system. Fine-grained facies: mudstone is mainly deposited by suspension. Coarse-grained sediments were dominantly divided from weathered granitoid provenance and rapidly deposited, close to the provenance.

Top E20-intra E20 subsequence
The main facies of this subsequence is dominated by medium-to coarse-grained sand. Sand stacks of different facies, with various thicknesses, are often alternated and separated by thin mud beds or mud drapes. Individual sand bodies within the sand stacks are often separated by thin mud beds or mud drapes. These sand facies combined to form stacked sequences of 5 m up to 50 m thick, with gamma-ray logs representing a subtle cylinder shape. This subsequence was deposited in lacustrine and braided plains.

Top E10-Top E20 subsequence
This subsequence has very fine to fine sand facies, minor of very fine to fine sand facies, and low sand/shale ratio. In addition, fine-grained facies-mudstone-was mainly deposited by the suspension. Coarser-grained sediments were dominantly divided from weathered granitoid provenance and transported over a long distance and slowly deposited in a low-energy environment.

Methodology
The data used for this research are the well-log data (porosity, permeability, and volume of clay), depth fault stick, depth surfaces, well markers, and oil-water contact data; more details will be introduced on each data type later. First, structural modeling was conducted for fault modeling, pillar girding, vertical layering, and contacts (Caumon et al. 2009). Property modeling is the process of filling the cells of the grid with petrophysical properties such as facies, porosity, permeability, and water saturation (Ali et al. 2020). Multiple geological realizations were conducted in the modeling process to reduce the geological uncertainty of the facies model (Al-mudhafar 2018). Thus, petrophysical models have also been generated according to each lithofacies realization (Al-mudhafar 2017b). Finally, the volumetric calculations for all realizations were calculated and ranked to determine the potential of oil reserves in the marginal field of the Cuu Long Basin, offshore Vietnam. Figure 3 illustrates the modeling loop used in this research.

Structural modeling
First, fault modeling is the initial step in the structure model process. Fault sticks were interpreted and imported into a geological package (Petrel Version 2012). Twelve fault stick sets were used for the fault modeling. In general, the Lower Oligocene faults were simple. Almost all fault systems are normal faults and fault orients in the NE-SW direction. Generally, fault modeling generates fault pillars, known as key pillars, which define the slope and shape of the fault. There are up to five so-called shape points along each of these lines to adjust the form of the fault to match the input data perfectly. The key pillars are generated based on input data, such as depth fault sticks and structural maps (Islam et al. 2021).
Pillar girding was the next step in the structure model process. This process generates a 3D grid from the fault model. The grid is represented by pillars (coordinate lines) that define the possible positions for the grid block corner points. Thus, we can define directions along faults and borders to guide the girding process. The result from the pillar girding is the "Skeleton Framework". The skeleton is a grid consisting of a top, middle, and base skeleton grid. The grid has no layers and only a set of pillars with the user given X and Y increments. For example, in the marginal Oligocene model, X and Y were set at 50 m × 50 m, and fault systems were made under force grid cells to be equally spaced along the faults. An illustration of the fault modeling and pillar girding is shown in Fig. 4. The next step in 3D geo-cellular modeling is to create horizons and layering. The make-horizon process step is the first step in defining the vertical layering of the 3D grid in Petrel software version 2012 (Fig. 5). The vertical layering of the 3D grid is defined in two process steps: make the main reservoir layers (from seismic data) and make a vertical resolution (determined by cell thickness or several cell layers).
The depth surface data of the top E10, base E10, top E20, and base E20 in surface type were transferred to Petrel with well tops data. Using the Petrel Make Horizon module, these surfaces were used to create three zones, which combine with the fault system to construct a raw structural model or structural framework as "house of properties". The geological horizon Top E10 is regarded as an erosion, while Base E10 and Top E20 are conformable. Base E20 is regarded as the base of the model.
The layering process is the last step in defining the vertical resolution of a 3D grid. The number of layers in each zone was considered to optimize the number of cells but was thin enough to capture the thinnest sand body. The layer modeling was focused on the main potential productive zones such as Zone E10 (Top E10 to Base E10) and Zone E20 (Top E20 to Base E20); the thickness of layers is 0.5 m. Under the scale of geological nature to be realized and runtime of the multiple simulations required, the grid configuration was designed to optimize the number of cells and computing speed concerning reservoir heterogeneities. In general, the grid size and layering for the fine model are reasonable for representing the geometry and geological heterogeneities of the Oligocene reservoirs. Table 1 describes the information of the modeling parameters.

Property modeling
Property modeling involves filling grid cells with petrophysical properties such as facies, porosity, permeability, and water saturation. Therefore, these processes are dependent on the geometry of the existing grid. When interpolating between data points, the geological modeling process propagates property values along with the grid layers. The second period in building geological models is property modeling (He et al. 2020). The data input for property modeling includes well-log data (facies, porosity, permeability, volume of clay) and seismic trend maps. Property modeling in the marginal field is split into three separate stages (Zhong and Carr 2019): • Scale-up well logs: The process of sampling values from well logs into the grid, ready for use as input to facies modeling and petrophysical modeling. • Data analysis: The process of applying transformations to input data (normally upscaled well logs) identifies trends and defines variograms describing the data. This is then used in the facies and petrophysical modeling to ensure that the same trends appear in the results. • Facies modeling: Interpolation of discrete data, e.g., lithofacies. • Petrophysical modeling: Interpolation of continuous data, for example, porosity, permeability, net to gross, and saturation.

Scale-up well log
The facies was defined on base cases (such as reservoir/ non-reservoir) as codes 1 and 0, respectively. These facies codes were created and applied for three wells based on volume shale and effective porosity cutoff of 0.35 and 0.09, respectively. Figures 6 and 7 illustrate the scale-up results of this study. Figure 6 presents the scale-up for different layer thicknesses for the E10 and E20 zones.  1 3 The layer thickness of 0.45 m is the best for both E10 and E20 reservoirs. Thus, this type of thickness will be utilized for the entire modeling process.
In addition, Fig. 7 illustrates the correlation between zones (E10, E20) and facies classification. As shown in this figure, the reservoir facies are yellow, and the non-facies are represented by the brown color.

Data analysis
The data analysis process (DAP) of applying transformations to input data (normally upscaled well logs) identifies trends and defines variograms describing the data (Gringarten and Deutsch 1999). The DAP can be accessed from the process diagram and performs a detailed property analysis. In addition, depending on whether a property is discrete (e.g., facies) or continuous (e.g., porosity and permeability), we can analyze the facies proportion, thickness, calibration between continuous properties (e.g., sampled seismic), and facies within each zone and create a variogram of the discrete property. The DAP settings defined for discrete properties are saved for the current property and are accessed directly during facies modeling. Moreover, the continuous property can be defined by data transformations and generate variograms in the main, minor, and vertical directions.
In addition, data transformation enables the user to make the data stationary and standard normally distributed, which are requirements of many of the standard geostatistical algorithms (Gringarten and Deutsch 1999). After the scale-up well-log process, the data analysis tool was applied to analyze the facies proportion and thickness and create a variogram model (Gringarten and Deutsch 1999). Owing to constraints on the number of exploration/appraisal wells in the marginal field, the variogram model was created only for the vertical direction. The horizontal variogram was created based on the compressional wave velocity to shear wave velocity (Vp/Vs) variogram maps. The results of the facies and porosity data analysis are shown in Figs. 8, 9, 10, and 11.

Facies model
Facies modeling is a means of distributing discrete facies throughout the model grid . We usually have upscaled well logs with discrete properties into the model grid and possibly defined trends within the reservoir by analyzing this data analysis process. In this research, the "Sequence Indicator Simulation" (SIS) method was applied to build the facies model (Seifert and Jensen 1999). The SIS method allows a stochastic distribution of the property using a predefined histogram. In addition, directional settings such as variograms and extensional trends were also honored. Based on the facies analysis results in Fig. 12, these results provide the facies volume proportion (histogram analysis), facies thickness distribution, and variation.
In the volume proportion, the probability curves were fitted from the original facies proportion at the well. The variogram was also checked carefully for each zone in the primary, minor, and vertical directions. To evaluate the geological uncertainty of the facies model, 20 geological realizations of facies distributions were created for the facies model. Then, the best one was selected to capture reservoir heterogeneity. Figure 13 represents the generation realizations and cross sections for the facies model of the E10 reservoirs. As shown in Fig. 12a, the effects of geological realizations were explored in this study. The Fig. 8 Vp/Vs variogram map is utilized for E10 reservoirs consistent case of facies was chosen to capture the characterization of the geosystem. Uncertainty would occur in geological modeling because the subsurface data are not easy to obtain. Nevertheless, geological modeling is an excellent tool for deeply understanding the geology aspect. For better visualization, Fig. 12b depicts the cross section of the best facies model for E10 reservoir.  Similarly, the generation of facies realizations and cross sections for the E20 zone is shown in Fig. 13. As can be observed in this figure, the amount of sand distribution is higher than that in the E10 zone. Therefore, it seems that the E20 zone has more potential reserves than the E10 zone. In addition, the thickness of E20 is greater than that of the E10 reservoir. However, this is just a facies distribution that cannot determine which area is of interest for future development. Thus, the petrophysical modeling and volume calculation will be discussed in the next section.

Petrophysical modeling
Petrophysical results can be achieved in several ways, such as well-log interpolation (Attia et al. 2015;Abudeif et al. 2016a), petrographic image (Abudeif et al. 2016b), and petrophysical fingerprint techniques (Radwan et al. 2020). These techniques are useful for hydrocarbon-type detection in oilfields (Abudeif et al. 2018). Petrophysical interpretation is also helpful for classifying rock types during the waterflooding process in reservoirs . However, petrophysical modeling is missing from these interesting works. Petrophysical modeling provides the distribution of porosity and permeability properties in the target field. Several steps to populate the model are required for the porosity (PHIE) modeling. First, the PHIE of each wellbore is derived from the log analysis result. Then, the PHIE was scaled up from the well log within the zones using the arithmetic average method. Next, the facies model is used as a filter to define the reservoir/non-reservoir in the PHIE model. The porosity values for all non-reservoir cells were always lower than the porosity cutoff value, whereas the porosity of all reservoir cells cannot be less than the cutoff value. Finally, PHIE modeling is run utilizing the sequential Gaussian simulation (SGS) method, variogram from data analysis, and standard distribution method from the upscaled log distribution (Pyrcz and Deutsch 2014). In addition, SGS could provide an understanding of the geometry and continuity of petrophysical properties that have a direct impact on the reservoir flow behavior (Al-mudhafar 2021).
In this work, PHIE was upscaled bias to the chosen facies distribution, and 20 realizations were run.
Assuming that the porosity at well is the best and the well, the decreasing porosity so that the porosity distribution for the best case is chosen based on this assumption. Figure 14 illustrates the porosity selection after the geological realizations. The best case is selected for the E10 reservoirs. The observation of this process shows that the porosity is consistent after the 17th realization. After this realization, the porosity distribution remains the same as in the 17th realizations. Regarding the E20 reservoirs, the selection of realization is similar to that of the E10 zone. The best-case porosity was achieved at the 20th realization. This scenario is shown to be the most consistent when compared to other realizations. Consequently, a deterministic method was applied to build the permeability model. Based on the porosity-permeability relationship provided from the core data of TL-3X, the permeability model is generated using the following equation: The J function was applied to gather the water saturation spatial distribution within the reservoir grid for water saturation modeling. The J function was obtained from the core data of TL-3X. The functions are as follows: (1) Perm_E10 = 2E − 05 * pow(PHIE, 5.0849) (2) Perm_E20 = 3E − 08 * pow(PHIE, 7.8243)     Figure 15a illustrates the permeability model using Eqs. (1) and (2). The co-kriging technique was used to (8) J_fun_20 = (Pc_20∕26.3) * RQI_20 (9) Sw_10 = 1.5975 * pow(J_fun_10, −0.383) (10) Sw_20 = 2.0544 * pow(J_fun_20, −0.417) Fig. 15 The result modeling is highlighted for permeability (a) and water saturation (b) model for E10 and E20 reservoirs distribute the 3D permeability model. Co-kriging could support the limitation of well data during the petrophysical modeling process (Thanh et al. 2020). This technique is a correlation factor between the primary and secondary data to calculate the distribution of the secondary variable at each point (Esmaeilzadeh et al. 2013). Thus, the co-kriging technique is suitable for the modeling process in this study. The illustration in Fig. 15a indicates that the target reservoir has a good distribution permeability in the marginal field. In addition, Fig. 15b presents the water saturation using the J function for both the E10 and E20 reservoirs.

Model verification
Data and model validations were performed at every significant modeling step. First, when data have been imported to the geological package (Petrel Version 2012), they should be under quality control (QC). Typical methods of data QC are to display them and to check statistics, histograms, etc. The QC process uses a general intersection to view the data in the cross section, and playing through the data set is an additional experience. Then, various cross sections were generated for QC, structural framework, facies modeling, and petrophysical modeling. Second, the skeleton grid was for quality control. The essential steps during QC also include checking for crossing pillars because possible negative cell volumes are generated. Finally, the static method checks the matching between the input and output during the reservoir property modeling processes.
In addition, the synthetic data generated in a global well is a critical parameter for checking model quality. Depth synthetic information was used for comparison with the raw data to ensure no-depth mismatch. Synthetic data were also checked to ensure that the generated parameters matched the raw well data while generating facies and petrophysical realizations.
A statistical value check was also conducted to examine the facies volume proportion and trend to determine the matching status between generated realizations and input well data. Figures 16 represents histograms showing the facies and PHIE distribution of the input and output data  for zones E10 and E20. The visualization analysis from these histograms indicated an excellent model distribution for porosity in both the E10 and E20 reservoirs. In addition, the statistics of the two zones were also conducted for a better comparison. Table 2 shows the statistical analysis of the E10 and E20 reservoirs. Surprisingly, results were obtained between the raw log and 3D model. This indicates an excellent fit between the well-log and geo-cellular modeling approaches. This result demonstrates that our model successfully represents the subsurface conditions for marginal fields in offshore Vietnam. The following section elaborates on the oil reserves for the E10 and E20 reservoirs for future development plans.

Volumetric calculation
The volumetric was computed based on the 3D geological reservoir model (Radwan and Chiarella 2021). Because of the facies definition, cutoff Vcl < = 35% and PHIE > ≥ 9%; therefore, the net to gross (NTG) was taken as 1. The NTG model was generated for volumetric calculation. In addition, the oil saturation (So), oil formation volume factor (Bo), and oil-water contact (OWC) were adapted from the Reserve and Resources Report (RAR) for the volumetric calculation. Volumetric sensitive processes were also carried out to choose the most appropriate OIIP value for the model (closer to the Conventional Volumetric Method shown in the latest RAR). The OIIP is calculated using the following equation (Petrel Software Manual 2012): To evaluate the uncertainty of the geo-cellular models, 200 cases were conducted for each reservoir. Figure 17 illustrates the volumetric calculations for the E10 and E20 reservoirs. The comparative results of the RAR and model calculations are shown in Table 3. There is a slight difference between the report and the model calculation. These results indicate that the proposed model is acceptable for further assessment.

Conclusion
This study applied geo-cellular modeling to evaluate the potential of reservoirs in the marginal field offshore. The following key points can be drawn from this research: OIIP = BRV * NTG * PHIE * (1 − Sw)∕Bo * C • As the maps of reservoirs are not being interpreted, the top and base maps of reservoirs generated from the top E and top basement lead to different volume reservoirs of rock in the Reserve and Resources Report and 3D model. Therefore, the structural model should be interpreted as a map of reservoirs calibrated with seismic attribute cubes. • The facies model helps determine the sand distribution of reservoirs to support petrophysical modeling. This study used 20 facies geological realizations to reduce the geological uncertainty in porosity and permeability modeling. The porosity realization distributes the porosity properties for each facies realization. • The petrophysical model was determined by a relationship adapted from the core analysis. Thus, deterministic modeling was used for permeability, while stochastic modeling was used to ensure the suitability of the target reservoirs. • Sequential Gaussian simulation (SGS) was used for porosity realization to preserve the characteristics of the lithofacies model. In addition, the permeability distribution was adapted to the co-kriging technique to achieve the minimum variance of the estimation error in the modeling process. • The volume of the reservoir rock is one of the critical factors that affect the oil's initial in-place calculation. It should be taken from the 3D model to eliminate the uncertainty. • Geological uncertainty is necessary for 3D geo-cellular modeling in the exploration and production of hydrocarbons. • This model provides the static geological model for future development plans by dynamic simulation The potential oil reserve of the study reservoir was demonstrated using 3D geological modeling. However, uncertainties still exist in the Oligocene reservoir model for the marginal field. Some uncertainties are highlighted below: • Structural maps are still unconfident due to the modeling process. • While not yet studying the depositional conceptual model, the intermediate conceptual model has been referred to as for the assumption of direction, location, and size of geological feature distribution. • Log upscaling and grid upscaling themselves lead to model uncertainty, particularly the hydrocarbon pore volume. Fig. 17. Two hundred cases run for calculation the STOIIP for E10 reservoir (a) and E20 reservoir (b)