Numerical simulations of geothermal systems with an unsaturated zone

Abstract

Unsaturated zones are large and important in geothermal systems in mountainous terrain or semi-arid regions. Geothermal features related to the unsaturated zone, such as fumaroles, steaming ground, mud pools, acid springs and acidic hydrothermal alterations, are widespread over all high temperature systems. The unsaturated zone in a geothermal system under Iong term production is enlarged as a result of pressure drawdown in the reservoir. However in general the unsaturated zone has been neglected in all past numerical geothermal simulation studies. The present thesis is a numerical modelling study of the unsaturated zone and related features in geothermal systems. The numerical simulators used in the present study are MULKOM and T0UGH2 (Pruess, 1983, 1991). Several numerical experiments were carried out in order to search for new techniques to improve the computational efficiency of steady state simulations of geothermal reservoirs. It was found that the default parameters of the MULKOM program are at their optimal values. A cascade technique, in which two or more grids of different discretization sizes are used, was found to result in substantial saving in computing time over the conventional one grid approach for steady state simulations. A local phase change procedure which utilizes the local time stepping technique and the box Gauss-Seidel method (BGS) (Brandt, 1982) was found to be very efficient in assisting phase changes. The applicability of multigrid methods (MG) to geothermal reservoir simulations were investigated. MG methods were tested on steady state simulations. A MG implementation, employing two grids, BGS as a relaxation scheme (or a smoother) and the modified FAS scheme, achieved a better convergence rate than the relaxation scheme alone. However, the current MG implementation is still inferior to other methods in terms of computing time. In order to simulate the redistribution at shallow depths of the Upflowing fluid in the presence of an unsaturated zone, the TOUGH2 program was enhanced by including a new thermodynamic module which handles a mixture of water, air and tracer. In order to speed up simulations involving the tracer, a decoupled version was developed for the TOUGH2 program and also for the MULKOM program. The decoupled versions of both programs work much faster than the coupled versions. Two 2D models was set up to represent liquid-dominated geothermal systems with the unsaturated zone of significant extent: those in mountainous terrain (Tongonan type) and those in semi-arid regions (Olkaria type). These two models reproduce the main features of their corresponding conceptual models, such as upflows and outflows, the unsaturated zone and associated surface manifestations (fumaroles and perched springs). Simulation results show strong interactions between the groundwater and geothermal systems. A 2D unsaturated model was set up for the Wairakei geothermal system where the shallow unsaturated zone in the natural state is generally insignificant, but is important in terms of production induced behaviour. The 2D model was able to reproduce the spectacular enhancement of steam heated features (steam flow to the surface) and therefore the associated increase in surface heat loss. The increased steam flow is caused by an increase in the relative permeability of the vapour phase which is in turn caused by the increase in vapour saturation related to pressure drawdown. The increase in steam flow to the surface increases fluid pressure at very shallow depth and therefore the frequency of hydrothermal eruptions. Other factors such as permeability reduction, atmospheric pressure lowering and decrease in infiltration, were investigated to see if they affect production induced shallow hydrothermal eruptions. A coarse 3D unsaturated model was set up for the Tongonan geothermal field. Comparisons were made between the approximate saturated model and a model which includes the unsaturated model for both natural state and production simulations. In the natural state the approximate fully saturated model tends to give much lower temperatures at shallow depth than the equivalent unsaturated model. Under production, the two models predict a similar pressure response in case of a two-phase reservoir. However for a liquid-only reservoir, the fully saturated model underestimates the pressure drop and therefore the unsaturated zone should be included in the model.