Building of Turbiditic Gas Field Dynamic Model with a Simplified 3D Simulation Software

This study provides a novel approach of building 3D simulation model with extremely shorter time needed using Rubis simulation software from Kappa Engineering. The study focused on X Field that is located in a turbiditic setting, mainly consisted of separated channel bodies filled with gas, located in a slope apron or passive continental margin of Mahakam Delta. Methods of the study is quite contradictive with common reservoir simulation where it includes data integration, data quality control, model geometry building, reservoir properties distribution, and is followed by wells definition to build the 3D simulation model. Afterward, the reliability of the structural model was checked by the volume calculation for each segment from GeoX model where all dynamic and static data used in the simulation were checked using history matching data derived from well-testing. In conclusion, simulation was run and X Field will be producing for 23 years with 3 years and 10 months plateau rate. Where the static and dynamic data are already provided, the simulation conducted here was very beneficial during the exploration phase of a gas field where the whole process of modeling and simulation could be done only for 3 to 6 months.


Introduction
X Field is composed of stacked gas reservoirs located 75 kilometre of Borneo Island ( Fig. 1) (Internal report, 2011). It is located within deep-water continental slope setting on the tip of a delta front. The field is currently under early development phase and the first gas production has been expected in the next few years (Internal report, 2011).
This field consist of 2 main segments: X North East (NE) Field and X Main Field where it is separated in the distance of 5 kilometre. This study would be focused in X Main Field where it has more gas in place during the previous exploration phase study. X Main Field was indicated as turbiditic reservoir deposited in the slope of passive continental margin (Internal report, 2011). This type of depositional environment would have the characteristic of intercalation between sand and shale with the shape of channel-like geometry. In the X Main Field, it consisted of mainly 10 separated independent segments that would be modeled and simulated during this study. It will compile the seismic data, core, special core analysis (SCAL), well logs, and well test data to build reliable yet quick 3D simulation model using Rubis software.

Aim
The aim of this study is to use all the data available during the exploration phase of this field to build simplified 3D reservoir model using Rubis software and to see if this model is quick and reliable enough to be compared with simpler model such as material balance model.
In order to build the simplified reservoir, it would be started by building structural model using the seismic interpretation available (Internal report, 2011). Afterward, all well log data will be used to give the value of net to gross (NTG) and average porosity for each reservoir segment (Internal report, 2011). While for the permeability value, porosity versus permeability values derived from core plug measurement from one exploratory well will be distributed in all segments. This might be representative enough for the homogeneous lithology deposited under the same turbiditic condition. SCAL data will provide the extended capillary pressure (Pc) curve to better estimate the value of irreducible water saturation (Internal reservoir study, 2011).
At the end, the ultimate objective of the study was to build a simplified 3D reservoir numerical model using Rubis software in a much faster way than using Petrel and Eclipse. Rubis (Rubis Software (Kappa Engineering), 2015) is somewhere located in between 1D material balance and the massive 3D simulation models. It replaces neither but does give reliable perspective of reservoir performance of a dry gas field. This study would highly benefit the reservoir engineer and geoscientist during the exploration or early development phase by building reliable simulation model in very short amount of time.
2. Geological Framework 2.1 Regional Geology X Main Field is located in Kutai Basin which was formed as a product of the interaction between three main plates: Eurasian plate, Indo-Australian Plate, and Pacific Plate. Located in the top east of Kalimantan/Borneo, Kutai Basin covers about 165.000 sq. km, and is one of the deepest of the tertiary basins in Indonesia which contains up to 12.000 m of Tertiary sediment (Alam et al, 1999). The similarity between the Eocene strata of the Barito, Asem-Asem and Kutai Basins has led many workers to suggest the three basins formed a single Eocene depocentre (e.g. Heryanto, 1993;Mason et al., 1993;Panggabean, 1991;Pieters et al., 1987;van Bemmelen, 1949;van de Weerd and Armin, 1992 Kalimantan Mega Ba e.g. Heryanto et al., 1996). The Kutai Basin became separated in the Early Oligocene due to down-throw on the northern side of the Paternoster Fault (Witts et al, 2012).
X Main Field was deposited in a continental slope where it was characterized by numerous stacked gas charged channels in the Pliocene succession (internal report, 2011). High energy deposition of turbidite complex produced sand-shale intercalation channel shape reservoir with Bouma sequence as a common sedimentological feature (Bouma, 1962). The channels can be up to several kilometres wide and hundreds of kilometers long, and provide the transport pathways for large quantities of sediment, nutrients and carbon Talling et al, 2007). X Main Field is located in a confined channel geometry, where the X NE Field is located in the downstream part of a submarine slope where lateral distribution is more extensive in this particular field ( Fig.1) (internal report, 2011). This field was composed of 10 separated segments that was based on its depositional environment. Each segment was filled by a single deep-water turbidite channel. It would be modelled as 10 different layers surrounded by shale (zero flow capacity lithology) based on its natural behaviour. These 10 layers were named as Segment B2, A2, A, C, D, E, F, G, I1, and I2 from the shallowest to the deepest. Reservoir would have channel geometry with around 762 meter wide and 6,7 km long based on the seismic interpretation. With the intercalation of sand-shale lithological features inside the reservoir caused by the deep water turbidity current, it will make a challenging reservoir to be modelled and simulated either by geologist or reservoir engineer.
Theoretically, turbidity currents tend to make their way through channels and bound themselves by many natural levees (Salaheldin et al, 2000).These channels open paths for the turbidity currents and directing sediments to the lower regions of submarine fans. With the gradual decrease of both slope angle and gravitational potential energy, the rate of sediment deposition decreased, forming a concave profile. Therefore there is expected loss of momentum to bring particle of sand onto far distances. The fine materials such as silt and clay are essential factors in turbidity currents to transport sand-sized particle (Reading and Richards, 1994). In the presence of the fine particles, turbidity currents able to keep enough momentum to produce suspension so the sands are able to travel along the channels (Mutti and Lucchi, 1972). Yet in the reservoir geometry modelling process, the channel would have the shape of elongated separated lobes with a concave base for each lobe.

