The spatial distribution of the geothermal fields in the Taupo Volcanic Zone, New Zealand

Journal of Volcanology and Geothermal Research 145 (2005) 136 – 150 www.elsevier.com/locate/jvolgeores

The spatial distribution of the geothermal fields in the Taupo Volcanic Zone, New Zealand W.M. KisslingT, G.J. Weir Applied Mathematics Group, Industrial Research Limited, PO Box 31310, Lower Hutt, New Zealand Received 15 April 2004; received in revised form 26 August 2004; accepted 17 January 2005

Abstract A numerical model of the large-scale water and heat flows within the Taupo Volcanic Zone (TVZ), New Zealand is constructed. In the model the regions of high permeability and heat flow in the TVZ coincide, and lie within the TVZ envelope, defined to contain the Taupo Fault Belt and the most recent caldera collapse structures. The model correctly predicts the absence of geothermal fields in the central region of the TVZ, where cold surface groundwater descends to depths close to 8 km before being heated and ascending in discrete plumes at the permeability barrier which occurs at the boundary of the TVZ envelope. The locations of these plumes can be identified with most of the major geothermal fields in the TVZ (Tokaanu, Lake Taupo, Wairakei and Tauhara, Rotokawa, Ohaaki, Waiotapu and Waikite, Tikitere, Rotoiti and Rotoma, Rotorua, Mokai). In addition, the model predicts two bands of convective upflow across the TVZ, both of which correspond roughly to the known geothermal features of Ngatamariki and Orakeikorako in the south, and Lakes Rotomahana and Tarawera in the north. D 2005 Elsevier B.V. All rights reserved. Keywords: Taupo Volcanic Zone; geothermal; hydrology; numerical modelling

1. Introduction The Taupo Volcanic Zone (TVZ) is a region of enhanced volcanic and geothermal activity in the North Island of New Zealand (Fig. 1). The TVZ covers an area approximately 30 km wide by 150 km long and has an estimated total heat output of 4200 F 500 MW (Bibby et al., 1995). Volcanic activity

T Corresponding author. Tel.: +64 4 9313351; fax: +64 4 9313003. E-mail address: [email protected] (W.M. Kissling). 0377-0273/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jvolgeores.2005.01.006

is believed to have commenced in the TVZ about 2 Ma ago, but the major rhyolitic volcanism (~15,000–20,000 km3) for which the TVZ is best known began approximately 1.6 Ma ago (Wilson et al., 1995). The geothermal and volcanic activity in the TVZ is associated with the subduction of oceanic crust on the Pacific Plate beneath the Indian–Australian Plate. As the crustal rock descends it melts, and the molten material (magma) rises due to buoyancy forces toward the surface (e.g. Elder, 1976; Stevens, 1980). The TVZ is presently undergoing extension at a rate of approximately 8 mm/year (Darby et al.,

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

137

Fig. 1. Location of geothermal fields (shaded areas), Taupo Fault Belt (grey lines) and known volcanic calderas (dashed black lines) within the TVZ. The geothermal fields are labelled in smaller text, and the caldera names are in capitals. The numerals indicate the natural heat output of each geothermal field, in MW. The inset in the top left corner shows the location of the mapped area in the central North Island of New Zealand.

2000). It is generally accepted that the volume created by this extension is being replaced by the ascending magma at a rate of about 1 m3/s (Stern, 1987). This magma supplies the heat for the geothermal and volcanic systems, although there is some debate (e.g. Hochstein, 1995) whether this volume of magma can in fact supply heat at the required rate.

The TVZ is unique on Earth in that it has a very high heat flow per unit length when compared to other volcanic arc settings. Hochstein (1995) gives a figure of 2600 MW/100 km for the TVZ, and compares this with other areas with well determined heat flows: NE Japan (1400 MW/100 km), Kyushu in southern Japan (800 MW/100 km), Sumatra (~600 MW/100 km) and Oregon (300 MW/100 km).

138

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

Hochstein argues that in most volcanic arc settings, the heat flow due to andesitic and dacitic volcanism is linearly related to the subduction rate, thence determines the danomalousT heat flow in the TVZ associated with rhyolitic volcanism to be 2000 MW/ 100 km. The large volume of rhyolite erupted from the TVZ has averaged 0.28 m3/s in the last 0.34 Ma and this rate is approached only by the Yellowstone area (0.13 m3/s), where the total heat flow is 5300 MW (Wilson, 1996). The role of groundwater convection is recognised as crucial for geothermal heat transport from depth to the surface in the TVZ (e.g. Stewart, 1978). However, so far there have been no detailed quantitative models of this convection. This paper represents the first attempt at a large-scale quantitative description of the TVZ convective system. Previous models of the TVZ hydrological system have been based on linear convection theory (Horton and Rogers, 1945; Lapwood, 1948; Wooding, 1978). These models predict that the geothermal fields are associated with upflows in regular convection cells, which in a homogeneous porous medium are separated by about twice the depth of the heat source. In the TVZ the typical spacing of the geothermal fields is about 15 km, implying a heat source at about 8 km depth. However, these models cannot explain the absence of geothermal activity in the central TVZ, or the west to east distribution of the geothermal heat output. The purpose of this paper is to propose a mechanism which controls the spatial distribution of the TVZ geothermal fields. To do this, a numerical model of the TVZ hydrological system has been created which includes the flows of heat and water, an idealised large-scale permeability distribution and a heat source at a depth of 8 km. To carry out this study, a super-critical equation of state module has been developed for the geothermal simulator TOUGH2 (Pruess, 1991), which can describe the flow of water at the conditions expected deep in the TVZ. Section 2 of this paper outlines the physical setting of the TVZ and the principal geological units of the area. Section 3 describes the heat flows, spatial distribution, and other properties of the TVZ geothermal fields. Section 4 describes the TVZ envelope, which is an essential component of the model proposed here. With this background, section 5 describes the model in detail, and in particular the

