Improved natural state simulation of Arjuno-Welirang Geothermal field, East Java, Indonesia

Arjuno-Welirang Geothermal Field is a green field located across Mojokerto Regency, Malang Regency, Pasuruan Regency and Batu City in East Java Province, Indonesia. Geothermal system in Arjuno-Welirang Geothermal Field is estimated to be liquid dominated and has a reservoir temperature of about 260°C based on CO2 gas geothermometry method. Several manifestations are found at the surface; consist of hot springs, fumaroles and alterations. The field is now managed by Geo Dipa Energi. Geological, geophysical, and geochemical studies have been conducted in the field but exploration drilling has not been carried out yet until recently. For reservoir characterization and resource estimation purposes, several conceptual models have been developed based on available geoscience data. Previously, a natural state modelling attempted to match existing conceptual models has been carried out by Wardana. An entirely new natural state numerical model is made in this study to improve the result from previous model by using smaller grid blocks (40% finer). In return, the result shows more accurate temperature matching of surface manifestations than the previous model. The numerical model is developed using TOUGH2 simulator to obtain natural state condition. The objectives of this study are to develop an improved natural state numerical model; to validate temperature of the manifestations in the model with measured manifestation data; to characterize geothermal reservoir including area, thickness, pressure, temperature and permeability distribution; and to determine resource potential using probabilistic method. All data used in this paper are based on published papers.


Introduction
Southern part of East Java is passed by the volcanic route, which is known as the ring of fire which consists of active volcanoes. According to Abadi Poernomo, Chairman of INAGA (Indonesia Geothermal Association) in his speech at Institut Teknologi Bandung in 2016, East Java has 1,012 MW of total possible geothermal reserve. Arjuno-Welirang is an integrated investigation area located across Mojokerto Regency, Malang Regency, Pasuruan Regency and Batu City in East Java Province, Indonesia ( Figure 1). This working area was handed by Government of Indonesia over Geo Dipa Energi in early 2017 to be developed further. This working area is now still in the preliminary survey phase and has not been entered exploration phase.
Geological, geophysical, and geochemical investigation have been carried out in this area, which produce result of hypothetical resource estimation around 265 MWe [4]. A natural state modelling attempt has been done by [10] based on an improved conceptual model, which is developed by synthesizing and enhancing the existing conceptual models. Three scenarios were simulated based on  Numerical reservoir simulation is used to imitate flow process with real parameter of reservoir condition using Darcy's law and heat-mass balance. Natural state matching is the first step in calibrating the model [3]. It is matched with relevant data such as the natural temperatures and pressures and the amount of surface discharge as both heat and mass. The reservoir model is constructed with mass and heat input at the bottom, infiltration from the surface, and leaks at sites of surface or subsurface discharge. The structure of the model should be as simple as possible but still having the mechanisms that affect reservoir processes. The model later can be used to predict the performance of the real reservoir over time. This simulation study particularly incorporates finer grid blocks than the previous study done by [10] with some improvements in temperature matching of surface manifestations.