Method
This study would compile seismic data, cores, SCAL, well log, and well test data in order to build 3D simulation model. Based on the seismic, log, core, and well testing data on Well KC-2, the reservoir model would be built based on this following workflow: 1. Geometry building, which is integration of depth structure maps (either top or bottom of an interval), isochore maps, and reservoir segments extensions into a structural model. Nevertheless, there were several assumptions during the building of this simplified gas reservoir. These assumptions were extracted based on the software capability to model the reservoir and due to the data limitation (Table 1).

Availability Data
Well data available in the X Main Field that is summarized in Table 2. KC-1, KC-2 DIR, and KC-3 DIR were the exploratory wells and the other wells were development wells that will be used for the gas production later on (Internal report, 2011). Interpreted seismic lines crossing all reservoir segments were also available during this study.

Seismic Data
Seismic interpretation has been done in order to understand the shape of the reservoir geometry, connectivity between the reservoirs, reservoir distribution, and fluid contacts. As mentioned before, the reservoirs are separated by shale-dominated zones in between. Based on amplitude anomalies on seismic line from North to South, seismic interpretation produced non-continuous top and bottom boundary of each separated reservoir segment. The seismic section was then used to create depth structure map of the whole field (Fig. 3). Submarine canyon and channel geometries with their confined boundaries could be clearly seen from the map.  This seismic interpretation would then be used to build basic structural model of the field. First of all, all the interpreted lines will be a guide on building the depth structural maps for each segment by running a time to depth conversion using check shots data available in the wells. Afterward, isochore map was built based on the true vertical depth (TVD) subtraction of the top to bottom structural maps (Fig 4a). This thickness map would be used to distribute the thickness of the reservoir in each segment (Fig. 4b).