heat source and permeability distributions used. Section 6 presents results from the numerical model, and a comparison is made between the modelled geothermal areas and those which actually exist in the TVZ. Finally, we summarise our conclusions, emphasising the essential mechanisms which control the locations of the geothermal fields: the existence of a high-permeability region within the TVZ envelope, the presence of strong downflows of cold meteoric water within this region, and the heat sweep mechanism which gives rise to the geothermal plumes.

2. Physical setting and geology of the TVZ The main features of the TVZ relevant to this paper are shown in Fig. 1. There are 23 geothermal fields, which transport a total of 4200 F 500 MW of heat to the surface by convection of groundwater (Bibby et al., 1995). The locations of these fields are the main focus of this paper, and this and other properties are described in more detail in section 3. The shallow (~100 m) hydrology of the TVZ region is discussed in Rosen and White (2001), but is not relevant to the main theme of this paper, which deals with convection of groundwater to depths of 8 km. The Taupo Fault Belt (grey lines in Fig. 1) is the result of the extensive rifting that has occurred in the TVZ in the last 2 Ma, and is associated with the intrusion of magma into the shallow crust which provides the heat source for the TVZ (Stern, 1987; Wilson et al., 1995). The Taupo Fault Belt is thought to be the primary conduit for the downflows of cool surface waters which ultimately supply the fluid for the hydrothermal convection in the TVZ (Bibby et al., 1995). This conclusion is supported by the observation of very low-temperature gradients with depth within the TVZ compared to more typical conductive gradients outside the TVZ (e.g. Thompson, 1977). Wilson et al. (1995) describe in detail the history of volcanism in the TVZ. The primary form of volcanism in the TVZ is rhyolitic, with more than 15,000 km3 having been erupted in the last 1.6 Ma. Most of this material is thought to be associated with eight calderas—Rotorua, Okataina, Kapenga, Reporoa, Mangakino, Maroa, Whakamaru and Taupo (dark

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

dashed lines in Fig. 1). The oldest of these are Mangakino, in the western TVZ, which was active c.1.6–1.0 Ma ago, and Kapenga, where the earliest activity is dated at 1.0 Ma. Recent activity has been confined to the Okataina and Taupo calderas, and all others were active in the dyoung TVZT period between 0.34 and 0.1 Ma ago (Wilson et al., 1995). In the context of this paper, these volcanic centres provide clues to the nature of the permeability within the TVZ and its underlying heat source, as will be described in Section 4. In the TVZ heat and chemical constituents from the deep magma are transported to the surface by convection of the groundwater in the upper few kilometres of the Earth’s crust. Isotopic evidence (Stewart, 1978) shows that most of this water is ultimately derived from rainfall. The depth to which the groundwater circulates is determined by the transition from brittle to ductile behaviour in the crustal rocks, as ductile material is unable to support the fractures and cracks which produce permeability. Bibby et al. (1995) show the depth distributions for a number of earthquake swarms with reliably determined depths, including the Edgecumbe aftershocks and the 1983 Taupo swarm. A sharp drop off in the number of swarm events occurs below 8–9 km depth in most of these distributions, including that at Taupo. This is interpreted as being due to the transition between brittle and ductile behaviour in the rocks at this depth. This is supported by inferences made by Villamor and Berryman (2001), who conclude, based on fault geometry and extension rates, that extension of major faults in the TVZ occurs to depths of 6–10 km. It is widely held that ductile rock cannot sustain long– term permeability (Fournier, 1991), and so this depth is considered to be the lower limit for hydrothermal fluid circulation within the TVZ. Convection of groundwater above this level gives rise to the TVZ geothermal fields which are the main focus of this paper. The idea that the heat source and permeable regions are the same size is consistent with models which relate the spreading of the TVZ to the heat generation process (e.g. the TVZ model of Stern, 1987). In these models the magma which provides the heat source for the TVZ is intruded into the same volume where permeability is being created by tectonic extension. Recently Bibby et al. (2000) have used magnetotelluric data to investigate the deep

139

structure of the TVZ. This study shows a dsignificant conductive structureT at 10 km depth below the Taupo Fault Belt, which is interpreted as possibly being due to magma in the crust. Another conductive body below 20 km depth lies under the eastern TVZ, where present hydrothermal activity is highest, and is interpreted as a zone of dconnected meltT in the upper mantle. These structures are too deep to be of direct importance to the modelling in this paper, which extends only to a depth of 8 km. 2.1. Infill region Long-wavelength (tens of kilometres) gravity anomalies in the TVZ are associated with a broad region of low-density volcanic infill, overlying denser basement material. The magnitude of this anomaly is about 350 AN/kg, and leads to estimates of the depth of the infill region of 1–2 km (Modriniak and Studt, 1959), or 2.5 km (Bibby et al., 1995). This picture is consistent with the limited drilling data available. Basement rock (greywacke) has been penetrated to depths of about 2 km in the eastern TVZ at Rotokawa, Ohaaki, Kawerau and Ngatamariki, but is not found in wells deeper than 2 km in the western TVZ at Mokai and elsewhere (Wood, 1996), suggesting, within the limitations of the data, that the shallow infill thickens toward the west. 2.2. Exterior region The region outside the TVZ infill region comprises greywacke ranges. To the east, the Kaingaroa Plateau extends to approximately 1000 masl and to the west, greywacke outcrops west of the Mangakino caldera.

3. The TVZ geothermal fields The spatial distribution of the TVZ geothermal fields is the main theme of this paper. In this section the properties of the geothermal fields are described, with particular emphasis on those aspects which later numerical modelling in this paper may be able to predict. A good summary of the principal characteristics of the TVZ geothermal fields can be found in Bibby et al. (1995).

140

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150 1200