Objectives and limitations
The objectives of this study are: • To develop natural state numerical model of Arjuno-Welirang Geothermal Field based on conceptual model generated from geological, geophysical, and geochemical data using TOUGH2 simulator. • To improve the result of previous numerical model by [10].
• To validate temperature of the manifestations in the model with measured manifestation data of hot springs and fumaroles. • To characterize geothermal reservoir including area, thickness, pressure, temperature, and permeability distribution from the model. • To determine geothermal energy resource of Arjuno-Welirang Geothermal Field numerical model using probabilistic method.
This study also has several limitations: • Exploration wells has not been drilled.
• Heat loss and discharge rate of manifestations data is unavailable.
• The use of EOS1 to model heat and mass flow in the model which only covers pressure, temperature, and gas saturation.
Manifestations in Arjuno-Welirang are distributed over 5 locations. Three of them are hot water springs (Padusan, Coban, and Cangar), while the rest are fumaroles and alteration located at the peak of the mountain. These manifestation data below were collected by [4].
• Padusan Hot Spring. Hot water emerges from pyroclastic lava stones and andesite lava chunks from Mount Welirang. The water is tasteless, colorless, and it contains iron oxide. As written above, all of the hot springs are situated at lower elevations and have lower temperature than the fumaroles. These hot springs are classified as bicarbonate water based on Cl-SO4-HCO3 ternary diagram. They are also classified as immature water zone based on Na-K-Mg diagram. The classification indicates an influence of surface water to the hot liquid from geothermal system resulting hot springs at the surface. The hot springs are also estimated to be the outflow of the geothermal system [4]. Based on those evidences, the geothermal system is estimated to be a single phase liquid.
Reservoir temperature of the Arjuno-Welirang geothermal system is determined using three geothermometry calculations method: SiO 2 , Na/K, and CO 2 , which result 176 o C, 313 o C, and 260 o C respectively. CO2 method shows the most reasonable result. The high concentration of SiO2 at Padusan Hot Spring indicates a correlation between the fumaroles and the hot springs, hence CO2 gas geothermometry is the most suitable choice which gives a temperature of 260 o C [4].
In 2010, the Center for Geological Resources of the Ministry of Energy and Mineral Resources of Indonesia has carried out magnetotelluric (MT) data acquisition in the northern and north-western part of Arjuno-Welirang prospect area. One of the characteristics of a geothermal system found in the resistivity model was the presence of conductive layer. The geometry of the resistivity structure has an up-doming structure for the base of conductive layer (BOC) centred below Mount Welirang.

Conceptual models
Conceptual model is the first guide to the numerical model. It is a comprehensive analysis from geological, geophysical, and geochemical data. Geological data gives the raw structure of the model. Geophysical data provides the reservoir boundary estimation, through geophysical anomaly. While geochemistry data estimates upflow and outflow zones across the area. The latest conceptual model was made by [10] by synthesizing and enhancing the existing models by [4], [9], and [2]. Two cross section of the conceptual model were made to show the outflow to the west and northwest of Mount Welirang. Figure 2 depicts the A-B cross section showing estimated clay cap and reservoir position and iso-temperature curves.
The upflow zone lies beneath the summit of Mount Welirang. Due to the absence of thermal manifestations in the east and south of Mount Welirang, it is estimated that there is no outflow towards those directions. The possible reason is the influence of the low permeability lava body from Mount Arjuno and the elevation is higher than Mount Welirang.
As described above, reservoir and clay cap of Arjuno-Welirang geothermal system is interpreted from the MT 3D inversion data and the magnetic interpretation. The temperature profile is based on manifestation temperatures and a reservoir temperature of 260 o C, which is estimated from CO2 gas geothermometry method. Unfortunately, drilling data is not available to confirm the temperature profile.

Methodology
The method in developing natural state reservoir model consists of several steps. It starts with collecting geoscience data (geological, geophysical, geothermometrical, and geochemical data). Then conceptual model is built according to the interpretation of these data. A computer numerical model based on the conceptual model is made using TOUGH2 simulator to simulate fluid and heat flow process of the reservoir. This numerical model consists of grid blocks in which each block is assigned with material properties. Top, bottom and boundary conditions of the model are set based on the conceptual model. The model is run to get the result of the model. The material properties, top, bottom, and boundary conditions are calibrated iteratively to obtain the best match between the model and the conceptual model which indicates that the model successfully represents the natural state condition of the reservoir. The workflow of this simulation is shown in Figure 3.  Figure 4, which is smaller than the previous model made by [10], in order to emphasize more on the upflow and manifestations area. This three-dimensional model covers total area of 10.8 km × 14.4 km which equals to 155.52 km 2 . The lowest elevation of this model is -1,500 masl and as described above the surface elevation varies from 636 masl to 3,208 masl, which makes the total thickness is within range of 2,136 m to 4,708 m.

Grid model
In order to simulate flow in reservoir better, the grid system around the reservoir area is made finer than the other area (200 m × 200 m). Grid sizes in other areas are 200 m × 800 m and 800 m × 800 m. Figure  4 shows the upper view of the model. It can be seen that the total lateral grid blocks in x-axis direction is 27 and 29 in y-axis direction. These lateral grid sizes are smaller than the one in the previous model which is very coarse (1 km × 1 km). A detailed information of lateral gridding system of this numerical model is described in Table 1.