Core and Well Log Data
Core data is available in Well KC-2 DIR in segment F where core description, routine core analysis (RCAL) and SCAL were available. These data would be used as the reference to populate the reservoir dynamic properties in other segments.
In general, segment F consisted of mostly intercalation of sand and shale with fining upward sequence as the main characteristic of turbiditic flow product (Fig. 5). Overall, 3 main lithofacies were recognized. They are described in the following lines from finer to coarser grain size. Facies 1 consisted of very thin graded beds of silt, silty clay, silty sand, and coaly flakes. The presence of grading and lack of traction features indicate that these are pure fallout from a sediment-water suspension or may be a deposited pelagic sediment, representing the surrounding reservoir lithology (zero flow capacity) or outside the main channel bodies. Brownish siderite (FeCO3) nodules and bands are often present. Regarding the Bouma Sequence, it may be included as Td or Te of the sequence.
Facies 2 consisted of intercalation of gradedlaminated fine sands and graded siltstones with almost equal proportions. Sand Layers show pure traction features (traction ripples). Thicker sand beds are often characterized by the occurrence of laminated intervals that are very rich of coaly flakes. The laminated intervals are describing cycles of current velocity increase and decrease. This was associated to Tc of Bouma Sequence.
Facies 3: mostly medium and thick-bedded layers of medium to coarse grain sand, containing intervals rich in disorganized coaly flakes and small mudstone clasts. Beds display a sharp, locally erosive base and wavy tops. These beds are the products of moderately concentrated sediment gravity flows. These layers represent the sedimentation within the channel bodies and might be considered as Ta of Bouma Sequence.
As a consequence, during the modelling, facies 1 would be considered as a non-reservoir facies. Facies 2 and 3 would be distributed inside the channel reservoir bodies where these facies have more flow capacity compared to facies 1. While well log characteristics of these 3 lithofacies were noticeable from the value of gamma ray (GR), resistivity (Res), sonic, neutron and density log (NPHI-RHOB) (Fig. 6).