Bibby et al., (1995) list 23 geothermal fields within the TVZ. Heat flows from individual fields range from b1 MW (Motuoapa) to 540 MW (Waiotapu). Ten of the smaller fields (b100 MW) constitute less than 10% of the TVZ heat flow, whereas over 40% of the thermal output of the TVZ comes from just four geothermal fields with heat flows greater than 400 MW, Wairakei, Waiotapu, Mokai and Rotorua. A histogram of the heat flows from the TVZ geothermal fields is shown in Fig. 2. Fig. 1 gives the most recent heat flow figures for each geothermal field, as listed in Bibby et al. (1995). The geothermal fields are distributed in two bands, roughly 20 km apart, and parallel to the eastern margin of the TVZ, as shown in Fig. 3. Currently 75% of the heat output occurs in the eastern band and 25% in the west (Bibby et al., 1995). The only known extinct geothermal fields, Ohakuri and West Rotorua, and six low-power (b20 MW) fields, all lie in the western TVZ. Careful examination of Fig. 4 shows that many of the geothermal fields lie on the margin of the Taupo Fault Belt. High electrical resistivity in the regions between the geothermal fields indicates that these areas have never undergone hydrothermal alteration (Bibby et al., 1995). The implication is that the geothermal systems have been confined to the present regions of activity

1000

10

Number of geothermal systems

9 8 7 6 5 4 3 2 1 0

0

100

200

300

400

500

600

Heat Flux (MW) Fig. 2. Histogram showing the number of geothermal fields in the TVZ, when grouped into 100 MW bins.

Heat Output (MW)

3.1. Heat flows and spatial distribution

800

600

400

200

0 -40

-30

-20

-10

0

Distance from Eastern Margin (km)

Fig. 3. The spatial distribution of heat flows (MW) from geothermal fields in the TVZ, from Bibby et al. (1995). The heat flow is bimodal, with approximately 75% of the total flow occurring in the eastern TVZ.

during the lifetime of the TVZ. It is thought that hydrothermal plumes, once formed, are stable, long lived features. 3.2. Temperatures The highest temperatures encountered in geothermal wells in the TVZ are over 300 8C—for example at Mokai (326 8C), Kawerau (315 8C) and Rotokawa (330 8C) (e.g. Ellis and Mahon, 1977). At Wairakei, the maximum temperature is about 265 8C (Grindley, 1965). The highest temperatures tend to occur in the eastern geothermal fields (typically 310 8C), while 265 8C is more common in the western geothermal fields. Low-temperature gradients (~20 8C/km) are observed in the deeper geothermal wells, for example at Wairakei (e.g. McNabb, 1992). Extrapolating these temperature gradients to depth suggests temperatures of the order of 350–400 8C at 8 km depth, which is consistent with the brittle–ductile transition occurring at that depth. Thus, large-scale magma-like temperatures (~800–1200 8C) are not expected within the region where convection of groundwater is the dominant form of heat transport. However, it should be remembered that no temperature data has been obtained from deeper than about 2.5 km.

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

141

Fig. 4. Relationship of the geothermal fields, Taupo Fault Belt and known volcanic calderas to the TVZ envelope (shown as heavy line surrounding these features). In the model presented in this paper, the features inside this envelope are associated with enhanced permeability. The computational domain of the model (inclined rectangle) is aligned with the Taupo Fault Belt, and is inclined approximately 408 to the border of the map, and extends partially beyond the boundaries of this figure. In later figures in this paper, the computational domain is presented vertically, that is independently of the geographical coordinate system.

3.3. Areas of the geothermal fields Shallow electrical resistivity studies show typical areas of 5–25 km2 for the TVZ geothermal fields (Bibby et al., 1995). These areas are measured at depths between 0.5 and 1 km, and result from the

interaction of the rising geothermal plumes with surface groundwater and shallow geological capping structures. A good example of this is at Mokai, where entrainment of geothermal fluid into the shallow groundwater system results in a resistivity boundary at least 10 km2 in extent (Bibby et al., 1981).

142

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

These areas do not reflect the areas of the rising plumes at greater depths. Data from deep wells at a number of geothermal fields (e.g. at Wairakei; Grindley, 1965) show temperature inversions below a few hundred metres depth, and the highest temperature wells (and therefore presumably the upflow regions) are often confined to regions of the order of 1 km across. 3.4. Ages of the geothermal fields There is evidence that some geothermal fields have lifetimes of several hundred thousand years. Grindley (1965) gives 0.5 Ma as the minimum age of Wairakei, based on fragments of hydrothermally altered ignimbrite found in younger sediments. Browne (1979) suggests that hydrothermal activity has occurred at Kawerau for at least 0.2 Ma. This estimate is based on the presence of greywacke fragments containing veins of hydrothermal minerals in strata at least this old. At Rotorua, the presence of altered rock at the old, highlevel, lake shore shows that thermal activity has existed in the area for at least the last 22–50 ka (Wood, 1992). Because there have been changes in volcanic activity in the TVZ in the last 2 Ma it is not known if the present thermal outputs of the geothermal fields are typical of those over a longer period. Wilson et al. (1995) suggest on this basis that geothermal activity is probably also quite variable. However, the geothermal fields do exhibit some stability. For example, the geothermal field at Wairakei appears to have been active through the period of the Taupo eruption in AD 186, and more recently, the Waimangu geothermal field did survive the 1886 eruption of Mount Tarawera (e.g. Keam, 1988), although with large changes in surface activity.

4. The TVZ envelope Fig. 1 shows the locations of the geothermal fields in relation to the Taupo Fault Belt and the known caldera structures in the TVZ. There is some association between caldera and fault zone boundaries and the locations of the geothermal fields. All of the dmajor Tgeothermal fields in the eastern TVZ, from Wairakei/Tauhara to Waimangu are associated with the eastern margin of the Taupo Fault Belt. Similarly, the northern fields, Rotorua, Tikitere and Rotoma lie