Wells
Five dummy wells are added in the model to match the reservoir temperature from model to geothermometry estimate (around 260 o C) and match the temperature slope to the iso-temperature line in the conceptual model made by [10]. Dummy well ARW-1 is located right at the centre of upflow area. Dummy wells ARW-2 and ARW-3 are embedded at the west of upflow area following the cross section A-B of the conceptual model. The other dummy wells, ARW-4 and ARW-5 are installed at the northwest outside the upflow area following the direction of Padusan Fault. Coordinates of the dummy wells are presented in Figure 6.

Material properties
Materials are used to represent rock properties in each zone. The properties are porosity, permeability to x, y and z direction, density, wet heat conductivity and specific heat. Every single change in property will affect the simulation result, but permeability has the greatest influence. These properties later should be calibrated to meet the best simulation result.

Initial and Boundary Conditions
Initial condition is used to define the initial state (pressure and temperature) of every grid block in the numerical model. The pressure and temperature are defined by equations (5) and (6) The top boundary is defined by ambient atmosphere temperature of 25 o C and pressure of 101,325 Pa. This condition is kept constant during the simulation time by applying volume factor of 1×10 20 . Annual rain fall of 2,000 mm/year (BMKG, 2010) and an infiltration rate of 10% are represented by water injection into the top of the model.
The side boundary is assumed to be far from the model. Its permeability is set low (0.001 to 0.01 mD) to ensure no flow of heat and mass coming into or going out of the system during natural state condition.
The bottom boundary is embedded with a heat source with volume factor of 1×10 20 , pressure 3.5×10 7 Pa and 3.1×10 7 Pa at the lowest and the second lowest layer respectively, and temperature of 300 o C. The constant flux of heat-in of 80 mW/m 2 and injection with enthalpy of 1.4×10 6 J/kg are also applied at the heat source [10]. It is also assumed that the heat source has a constant flux of water rate injection of 2.5×10 -5 kg/s-m 2 which is assumed based on PetraSim User Manual (2007).

Pressure and temperature
After some calibration efforts, the numerical model finally matches the conceptual model. The numerical model was run until it reaches steady state condition with simulation time of 5.24×10 15 seconds (166,175,799 years). As described above, due to the absence of well data, the matching process is only feasible for these parameters: • Temperatures in manifestations at the surface of the model to the real measured data.
• Reservoir temperature in the model to the estimation of reservoir temperature using CO 2 gas geothermometry method. • Iso-temperature lines in the model to the ones from the conceptual model of [10].
• Heat and mass flow in the model to the conceptual model. There is no measured pressure data, thus pressure matching cannot be done. In the reservoir, the heat spreads from the center of upflow zone beneath the peak of Mount Welirang towards northwest and west direction. The heat is channeled through Padusan Fault to the northwest and Ledug Fault to the west.
The well is right above the heat source (area with temperature of 300 o C). Dummy well ARW-1 is located exactly at the center of the upflow zone, which is also the location of fumaroles. As seen in Figure 7, surface temperature of dummy well ARW-1 is 160.65 o C which is 23 o C higher than measured fumaroles temperature data (94.1 -137.5 o C). Reservoir under dummy well ARW-1 is ranging within 1300 masl until -900 masl. The temperature is around 240 -260 o C, matched to the result of CO2 geothermometry estimate (260 o C).
The temperature profile in reservoir layer under dummy well ARW-2 ( Figure 8) is lower than the estimated temperature due to the location of the well which is at the outside of upflow zone, where the temperature has decreased.  Dummy well ARW-3 is located right near Cangar Hot Spring at the west of outflow area. The temperature profile (Figure 9) shows the hot spring temperature at the surface (50 o C) which is close to the survey temperature (54 o C & 48.3 o C). The subsurface temperature is much lower than the one in dummy well ARW-2 due to the farther distance to the reservoir at the upflow zone.
Dummy well ARW-4 is located outside the upflow area, thus the temperature must be lower than the reservoir temperature. Simulation shows that temperature beneath ARW-4 does not exceed 200 o C. The temperature rises from the surface down to the elevation of 0 masl due to the presence of Padusan Fault which channels heat from reservoir. Then the temperature drops again at elevation below the fault.   Surface temperature distribution is presented in Figure 10, also showing manifestation temperatures. Table 2 compares manifestation temperatures from the model and measured data. Iso-temperature line from the numerical model is matched to conceptual model of [10]. Figure 11 shows the heat source temperature at the bottom of the geothermal system is around 300 o C. Reservoir temperature is ranging within 200 -300 o C. The 100 o C iso-temperature line approaches the surface at the peak of Mount Welirang producing outflow of fumaroles.