Special Core Analysis
SCAL data was available in Well KC-2 DIR with 8 samples used for capillary pressure measurement and relative permeability measurement. Pc curve was obtained using air-brine in centrifuge at overburden pressure, while Kr curve was obtained using brine acting as the reservoir liquid (wetting phase) and nitrogen acting as the non-wetting phase. X-ray was also used in order to monitor the saturation changes during steady state testing. All core plugs from SCAL were used to obtain dynamic rock typing. Dynamic rock typing was defined using the Leverett J-Function (Leverett, 1940). This function had been used for correlating capillary pressure data for rocks with similar pore types and wettability. This technique was appropriate for correlating capillary pressure data from samples having a similar pore size distribution; as a consequence, this was not appropriate where there was heterogeneity of rock types within the reservoir (Eqn. 1). Plotting the Leverett J-Function resulted in 3 dynamic rock types that are easily distinguished (Fig.  7).
(1) Where: Pc : capillary pressure  : interfacial tension of the fluid (72 dynes/cm was used)  : contact angle (0 degree for perfectly water-wet system Well testing had been done in one exploration well, testing particularly segment F. This deviated well was penetrating more or less in the central part of high amplitude channel shape reservoir (Fig. 9). Thus well test interpretation could be expected to show two parallel boundaries. This segment consists of mostly sandstone with average porosity of 0.28 and wide range of permeability from 85 to 1500 mD based on core plugs measurement in ambient stress condition (800 psig). This interval was perforated from 7272-7377 ft mD and 7386-7402 ft mD (Fig. 10). From the well log data, this interval consists of 2 different sandstone intervals separated by one shale layer in between. Reservoir properties are very good and lead to 70 feet of net pay. The interval is dry gas bearing. It has 0.6 specific gravity and consists mostly of methane with no H2S. MDT data was also available in this interval (Fig.   11.). The MDT data showed that segment F consists only of 1 connected segment, despite the presence of 2 sand packages on the logs. This is coherent with the geological environment of the reservoir where in a turbiditic setting, a reservoir body that already settled before could be cut by the next sedimentation event on top of it. This would form stacked-channel reservoir with shaly thin shale layer in between the sand bodies. Therefore it is possible to have vertically connected sand bodies like in this segment with non-continuous shale layer in between. This data then will strongly drive the well test interpretation. Fig. 11. MDT data from segment F Main build up analysis showed that several models could be used for matching. There are homogeneous single layer model with parallel boundaries, dual permeability model with parallel boundaries, and even 2 layers with parallel boundaries. But the most reliable model that should be used was homogeneous single layer model. The data plot, log-log plot, and semi-log plot were well matched especially in late time (after 1 hour) (Fig. 12).    Well test interpretation resulted in early time small wellbore storage (due to downhole shutting), infinite acting radial flow and homogeneous reservoir during middle time and two parallel boundaries during late time. This was coherent with overall geological model and MDT data which indicated single layer reservoir. Extrapolated final build-up pressure (P*) was plotted with MDT pressure data, showing a 4 psia difference, which is acceptable (Fig. 13). Based on the detail result of well test interpretation, the permeability value had a mismatch between values from well test and from core plugs measurement (Table  3). Arithmetic average permeability obtained from the core plugs was 395 mD. This might be due to cores not being fully representative of the whole sand package and therefore by-passing some lower quality intervals. This might be also due to different (P, T) laboratory conditions compared to reservoir conditions. Nevertheless, this difference was still within the acceptable range (same degree of magnitude), thus permeability value from well testing would be applied for the model.
As discussed before, the interpretation strongly respected the geological setting of the field where the 70 feet net pay corresponded to a stacked channel body. Even though there was 9 feet shale layer in between, this shale layer was probably not continuous enough to be an efficient vertical flow barrier between those two sand bodies. This might be a typical noncontinuous shale layer in turbiditic environment representing fall-off or pelagic sediment. There were 2 boundaries interpreted to be located 500 feet and 2000 feet from the well, which was consistent with the seismic interpretation presented before. These boundaries were interpreted as sealing boundaries, with typical 1/2 slopes in derivative plot. In the turbiditic environment mostly consisting of intercalations between sand and shale, formation heterogeneity was represented by 0.01 vertical anisotropy with overall k*h of 8000 mD.ft. This well test interpretation would then be used to support the reservoir properties and behavior of the 3D reservoir model.

Reservoir Simulation
As mentioned earlier regarding the workflow of building the reservoir simulation (see Fig. 2), reservoir simulation would be started by building the structural model or geometry building, followed by reservoir properties distribution across the field, placing the wells, fine tuning, pressure and saturation initialization, and running the simulation which will create a production forecast.
Step by step details on building the reservoir model will be discussed in the following sections.

Geometry Building
During this step, the reservoir lateral extension and the number of geological layers were defined. The X Main Field consisted of 10 separated segments with specific reservoir properties in each segment. Using the data from well log correlation and seismic amplitude, layers and the extension of reservoir segments could be defined.
Lateral extension of the reservoir segment were defined by the depth structure map based on the seismic interpretation (Fig. 14). In Rubis software, simplified segment polygons are necessary in order to optimize CPU time without jeopardizing model accuracy in terms of in place volumes. Where quality control (QC) of in place volume will better guide the robustness of the reservoir model.
Regarding the 3D extension of the reservoir, depth structure map from Petrel was imported to the model and the depth contour points were determined manually based on the depth structure map available. Beside the reservoir geometry or boundaries laterally, depth structure map could be also extracted during the contouring step (Fig. 15a). The depth structure map should be traced in Rubis software for some major contour points (i.e. every 250 feet). Afterward, Rubis would interpolate points in between using several geostatistical methods that could be independently chosen. In this case, simple linear interpolation was chosen to create the top segment or boundary of the reservoir layer. As long as the interpolated map in Rubis was still in accordance with detail depth structure map from the seismic interpretation, any methods to interpolate the contour points would remain considerable.
The next step after the contouring was to define the reservoir base layer. At that time, the isochore maps were available to be extracted as the thickness for each reservoir segment. It was started by using the isochore maps based on seismic interpretation to be traced in every point in order to mimic the thickness heterogeneity inside the segment (Fig. 15b). This step is

KC-2 DIR
quite similar with the previous step, the only difference was only the data used during this step.

Reservoir Properties
During this step, properties in each reservoir segment was defined. In Rubis software, heterogeneity inside a reservoir segment could not be managed to be heterogeneous. Therefore, simplification of reservoir properties for each reservoir segment would be done. One value of permeability, porosity, NTG, saturation, and Kr-Pc derived from well data penetrated each reservoir segment would be spread inside a reservoir segment. Ideally, populating the reservoir properties inside the reservoir segment should follow the sedimentological concepts of deep-water turbidite system with the bunch of data. But due to the limitation of the software capability and well data, permeability property would be mitigated by using PHI-K law obtained from core plug measurement in well KC-2 DIR. The result of the reservoir properties in each reservoir segment are homogeneous for NTG and porosity, but for the permeability distribution, it follows the PHI-K law obtained in well KC-2 DIR (Fig.  16). Fig. 16. Distribution of each reservoir property in the model that was derived from well data penetrating each reservoir segment (a) Net to Gross property was homogenous in each reservoir segment (b) porosity distribution across reservoir segments were derived from well log data that was simplified having a same value for each reservoir segment (c) permeability distribution across reservoir segments followed the PHI-K relationship from sandstone interval in KC-2 DIR well.

Wells Definition
During this step, well trajectory, perforations, completion, and well control would be defined. This field consisted of 3 exploration wells and 5 development wells. All the development wells would be completed as producing wells while the exploratory wells would be used for well test history matching.
Trajectory for each deviated well was defined with the deviation, kick-off point, well length, and total depth from the drilling information. The actual trajectory was using coordinate system, while in the model, well tracjectory and position were simplified using the actual position found in Petrel model. Unfortunately, Rubis software encountered several errors where the deviated wells were created. The errors mainly due to the well that was too close to the segment boundary or too close to another well. In order to overcome this problem, several 3 deviated wells were represented as vertical wells.

Grid Building
After the grid was designed and filled with both static and dynamic properties, Rubis would automatically build the grid with minimum number of cells. However if needed, the user has the freedom to define the number of cells and the geometry of the cells manually both vertically or laterally. For the sake of running simulation time, 2 vertical cells were defined for each reservoir segment. Total cell worked in this simulation consisted of 23.627 active cells with hexagonal geometry having 100 meter length of each side of polygon.
Regarding the near wellbore area, the most important parameter was the upscaling factor with a range between 0-100% value. This parameter would control the grid refinement around the well during the simulation. With no upscaling (0%), the cells surrounding the well will be kept with very fine size. While 100% scaling would make the number cells around the wells to be highly reduced allowing for a faster simulation run. This well module option would be particularly useful in case of water coning or very heterogeneous properties.

Model Quality Control
In order to be able to build reliable model, quality control was done before running the simulation. Quality control consisted of both dynamic and static parameter. Dynamic parameter consisted of history matching in segment F using KC-2 DIR well test data while the static parameter would comprise of comparison of segment in place volume between Rubis model and GeoX model.
The history matching was done by inputting the Rate and pressure data from well test analysis and interpretation in Saphir was imported directly to Rubis. Simulated pressure in the model matched perfectly with the well testing pressure. There was 4 psia difference during final flow and 2 psia difference during the main buildup stage which were acceptable.
Afterward, in place volume for each segment was compared with P50 in place volume estimated in GeoX model. Some adjusments in geometrical modelling were needed in order to get a good match between in place volume from GeoX model and Rubis model. At the end, average total volume in place difference between GeoX and Rubis model was 3.47%.

Initialization and Simulation Setting
In this step, intial state was defined with reference depth, reference initial pressure, water gas level (WGL), and Kr-Pc relationship for each reservoir segment. Those data were gathered based on SCAL result in KC-2 DIR Well and would be populated throughout all reservoir segments. Rubis software would calculate the saturation distribution in each cell based on the capillary pressure, relative permeability, and WGL given.
During the initialization, Rubis will use Pc values available to get Pc curve as a function of height above free water level in equilibrium using hydrostatic pressure equation. Afterward, the relationship between Pc and Sw will be approached using J-function where it takes into account the reservoir characteristic (porosity, permeability, and wettability). Based on this Pc-Sw relationship, Rubis would have a height as a function of water saturation and would be spread across the reservoir segments.
Maximum gas saturation above FWL would be 0.72 in all segments, this value was based on the irreducible water saturation value derived from SCAL (0.28) with very thin transition zone. In a gas/water system, large difference between gas and water density quickly leads to high capillary pressure, limiting the extension of the transition zone. The simulation would be run from 2017 to December 2032. All the perforations in the development wells would be opened with the surface gas rate target of 57.600 Mscf/d and 35 bars minimum surface pressure. When simulation started, 3 scenarios were made in order to compare the result from 2 MBAL cases and base case from Rubis model (Table 5).

Results and Discussion
With the same field target rate, those models produced different results. The shortest plateau rate was obtained by the base case model, the MBAL case could provide longest plateau rate with the highest recovery factor, while MBAL trans case produced in between (Fig.17). In the base case, the field will produce 815 Bscf of gas with the recovery factor of 66%. The gas peak rate of 288 MMscf/d can be retained up to 3 years and 10 months ( Table 6).
The preliminary study of the gathered data from seismic, well, core, and well testing had been used in order to build simplified 3D simulation model from Kucing Field. The result of the study had been presented as follow: Fig. 17. X Main Field gas production forecast Using Rubis as 3D reservoir simulation software has several advantages and disadvantages that one should be aware of. Several advantages using Rubis could be explained below:  User friendly Rubis is very easy to use where everything has a simple figure to make the user understands very clearly. Simulation steps are also well-organized thus the one can follow from PVT, field geometry, to running the simulation very easily.  Time efficiency This software was not intended to build very detail model nor complex model. As well as simplification of the model and integration between static and dynamic model in one platform could save time efficiently. Building this X Main Field model took more or less 3 months, while running the full field simulation only took 30 minutes. Comparing with other software, Rubis is the efficient one in term of time consuming.  GRDECL input As another benefit, it is possible to load GRDECL that may contain net to gross, porosity, or permeability distribution that may have been built by another software. While some disadvantages from this software would be explained below:  Well control limitation Not like in other 3D simulator software, Rubis could not have any field control. That means the field should be control per well. For example in this case, where the field should produce with 288 MMscf/d plateau rate, it could be defined as each well in this field will produce with 57.6 MMscf/d surface gas rates where there are 5 wells available in X Main Field. This could lead into some difficulties for reservoir engineer to target and control the field production. The other effect of this limitation (particularly in X Main Field) was the rough decreasing rate of X Main Field.  Software instability Since this software was not intended to model a complex field, building X Main Field where there were several separated segments with more than 50 regions were created could encounter some errors. Also some deviated wells could not be assigned perfectly in this model due to the complexity of the field. In order to anticipate these errors and be able to successfully build 3D simulation model in Rubis, based on this case, several recommendations could be proposed: 1. Rubis is highly recommended to model a simple field where the segments were not separated like in X Main Field or other turbiditic channels field. 2. Rubis may also be useful to model the field that has no deviated well. This due to the design of the well that could not be perfectly represented in Rubis. 3. To check the reliability of the model, several quality controls could be done. Quality control might consist of history matching, compared volumetric with previously available in-house volume, well tops quality control, and when it is available, Rubis model should be compared with the more detail Eclipse model.

Conclusion
1. Well log data, core description, and SCAL result led to the determination of lithofacies, petrofacies, and electrofacies, while further analysis of SCAL result " Pc curve and visual inspection of SCAL result to help the determination of the dynamic rock typing. 2. From well test interpretation, Segment F behaved as a single homogeneous layer with 2 parallel boundaries indicating a channel reservoir body with k=114 mD. This result then was used as a reference to populate the reservoir properties in other segments. 3. The quality control had been done by doing history matching based on short historical data from well testing in Segment F. It gave only 2 psia differences between the model and historical data. While structural model quality control had been done as well by comparing the well tops and segment volume that gave 3.47% volume discrepancies. 4. The production forecast indicated that the field would produce in plateau rate of 288 MMscf/d for 3 years and 10 months. At the end, as a new approach in building 3D simulation model, Rubis could be a powerful tool to build reliable yet simple 3D simulation model with very short time needed to build one (more or less 3 months) compare to build a complex Eclipse simulation model where it will take approximately a year.