close to the boundaries of the Okataina, Rotorua and Kapenga volcanic centres, as do the Taupo and Tokaanu fields to the Taupo caldera boundary. Mokai is located on the western margin of the Taupo Fault Belt (TFB). Of the small geothermal fields, Ongarato, Atiamuri and Horohoro are also associated with the western boundary of the TFB. The region defined by the TFB and the volcanic calderas will henceforth be referred to as the dTVZ envelopeT, and will play a central role in the modelling in this paper. The TVZ envelope is shown in Fig. 4 as the dark blue line. There are three geothermal fields which do not fit with this general picture, and the positions of these fields may be governed by entirely different processes. The Kawerau geothermal field lies in the north-east of the TVZ, and is apparently outside the presumed boundary for the Okataina volcanic centre. One explanation for this is that the location of the field results from a shallow magmatic intrusion, as suggested by Christenson et al. (1998). A more mundane explanation is that there is unmapped volcanic or faulting structure extending north-east from the Okataina volcanic centre. Reporoa, in the eastern TVZ, is a small geothermal field which lies inside a caldera structure, but outside the TVZ envelope as it is presently defined. Reporoa is unusual in two respects. Firstly, it has an unusually low total heat flow (15 MW) in relation to the natural geophysical heat flow from a 10 km diameter caldera, which is approximately 7 MW. Secondly, there is an apparent connection with the Waiotapu–Waikite system, to which Reporoa is clearly linked by electrical resistivity surveys (Stagpoole and Bibby, 1998). Reporoa may represent an doutflowT from the Waiotapu system (Healy and Hochstein, 1973). Similarly, Mangakino geothermal field lies within, rather than on, its caldera boundary, and is outside the TVZ envelope as it is defined above. This geothermal field is very small and old (2 Ma), being associated with the earliest volcanism in the TVZ (Wilson et al., 1995). Because of its small heat output and western location, Mangakino will have little influence on the overall heat and mass flows in the TVZ. In addition, the natural geophysical heat flow from the Mangakino caldera (about 36 MW) is much larger than the geothermal output of 2 MW. There is no mapped geothermal heat flow within the Tihoi caldera. According to Soengkono (1995),

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

this has a similar magnetic signature to Mangakino, and was dprobably active at the same time, 1.62–0.87 Ma ago.T Tihoi does not appear on the maps of Rogan (1982) or Stern (1987). Because of its uncertain nature and the absence of geothermal heat flow, Tihoi caldera has not been used as part of the definition of the TVZ envelope. The above definition of the TVZ envelope depends on the positions of the faults within the TFB and the presently known caldera boundaries given in Fig. 4. The caldera boundaries are taken from Soengkono (1995) and the locations of the faults within the Taupo Fault Belt and the geothermal fields are from Bibby et al. (1995). These represent the latest available information on these structures. It is possible that other unmapped calderas exist in the central TVZ, but these are likely to be masked by more recent volcanism and therefore will be difficult to find (Wilson et al., 1995). Should such structures be discovered, the model described in this paper would need to be modified accordingly. The region within the TVZ envelope is assumed to be of high permeability relative to the that outside the envelope because fault movement and caldera collapse will generate permeability. The depth of this enhanced permeability within the Taupo Fault Belt may extend down to 8 km, but the depth of enhanced permeability in the caldera regions is likely to be somewhat less.

Kissling (2004) showed that discrete geothermal plumes occur at caldera boundaries when the heat source is coincident with the high-permeability region which defines the caldera. In this paper the same assumption is made, and the TVZ geothermal heat source is defined to exist only in the region defined by the TVZ envelope. The model has a spatially variable heat source with a total heat flow of 4289 MW. The heat is applied piecewise in separate regions across the lower boundary of the model. To do this, the model domain was first divided into ten 15 km wide strips in the SWNE direction along the TVZ, and within each of these strips the TVZ envelope was further divided into two approximately equal areas, as shown in Fig. 5.

8

16

7

15

6

5

14

12

3

11

5.1. Model heat source The TVZ heat source is highly variable in time and space, as evidenced by the sporadic nature of the TVZ volcanism (Wilson et al., 1995). Weir (1998) has modelled a spatially averaged heat source which links the rate of magmatic intrusions to the TVZ extension rate. However, a complete treatment of the variability of the heat source is not attempted in this paper. Instead, a simpler approach is taken, where contributions to the total heat flow into the TVZ from magmatic fluid, magmatic intrusions and conduction are lumped together and treated as a steady effective heat source at the base of the model.

13

13

4 5. The TVZ model

143

2

10

1

9

Fig. 5. Assumed heat source and geological regions within the TVZ envelope (see Fig. 4). The envelope is divided into 16 regions, as indicated. The geothermal heat source within each region is placed at 8 km depth, and its strength is proportional to the measured heat flux in that region. The permeable dinfillT region of the model is also defined by this envelope. The heat flows and permeabilities for each region are listed in Table 1. This figure must be rotated clockwise by 408 to match the TVZ maps in Figs. 1 and 4.

144

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

Heat flows were then assigned to each of these 16 regions. These flows were calculated by summing the observed geothermal heat flows (taken from Bibby et al., 1995) within each region. The division of the heat source in this manner is used to investigate possible reasons for the west to east distribution of heat flow in the TVZ. A summary of this calculation is given in Table 1. For all geothermal heat sources, there is a corresponding inflow of magmatic water. Using typical figures of 6% and 14% for the fractions of magmatic water in the western and eastern geothermal fluids (Giggenbach, 1995), a total flow of magmatic fluid of about 400 kg/s for the whole TVZ can be derived. This is roughly equivalent to 0.1 kg/s per MW of TVZ heat flow, and magmatic water is applied at this rate for all the geothermal heat sources in the model. 5.2. Model permeability Inside the TVZ envelope the permeability is defined piece-wise in the same 16 regions as the heat source. Details are given in Table 1 and Fig. 5. Permeabilities are set to be proportional to the heat flux for each region—the higher the heat flux, the higher the permeability. There are several reasons for Table 1 Summary of areas (km2), total heat flows (MW), heat fluxes (W/m2) and permeabilities adopted for the 16 heat source regions of the TVZ model (see Fig. 5) Region

Area

Heat

Flux

Perm.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

36 158 241 280 152 152 241 116 152 208 237 240 164 135 239 174

0 0 400 2 10 4 420 423 200 200 530 680 370 625 325 0

0 0 1.66 0.01 0.07 0.03 1.74 3.65 1.32 0.96 2.24 2.83 2.26 4.63 1.36 0