Heat and mass flow
Sliced numerical model is used to observe the heat and mass flow at several places. Figure 12 describes the heat and gas flow vertically from the heat source to the fumarole at the peak of Mount Welirang. The heat and mass successfully move through the cap rock by the support of a higher permeability conduit. As fluid with high temperature flows through the conduit, it meets the lower pressure near the peak of Mount Welirang, thus gas is formed. Gas saturation distribution of the model is shown in Figure  13.
Simulation result shows that heat and mass of Cangar Hot Spring are originated from Ledug Fault which connects the reservoir to the western part of the field. Coban and Padusan Hot Spring are the products of Cangar and Padusan Faults respectively.
Simulation result is able to show the high permeability of Padusan Fault which channels the heat and mass from the reservoir to northwest direction. It also shows the high temperature area decreases respectively from heat source, reservoir, up to cap rock. In the heat source and reservoir layer, there are channeling heat from the upflow zone to the Northwest and West direction through the high permeability fault. The reservoir is passed by Padusan (Southeast to Northwest) and Ledug Faults (East to West) which bring the heat from the reservoir to the outflows in Padusan and Cangar Hot Springs respectively as temperatures decrease along the path of the faults. There is also Cangar Fault (Southwest to Northeast) which also connects Ledug Fault to Coban Hot Spring. Permeability range of the reservoir is 20 to 80 mD.
The reservoir is assumed as single porosity system to simplify the model. It only contains one phase of liquid, thus the geothermal system is one phase reservoir where fluid temperature is less than the boiling point, except at the peak of Mount Welirang where gas is flowing out to the fumaroles as seen in Figure 13.

Model sensitivity
Monte Carlo simulation is simple probabilistic method used to calculate the production ability of a reservoir during. It utilizes heat stored method to determine generated electricity, which requires several parameters: area, thickness, porosity, rock density, rock heat capacity, initial and final water saturation, initial and final temperature, recovery factor, electricity conversion factor, and the contract lifetime. These data shown in Table 3 are collected from reservoir characterization and assumptions. Monte Carlo uses random numbers within the parameter range. This simulation uses 60,000 random numbers. Furthermore, sensitivity analysis is also carried out using deterministic method. The simulation result is  Figure 14, the field has resource of 45 MWe during 30 years of production based on P10 probability, 72 MWe (P50), and 98 MWe (P90).

Result improvement from the previous model
This numerical model has some significant improvement from the previous model by [10]: • The previous model used 1 km × 1 km size of grid blocks. This new numerical model utilizes finer grid (200 m × 200 m, 200 m × 800 m, and 800 m × 800 m) which produce more detail result. • This model covers smaller area than the previous one to emphasize more on the upflow and surface manifestations area. • As described in Grid Model section above, the surface topography is modelled smoothly in this numerical model while [10] only used ladder gridding in his model. • This numerical model also successfully matches the actual manifestation temperature and gas saturation of fumaroles at the peak of Mount Welirang while the previous model by [10] was unable to do such work. To be noted, the temperature values from [10] are not stated explicitly in the publication. The comparison is shown in Table 4. • Temperature profile result of this numerical model is more reliable than the previous model.
Temperature near the surface in the previous model barely rises as the depth increases up to 1,000 m.

Conclusions
Several conclusions can be drawn from this study: • The new model improves the previous one done by [10].
• Reservoir model shows an upflow to the peak of Mount Welirang and outflow to Padusan, Coban, and Cangar Hot Springs. • The model matches actual manifestation temperatures in Padusan, Coban, and Cangar Hot Springs. It is also able to match the actual fumaroles temperature of which the previous model was unable to do. The model matches reservoir temperature estimate from CO2 gas geothermometry method (around 260 o C). • According to the Monte Carlo simulation, the reservoir resource potential is 45 MWe for 30 years production based on P10 probability, 72 MWe (P50), and 98 MWe (P90).