5/0.5 5/0.5 25/2.5 5/0.5 5/0.5 5/0.5 25/2.5 75/7.5 25/2.5 5/0.5 50/5 50/5 50/5 100/10 25/2.5 5/0.5

The permeabilities are given in the form horizontal/vertical in units of 10 15 m2. All regions outside the heat source region are assigned a permeability of 1/0.1  10 15 m2.

this proportionality. Firstly, the numerical model becomes impractically slow if lower permeabilities are used. This is because conductive regions with very low-permeability form when the temperature exceeds 350 8C, due to the temperature dependence of permeability imposed in the model (Fournier, 1991). Forcing the fluid through these regions is difficult both physically and numerically. There are other arguments against the presence of these conductive regions—for example, it is difficult to explain the observed magmatic component of the fluid chemistry of the TVZ geothermal fields (e.g. Giggenbach, 1995) if the fluids are merely heated by contact with these conductive regions. Further, magneto-telluric observations (Bibby et al., 2000) show no evidence for their presence. Lastly, the TVZ is an actively rifting tectonic region and implicit in our model is the assumption that the heat source and rifting region are coincident. Modelling of TVZ-like systems in 2D and 3D (Kissling, 2004) shows this has to be the case to get the discrete geothermal plumes to form properly. Because of this, it is reasonable to suggest that permeability is created during this rifting process (the TVZ is a very faulted region), and that it is likely to be higher in regions of greater heat flow. The direct proportionality between heat flow and permeability was the simplest assumption consistent with this. The relationship between the heat flux and the assigned permeability is shown in Table 1. For the region outside the TVZ envelope, the horizontal/vertical permeability has been set to 1/ 0.1  10 15 m2. This exterior permeability provides an adequate permeability contrast between the interior and exterior regions for stable plumes to occur, and also ensures that the temperature in the exterior region remains close to its initial conductive distribution, thus avoiding the convective flows seen in this region when higher permeabilities are used (Kissling, 2004). As already discussed, the Mangakino, Tihoi and Reporoa calderas lie outside the TVZ envelope, and are assigned low permeabilities equal to the rest of the exterior region on the basis of their age, low heat flow and current activity. 5.3. Computational details The model presented in this paper represents the complete onshore section of the TVZ and the

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

150 km

surrounding region. The model domain covers an area of 150 km  80 km and is 8 km deep, consistent with the inferred depth of the brittle–ductile transition. The element size throughout the model is 1 km  1 km  0.5 km, and these are arranged in 16 layers. A plan view of the element layout for a single layer of elements is shown in Fig. 6. Excluding boundary elements, the total number of elements is 192,000. This model domain is aligned SW–NE along the TVZ, and is shown by the inclined rectangle in Fig. 4. The upper boundary of the model is held at datmosphericT conditions—fully saturated at a pressure of 1 bar and a temperature of 20 8C. A fully saturated upper boundary is used because the average annual rainfall is approximately fifty times the infiltration rate of water into the TVZ required to supply the fluid for the geothermal systems (Kissling, 2004). This is consistent with the observation of Stewart (1978), that most geothermal waters have a meteoric origin. Simulations of 2D slices across the TVZ (Kissling, 2004) show that topographically driven flows have little effect on the overall TVZ

145

convective system, which is dominated by the geothermal heat source, and for this reason the upper boundary of the model does not include any topographical features. The lateral boundaries are held at hydrostatic pressure corresponding to the geophysical temperature gradient of 30 8C per km. These conditions also serve as the initial conditions in the interior of the model domain. On the lower boundary, there are two heat sources. The first has a strength of 0.09 W/m2 and is distributed uniformly over the lower boundary. This is the conductive heat flow associated with the geophysical temperature gradient of 30 8C/km which is observed in regions adjacent to, but outside, the TVZ (Thompson, 1977). The second dgeothermalT heat source has been described in section 5.1. To carry out the simulations described in this paper, a super-critical equation of state module has been developed for the geothermal simulator TOUGH2 (Pruess, 1991), which can describe the transport of water at the conditions expected deep in the TVZ. TOUGH2 is a fully implicit integrated finite difference code which solves the equations of mass and energy conservation in a porous medium, supplemented with appropriate constitutive relationships for the density, enthalpy and dynamic viscosity of water. The transport of water is described by Darcy’s law. Except near the critical point of water, the thermodynamic description of water in this code is based on the United Kingdom Steam Tables in SI Units, 1970) formulation for pressures less than 1000 bar and temperatures less than 800 8C. For conditions outside this range, the NBS/NRC formulation of Haar et al. (1984) is used. In the immediate vicinity of the critical point ( P = 221.2 F 1 bar, T = 374.15 F 0.005 8C) an approximation to the density of water is used because of numerical difficulties encountered with the UK Steam Tables formulation. Kissling (2004) contains details of this approximation and other features of the code, including some test examples in the nearcritical region, and validation of the code against results published in the literature.

80 km Fig. 6. Computational grid for the TVZ model. Each layer of the model consists of 12,000 1 km  1km  0.5 km elements, and there are 16 such layers. The TVZ envelope is also shown. This figure must be rotated clockwise by 408 to match the TVZ map in Fig. 4.

6. Modelling results The model presented in this paper was run to simulate a period of 2 Ma, approximately the time

146

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

since volcanism started in the TVZ (Wilson et al., 1995). Fig. 7 shows the temperature at a depth of 0.5 km in the model at 2 Ma. In this case there are a number of stable high-temperature plumes around the perimeter of the TVZ envelope. These have a maximum temperature of about 220 8C. The cold downflows again appear in the eastern TVZ, in regions where the permeability is highest. Fig. 8 shows the temperatures at a depth of 7.5 km in the model at 2 Ma. The regions of strongest downflow (dark blue) are apparent. The higher temperature regions inside the TVZ envelope and about its margin correspond mainly to upflows of hot fluid, and correlate well with the positions of the hot plumes at shallow depths in Fig. 7. The exterior region of the model has an approximately uniform temperature distribution, with a mean temperature of about 230 8C. This is consistent with the undisturbed initial temperature at this depth, and shows that the flows of heat and mass to the exterior of the TVZ envelope are rather small. In this model most of

Fig. 7. Plan view of the temperature distribution at a depth of 0.5 km at 2 Ma for the TVZ model. The heat flows for each of the numbered plumes are listed in Table 2. This figure must be rotated clockwise by 408 to match the TVZ map in Fig. 4.

Fig. 8. Plan view of the temperature distribution at a depth of 7.5 km at 2 Ma for the TVZ model. This figure must be rotated clockwise by 408 to match the TVZ map in Fig. 4.

the fluid which comprises the TVZ geothermal systems originates from the surface within the highpermeability TVZ envelope. From the surface it descends to a depth close to 8 km, is heated, and rises in discrete plumes around the envelope boundary to form the geothermal systems. The temperature at 7.5 km depth within the TVZ envelope is mostly below 200 8C, in contrast to the expected temperature of ~350 8C. This is largely because of the assumption that the permeability in the envelope is constant with depth. At present, the model is so computationally demanding that it has not been possible to experiment with other permeability distributions. The assumption that the permeability is constant with depth is the simplest, and in the absence of other information serves the main purpose of this paper, which is to identify the mechanism which controls the locations of the geothermal fields. Higher temperatures would be expected if the permeability decreased (say exponentially) with depth, but this calculation has not yet been carried out. Fig. 7 shows a good correspondence to the TVZ. The geothermal fields occur on the internal boundary of the TVZ envelope, and many of these can be

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150 Table 2 Comparison of predicted and observed heat flows (MW) in the TVZ Plume

Identification

Actual heat

Model heat

1 2 3 4 5 6 7 8 9

Tokaanu Lake Taupo Wairakei/Tauhara Rotokawa Ohaaki Waiotapu/Waimangu/Waikite Tikitere/Rotoiti/Rotoma Rotorua Mokai

200 200 530 300 70 935 410 420 420

211 244 191 632 446 1125 406 725 385

147

system (3). The next plumes correspond to Rotokawa (4) and Ohaaki fields (5), with Ohaaki again appearing at the apex of two of the faces of the TVZ envelope. Moving northward a series of four or five small plumes correspond to the geothermal fields at Waiotapu, Waikite and Waimangu (6). The northern most modelled plume lies close to the fields at Tikitere, Rotoiti and Rotoma (7). Next, moving south, the double plume contains the heat of the Rotorua geothermal system (8). The last significant hot plume occurs close to the location of the Mokai geothermal field (9). A 3D rendering of Fig. 7 is shown in Fig. 9. This shows temperature isosurfaces for T = 220 8C, which delineate the structure and location of the geothermal plumes more clearly. The geothermal fields are labelled. In addition the temperature on a horizontal slice through the model at a depth of 3 km is shown. This is approximately the greatest depth to which wells have been drilled in the TVZ, and shows clearly the cold downflows of surface water in the central/ western TVZ area.

The plume numbers are given on Fig. 7.

identified with those in Fig. 4. The plumes are listed by number, and these are referred to in Table 2. Comparing these figures and starting at the southern most point of the TVZ envelope, the first geothermal field appears at the apex of the envelope. This corresponds to the Tokaanu geothermal field (1). Moving anti-clockwise, the next large plume corresponds to the Lake Taupo field (2). Two smaller plumes are at the location of the Wairakei/Tauhara

Tikitere/Rotoiti/Rotoma Rotorua Mokai

-2

Z (km)

-4 -6

150

-8 80

100 60

Y(

km40 )

50

Temp (C) 240 220 200 180 160 140 120 100 80 60 40 20

km)

X(

20 0

0

Waiotapu/Waimangu/Waikite

Wairakei/Tauhara/Rotokawa Tokaanu LakeTaupo

Fig. 9. 3D view of the TVZ model after 2 Ma of simulation time. The structure of the geothermal plumes are shown clearly by the T = 220 8C isosurface. The slice shows the temperature at 3 km depth, with the dark blue area showing the downflow of cold surface water in the western/ central TVZ. The geothermal plumes are indicated. Refer to Table 2 and Fig. 7 for more details. The vertical scale has been exaggerated for clarity.

148

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

The heat flows for each of these geothermal plumes can be calculated. These are obtained by summing the heat flow to the surface for each plume, for all elements with a heat flow of more than 5 MW. These are listed in Table 2, together with the didentificationT of the plume with its actual TVZ counterpart(s). The Table shows that the heat flow from many of the modelled geothermal fields (Tokaanu, Lake Taupo, Mokai), or groups of fields (Waiotapu, Waimangu, Waikite and Tikitere, Rotoiti, Rotoma), are in good agreement with measurements. The calculated heat flow from the Wairakei/Tauhara system is too low, and the heat flows from Rotokawa, Ohaaki and Rotorua are too large. Except for Ohaaki, the heat flows agree with observation to within a factor of three, and the agreement is usually much better. The model does not reproduce the geothermal fields with low heat flow in the western TVZ such as Atiamuri or Horohoro. Thus, this model is quite successful in predicting the location of the large geothermal systems, and the calculated heat flows from these areas are in reasonable agreement with observation. Inside the TVZ envelope, the model predicts a nearly continuous band of hot fluid at the surface which crosses the TVZ between the modelled geothermal fields dRotokawaT and dMokaiT. In the actual TVZ this fluid appears in discrete plumes as the Ngatamariki, Orakeikorako and (possibly) Te Kopia geothermal fields. These fields are slightly to the north of the dbandT, and have a total measured heat flow of 530 MW. The modelled heat flow is only 184 MW, although this figure depends on exactly how each of the Mokai and Rotokawa plumes are defined. In the northern part of the TVZ envelope, a similar but rather weaker band of hot fluid appears in the interior between the modelled dWaimanguT and dRotoruaT systems. This heat flow might be associated with the known heat flows in Lakes Rotomahana and Tarawera (e.g. Whiteford et al., 1996). The total upward flow of water associated with the geothermal plumes is calculated to be 4750 kg/s, and the total downflow of surface water inside the TVZ envelope is 3870 kg/s. The balance comprises magmatic water supplied with the heat source (430 kg/s) and lateral flows from the exterior region into the TVZ envelope (450 kg/s). The induced geothermal flow of 4750 kg/s implies an average magmatic fluid

content of 9%. This is consistent with Giggenbach’s (1995) values of between 6% and 14% for the western and eastern geothermal fluids respectively, and its implication that the majority of water from which the geothermal plumes form in the TVZ circulates to the depth of the magmatic source ~8 km. In the model, the temperature in the surface plumes is in the range 200–230 8C. This is compared with the 265–310 8C found in the TVZ geothermal fields. Reducing the permeability of the exterior region outside the TVZ envelope might reduce the water available to form the geothermal plumes and raise their mean enthalpy. Another possible reason for the low temperatures is insufficient spatial resolution in the model. Experiments with simpler TVZ models in 2D settings (Kissling, 2004) show that the temperature of the plumes increases as the mesh size is decreased, but such a reduction is not yet possible with models of the size presented in this paper.

7. Conclusions In this paper a model of the complete TVZ hydrological system has been constructed. This incorporates the idea of a dTVZ envelopeT, which is defined as the region including the most recent caldera structures within the TVZ and the Taupo Fault Belt. The region within this envelope has moderate permeability and the exterior region has a lower permeability. A dheat sweepT mechanism is proposed that leads to the formation of stable hydrothermal plumes. In this model, cool downflows from the surface in the interior of the TVZ are heated at depth, and the hot fluid is then swept outward to the permeability contrast where it ascends, rather than moves laterally, because of the low-permeability barrier. The model includes spatial variations in both the TVZ heat source and the permeability within the TVZ envelope as a means of explaining the highly variable heat output across the TVZ. The heat flows are varied across the TVZ in accordance with the observed geothermal heat flow, and the permeabilities within the envelope are proportional to the heat flows. This prevents the formation of low-permeability conductive regions above areas of high heat flux which geochemical and geophysical evidence suggest do not occur in the TVZ. Because of the contrast between the

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

low-permeability region exterior to the TVZ, and the higher permeability inside the TVZ envelope, this heat source distribution leads to downflows of cool surface waters in the central TVZ, and the formation of discrete geothermal plumes about the envelope boundary. Both of these are consistent with observed features of the TVZ. The model predicts that geothermal systems lie around the boundary of the TVZ envelope, and exist as stable, discrete plumes a few kilometres across, persisting for periods of the same order as the age of the TVZ. The positions of the larger groups of geothermal fields are well predicted. The heat flows from the major geothermal fields are usually predicted to within a factor of two to three (Table 2), but the small geothermal fields (~10 MW) such as Horohoro and Atiamuri are entirely absent. The model predicts large downflows of surface waters in the western and central TVZ which supply most of the water for the geothermal systems. This is in accord with the observation (Bibby et al., 1995) that there is a dnegative correlationT between regions of high heat flow and major faulting in the TVZ. The temperatures of the geothermal plumes are about 220 8C, compared with the 260–320 8C observed in the TVZ. These temperatures might be obtained by reducing either the permeability exterior to the TVZ envelope, or perhaps increasing the spatial resolution of the model. In the TVZ the deep hot upflow zones of the geothermal fields are typically only 1 km across, and so model elements smaller than this are required to fully resolve these features. At present the significant increase in the number of model elements required to do this is not practical, as run times with the present model are typically 4–6 weeks on the University of Auckland’s highperformance computer (IBM Regatta). In the model there is a weakly defined large-scale system of three convection cells superimposed on the TVZ. Downflows of cool surface fluid occur in the middle of these cells, and the cell boundaries are formed by linear dbandsT of hot fluid. These bands appear to be controlled by the geometry of the TVZ envelope, as they occur also in models with uniform permeability and/or heat source distribution. For the TVZ model described in this paper, these bands correspond approximately to actual geothermal areas in the TVZ (Ngatamariki and Orakeikorako in the

149

southern band, and the thermal areas in Lakes Tarawera and Rotomahana in the northern band), although the predicted heat output from these fields is not very accurate. Some geothermal areas, such as that at Kawerau, do not fit with the model presented in this paper. This may be because of unmapped faults within the TVZ extending north from the Okataina volcanic centre, or possibly because Kawerau is controlled by a heat source related to the nearby active volcano Putauaki (Mt. Edgecumbe). Mangakino and Reporoa geothermal fields are also dproblemsT, although they are small fields which cannot influence the overall TVZ hydrology significantly.

Acknowledgments This work was initially funded by the New Zealand Foundation for Research, Science and Technology. The authors also wish to thank the New Zealand Institute of Mathematics and its Applications for facilitating access to the University of Auckland’s high-performance computer.

References Bibby, H.M., Dawson, G.B., Rayner, H.H., Stagpoole, V.M., Graham, D.J., 1981. Geophysical investigations of the Mokai geothermal field. DSIR, Geophysics Division, Report 184. with appendix by Stagpoole,V.M. and Grantham, A.R, 41 pp. Bibby, H.M., Caldwell, T.G., Davey, F.J., Webb, T.H., 1995. Geophysical evidence on the structure of the Taupo volcanic zone and its hydrothermal circulation. J. Volcanol. Geotherm. Res. 68, 29 – 58. Bibby, H.M., Risk, G.F., Caldwell, T.G., Ogawa, Y., Takakura, S., Uchida, T., 2000. Deep electrical structure beneath the Taupo volcanic zone and the source of geothermal heat. Proc. 19th New Zealand Geothermal Workshop, pp. 81 – 86. Browne, P.R.L., 1979. Minimum age of the Kawerau Geothermal Field, North Island, New Zealand. J. Volcanol. Geotherm. Res. 6, 213 – 215. Christenson, B.W., Wood, C.P., Arehart, G.B., 1998. Shallow magmatic degassing: processes and PTX constraints for paleofluids associated with the Ngatamariki Diorite Intrusion, New Zealand. In: Water–Rock Interaction 9, Taupo. A.A. Balkema, Rotterdam, pp. 435 – 438. Darby, D.J., Hodgkinson, K.M., Blick, G.H., 2000. Geodetic measurement of deformation in the Taupo volcanic zone: the north Taupo network revisited. N.Z. J. Geol. Geophys. 43, 157 – 170.

150

W.M. Kissling, G.J. Weir / Journal of Volcanology and Geothermal Research 145 (2005) 136–150

Elder, J., 1976. The Bowels of the Earth. Oxford University Press, London. 222 pp. Ellis, A.J., Mahon, W.A.J., 1977. Chemistry and Geothermal Systems. Academic Press, New York. 392 pp. Fournier, R.O., 1991. The transition from hydrostatic to greater than hydrostatic fluid pressure in presently active continental hydrothermal systems in crystalline rock. Geophys. Res. Lett. 18, 955 – 958. Giggenbach, W.F., 1995. Variations in the chemical and isotopic composition of the fluids discharged from the Taupo volcanic zone, New Zealand. J. Volcanol. Geotherm. Res. 68, 89 – 116. Grindley, G.W., 1965. The geology, structure and exploitation of the Wairakei geothermal field, Taupo, New Zealand. New Zealand Geological Survey. Bulletin8 vol. 75. DSIR. Haar, L., Gallagher, J.S., Kell, G.S., 1984. NBS/NRC Steam Tables. Hemisphere Publishing Corp. 320 pp. Healy, J., Hochstein, M.P., 1973. Horizontal flow in hydrothermal systems. J. Hydrol., N.Z. 12 (2), 71 – 82. Hochstein, M.P., 1995. Crustal heat transfer in the Taupo volcanic zone (New Zealand): comparison with other volcanic arcs and explanatory heat source models. J. Volcanol. Geotherm. Res. 68, 117 – 151. Horton, C.W., Rogers, F.T., 1945. Convection currents. A porous medium. J. Appl. Phys. 16, 367 – 370. Keam, R.F., 1988. Tarawera: The Volcanic Eruption of 10 June 1886. R.F. Keam, Auckland. 472 pp. Kissling, W.M., 2004. Deep hydrology of the geothermal systems in the Taupo Volcanic Zone, New Zealand. Unpublished PhD thesis. University of Auckland. Lapwood, E.R., 1948. Convection of a fluid in a porous medium. Proc. Camb. Philol. Soc. 44, 508 – 521. McNabb, A., 1992. The Taupo–Rotorua hotplate. Proc. 14th New Zealand Geothermal Workshop, pp. 111 – 114. Modriniak, N., Studt, F.E., 1959. Geological structure and volcanism in the Taupo–Tarawera district. N.Z. J. Geol. Geophys. 2, 654 – 684. Pruess, K., 1991. TOUGH2—a general-purpose numerical simulator for multi-phase fluid and heat flow. Lawrence Berkeley Laboratory. Report LBL-29400. Rogan, M., 1982. A geophysical study of the Taupo volcanic zone, New Zealand. J. Geophys. Res. 87 (B5), 4073 – 4088. Rosen, M.R., White, P.A. (Eds.), Groundwaters of New Zealand. New Zealand Hydrological Society Inc., Wellington. 498 pp.

Soengkono, S., 1995. A magnetic model for deep plutonic bodies beneath the central Taupo volcanic zone, North Island, New Zealand. J. Volcanol. Geotherm. Res. 68, 193 – 207. Stagpoole, V., Bibby, H.M., 1998. Shallow resistivity of the Taupo volcanic zone, New Zealand. Proc. 20th New Zealand Geothermal Workshop, pp. 303 – 309. Stern, A.T., 1987. Asymmetric back-arc spreading, heat flux and structure associated with the central volcanic region of New Zealand. Earth Planet. Lett. 85, 265 – 276. Stevens, G., 1980. New Zealand Adrift. A.H. and A.W. Reed Ltd, Wellington. 442 pp. Stewart, M.K., 1978. Stable isotopes in waters from the Wairakei geothermal area, New Zealand. Stable Isotopes in Earth Sciences, Bulletin8 vol. 220. DSIR, pp. 113 – 119. Thompson, G.E.K., 1977. Temperature gradients within and adjacent to the north island volcanic belt. N. Z. J. Geol. Geophys. 20, 85 – 97. U.K. steam tables in SI units, 1970. Edward Arnold Ltd, London. 161 pp. Villamor, P., Berryman, K.R., 2001. A late Quaternary extension rate in the Taupo volcanic zone, New Zealand, derived from fault-slip data. N. Z. J. Geol. Geophys. 44, 243 – 269. Weir, G.J., 1998. Energy transport processes in a brittle–ductile intrusive model of the Taupo volcanic zone, New Zealand. J. Volcanol. Geotherm. Res. 84, 61 – 72. Whiteford, P.C., Graham, D.J., Risk, G.F., 1996. Thermal activity beneath Lake Tarawera, New Zealand, outlined by temperature measurements. Proc. 18th New Zealand Geothermal Workshop, pp. 255 – 260. Wilson, C.J.N., Houghton, B.F., McWilliams, M.O., Lanphere, M.A., Weaver, S.D., Briggs, R.M., 1995. Volcanic and structural evolution of Taupo volcanic zone, New Zealand: a review. J. Volcanol. Geotherm. Res. 68, 1 – 28. Wilson, C.J.N., 1996. Taupo’s atypical arc. Nature 37, 27 – 28. Wood, C.P., 1992. Geology of the Rotorua geothermal system. Geothermics 21 (1/2), 25 – 42. Wood, C.P., 1996. Basement geology and structure of TVZ geothermal fields, New Zealand. Proc. 18th New Zealand Geothermal Workshop, pp. 157 – 162. Wooding, R.A., 1978. Large-scale geothermal field parameters and convection theory. N.Z. J. Sci. 21, 219 – 228.