Heat Output from Spreading and Rifting Models of the Taupo Volcanic Zone, New Zealand

A conceptual model of the Taupo Volcanic Zone (TVZ) is developed, to a depth of 25 km, formed from three constant density layers. The upper layer is formed from eruption products. A constant rate of eruption is assumed, which eventually implies a constant rate of extension, and a constant rate of volumetric creation in the middle and bottom layers. Tectonic extension creates volume which can accomodate magmatic intrusions. Spreading models assume this volume is distributed throughout the whole region, perhaps in vertical dykes, whereas rifting models assume the upper crust is thinned and the volume created lies under this upper crust. Bounds on the heat flow from such magmatic intrusions are calculated. Heat flow calculations are performed and some examples are provided which match the present total heat output from the TVZ of about 4200 MW, but these either have extension rates greater than the low values of about 8 ± 4 mm/a being reported from GPS measurements, or else consider extension rates in the TVZ to have varied over time.


Introduction
The majority of geothermal activity in New Zealand is concentrated in the Taupo Volcanic Zone (TVZ). Hochstein (1995) has emphasised that the long term eruption rate of rhyolites from the TVZ is the highest on earth for a volcanic arc setting, and also that the total crustal heat transfer is arguably the highest on earth for an arc setting. The geological, geochemical and geophysical properties of the TVZ have been summarised in a Special Issue of Journal of Volcanology and Geothermal Research (Simmons and Weavers, 1995) on the TVZ. Large-scale reservoir engineering properties of the TVZ are summarised by McNabb (1975), and in the book by Elder (1976).
The location of the 23 known geothermal fields in the TVZ are plotted in Fig. 1. At least six different models (spreading, rifting, plasticity, hotplate, froth, and transtensional) have been proposed to describe some of Initial average TVZ width 9 km z Depth z 1 Lower infill surface depth z 2 Depth of upper Layer 3 surface z 1∞ Asymptotic thickness of Layer 1 4 km z 2∞ Asymptotic thickness of Layer 2 7 km z 3∞ Asymptotic thickness of Layer 3 14 km z 3 Depth of mantle 25 km z bd Brittle-ductile depth 8 kṁ A 1 Rate of increase of Layer 1 1.5 × 10 −6 m 2 /ṡ A 2 Rate of increase of Layer 2 3 × 10 −6 m 2 /ṡ A 3 Rate of increase of Layer 3 10 × 10 −6 m 2 /s C the special properties of the TVZ (see Weir, 1998a). The large scale physics of the TVZ have been summarised by McNabb (1975).

Eruptive Rates and Time-scales
The average magma flow from all of the TVZ is of the order of 1.5 m 3 /s, formed from roughly 0.3 m 3 /s (Wilson et al., 1995) of extrusive magma and roughly 1.2 m 3 /s (Wilson et al., 1984) of magma driving the geothermal fields in the TVZ. The highest rate of eruption reported (Wilson, 1993) for the TVZ occurred in the last 50 -70 ka, with time-integrated rates of 3.5 km 3 /ka for Okataina and 6.5 km 3 /ka for Taupo. While these are the highest individual rhyolitic time-averaged eruption rates in the world, to account for over 15,000 km 3 of rhyolites over 1.6 Ma (Wilson et al., 1995) requires that the present rate of rhyolitic eruptions has continued for the full 1.6 Ma of TVZ's existence. Because of these reasons, a constant rate of erupted volcanic volume in the TVZ is assumed, because present data cannot warrant more detailed information.

Assumed Depth Structure in the TVZ
Analysis of the seismic velocity structure under the TVZ indicates a three layer structure (Stern and Davey, 1987;Stern, 1987). An idealised sketch of these three layers is given in Fig. 2. The upper 2.5 km of about 3.0 km/s material is formed mostly from pyroclastic infill. The region between about 2.5 km to 15 km depth is occupied by 5.5 -6.1 km/s material. Adjacent greywacke has a similar seismic velocity. Material between depths of 15 and 25 km has a seismic velocity of 7.4 -7.5 km/s material. This region appears to be greater in width than the TVZ region above it. Mantle material adjacent to the TVZ has a seismic velocity of about 7.6 km/s.
The crust under the TVZ is perhaps only 15 km thick, whereas elsewhere is it typically 25 km thick. A major challenge for any mathematical model is to explain how the crust under the TVZ has thinned from 25 km to 15 km in the last 1.6 Ma. This is modelled in this paper, thereby extending the previous static model of Weir (1998).

Transient Spreading Model
The aim of this section is to develop approximations to the evolution of the three regions in Fig. 2 for a transient speading model. If z is depth, then z 1 is the base of the pyroclastic infill region, z 2 the top of the mantle rise, and z 3 the base of the mantle rise, assumed fixed at 25 km depth. The height of the upper (land) surface in the TVZ region will alter with time, but this surface (z = 0) is assumed to be fixed in space, because its motion is small relative to the motion of the other surfaces, and little is known about its motion. Cole (Fig. 14, 1990) suggests the middle layer (z 1 ≤ z ≤ z 2 ) region is formed by partial melting of volcanic and plutonic rocks, and high-level reservoirs of rhyolitic magma. The upper layer (0 ≤ z ≤ z 1 ) is the infill layer, although this region will comprise of both infill, mostly from rhyolitic volcanism, and intrusions. For example, a diorite intrusion was discovered at Ngatamariki (Wilson et al., 1995) in the infill region. The lower region (z 2 ≤ z ≤ z 3 ) is the mantle region, although it is possible that its composition may differ from the mantle to the west, due to the incorporation of greywacke and earlier intrusions into the rising mantle plume, and to the entry of magma from below.
The pressure at a depth of 25 km (z = z 3 ) is assumed constant. Then if where it is also assumed that the initial density of the rock is the density of the middle layer. These two assumptions are not exactly true, as a pressure difference above lithostatic is needed to drive rock upwards. This pressure difference may be relatively small, given the transient time scale for such processes is of the order of 0.01 Ma (Weir, 1998), whereas the rifting processes in the TVZ have been occurring for 1.6 Ma. Also, mixing of greywacke and igneous intrusions, for example, will alter the density in the middle layer.
Rearranging (1) shows that and so the thickness of pyroclastic infill is proportional to the height of the mantle rise. In Fig. 2 the width of the pyroclastic infill is assumed proportional to the width of the mantle rise, so (2) also implies the volume of total extruded volcanics in the TVZ is proportional to the volume of the mantle rise. This theoretical argument is then analogous, though different, to the observation that there is an approximate proportionality between magma volumes and extrusive volumes at some locations. It is shown below that the transient spreading model also has the property of a constant ratio between geothermal volumetric magma flows and erupted volumes.
The assumption above that the density of rock remains constant in each layer imposes strong constraints on the nature of the intrusions in the TVZ which can be considered in such a model. If a significant amount of rhyolite is intruded into the upper layer, then the average density there could vary over time. Below it is assumed that geothermal heating mostly results from magmatic intrusions into the middle layer, because most of the heat associated with the formation of the upper layer will be lost to the atmosphere on being deposited onto the upper surface after eruption.
Once material in the mantle rise region reaches the mantle-greywacke interface, the greywackes will be both contaminated and mobilised, and a fraction δ moves upwards, with the cooled remnant and contaminated mantle material tending to form the rising mantle plume, replacing volume lost by the mobilised greywacke.
If the average width of the TVZ is w, and this increases with time, starting from an initial average width of w 0 , then the increase in incremental mass (per unit length of the TVZ) incorporated into the TVZ region from the rise in the upper surface of the mantle plume equals the increase in mass in the middle and upper layers in the TVZ.
The factor f has been included in (3) because the average width of the mantle rise in layer 3 appears to be about 1.5 -1.7 times greater than the width of the middle layer (Stern and Davey, 1987). However, f could also be about one (Cole, Fig. 14, 1990), if the mantle plume is moving upwards, roughly parallel to the descending oceanic plate. The constant factor δ in (3) equals the fraction of Layer 3 mobilised into Layers 1 and 2. The left hand side of (3) is the mass (density ρ 2 ) initially outside, but about to be covered by Layer 3.
From (4), the final thickness of Layer 3 must be greater or equal to its present thickness, and using the density values above, It remains to determine f δ and w 0 , and the dependence of w on time. This is done by choosing a constant eruptive rate, and then considering constraints on f δ.
From Fig. 2, a constant rate of eruption is assumed, where (2), (4) and (5) have been used. From (7), w must increase linearly with time. The volume in the middle layer also increases linearly in time, since from (2), (4) and (7), where and the inequality in (9) follows from that in (6).
The left hand side of (8) is the increase in volume of the middle layer (per unit length), which must represent a flow of magma from the mantle rise region. The inflow rate of this magma is then a constant, since (8) shows a linear increase with time. Since Rt is the erupted rate of magma (per unit length) onto layer 1, R GV is the ratio of the volume of magma flow into the middle layer (which presumably fuels the geothermal fields in the TVZ) to the volume of magma extruded in the TVZ. R GV is independent of the initial width of the TVZ, and is less than 5, if the model is to be consistent with the present rise of the mantle under the TVZ.
The cumulative volume of magma in the middle layer per unit length of the TVZ is R GV Rt, in (8). Then the rate of mass flow into the middle layer, per unit volume, is which decreases monotonically with time.

Transient Heat Flow Model in the Ductile Layer
It is widely accepted that fracture surfaces in rocks with a temperature of over about 400 • C tend to creep together (Turcotte and Oxburgh, 1972). Such rocks are called ductile. The removal of fractures in ductile rocks dramatically reduces their permeability, and in this paper, it is assumed that convection of meteoric groundwater is insignificant in ductile rocks. It has been suggested (Bibby et al, 1995, McNabb andMcKibbin, 1998) that rocks become ductile in the TVZ at a depth of below about 8 km. It is assumed that any magma entering the brittle region will be cooled by circulating meteoric water to under 400 • C, and the heat transported to the surface. The equation describing conductive heat transport in the ductile region (where convective heat transport by meteoric water is assumed to be insignificant) in the presence of magmatic intrusions is where C is the heat capacity of the rock, T is temperature, t time, K thermal conductivity, T m the molten temperature of the magma, and L is the latent heat of the magma. In all calculations, K = 2 W/mK (Clauser and Huenges, 1995). Equation (11) is solved subject to an initial temperature distribution which is linear with depth, with a temperature of 400 • C at the brittleductile surface (z = 8 km), and a temperature of T m at the surface z = z 2 .
The mass flow per unit volume term, in the transient calculations, was taken as which corresponds to the TVZ increasing in width at rate v over a fixed height, where ρ 2 is the density of Layer 2, of width w. Here we assume because (13) implies maxima in volcanic activity every .1Ma, and a velocity of opening of v ± 10 mm/a. The corresponding total surface heat flows are shown in Fig. 3. The solid lines correspond to a constant opening rate of 13 mm/a, and the broken lines correspond to (13) The contribution to the energy flow from intrusivesĖ bd is constant, from (4), (5), (2) and (11), because the rate of volume creation above a fixed depth in Layer 2 is constant. In (14), z bd is the fixed brittle-ductile depth of 8 km. These values show that both conduction (1800 MW) and intrusives (1500 MW) contribute about equal energy flows at present, although earlier the majority of the heat flow resulted from intrusives. These heat flow calculations show that the variable opening rate (of 0.1 Ma period) does not significantly affect the total heat flows. If about 1000 MW of energy are added from an andesitic water flow (Weir, 1998), then total heat flows similar to those occurring in the TVZ (4200 MW) result. However, to achieve this result, it has been necessary to use an extension rate (13 mm/a) considerably greater than is measured at present (8 mm/a).

Transient Rifting Model
Chemical and isotopic analyses suggest (Graham et al., 1995) that TVZ rhyolite appears to be mostly of andesitic composition, with only a small amount of greywacke present. This suggests that the rock originally occupying Layer 3 was mostly andesite, as suggested by Cole (1990), and this requires an immense amount of andesite to initially underlie the TVZ, and also a large deep heat source. One way to circumvent these difficulties is by allowing the volume occupied by Layer 3 to be formed by rifting Layers 1 and 2 apart. Layer 3 would then be formed from andesite and other components, some of which would flow up into Layers 1 and 2, being slightly contaminated by greywacke, to form rhyolite. Some of this rhyolite must remain within Layer 2 to provide the heat and chemicals for the geothermal fields in the TVZ.
For simplicity, it is assumed that greywacke originally occupies the volume of Layer 3, and that rifting occurs, so that no greywacke is present in Layer 3. Then conservation of greywacke requires where is the mass fraction of greywacke in rhyolite and λ is the volume fraction of rhyolite in Layer 2. Layer 1 is assumed to be fully occupied by rhyolite, and so the c. 5% of other eruptive types have been ignored in Layer 1. Defining where θ is the volume ratio between stored rhyolite in Layer 2 to that in rhyolite. Similarly, the Layer areas increase asȦ 1 = 1.5 × 10 −6 m 2 /s,Ȧ 2 = 3 × 10 −6 m 2 /s, andȦ 3 = 10 × 10 −6 m 2 /s for f = 1.6. Thus volume is being created much faster in Layer 3 than it can be removed above, and so the displaced greywacke initially occupying Layer 3 must be moved out of the TVZ region, probably mostly under the North Island Shear Belt (NISB), which must move to the east in response to this flow of incoming rock. Thus uplift of the Axial Ranges of the NISB, erosion and or displacement is needed in this model. "Major uplift and dissection west, south and east of the TVZ between c. 1 Ma and c. 0.34 Ma" is reported by Wilson et al. (p.21, 1995).
Since this transient rifting model has the same extension rate (13 mm/a) and initial width (9 km) as the transient spreading model above, for θ = 2, the two models have identical thermal histories. The main difference between the two models then is the composition of Layer 3. The spreading model requires that about 56% of Layer 3 is greywacke, but in the rifting model, no greywacke is present in Layer 3.

Discussion and Conclusions
A three layer model between the surface and a depth of 25 km was formulated. Material erupts onto the surface, forming Layer 1, to balance the higher density material in Layer 3. The predicted heat output of the model contained a component from intrusions above 8 km which tended to produce a constant heat flow with time of about 1000 MW. Intrusives below 8 km produced an increasing conductive heat flow with time (because of the increasing surface area), which now could be about 2000 MW. Additional to this may be another 1000 MW from andesitic water.
Rising material which stops at a depth above about 8 km is cooled by circulating meteoric water to release its heat and chemicals into rising geothermal plumes. Material which stops below about 8 km depth can contribute some heat by conduction to the circulating groundwater above about 8 km depth. Additionally, andesitic water to the east will enter the geothermal plumes, adding chemicals and heat there.
A model valid throughout the 1.6 Ma of the TVZ was not found. No heat model for the basic heat source driving the TVZ was presented. However, if the mechanisms operating in the TVZ are those of spreading, or of rifting, then the simple models above showed that, to be consistent with the heat flow and seismic records, required spreading rates (13 mm/a) larger than are being currently measured (8 ± 4 mm/a). Another model suggested in this paper ignored the early formation of the TVZ, and considered only the latter stage of the TVZ, which was associated with the eastward motion of the NISB. Clearly, more research is needed to clarify these issues.

Call for Papers
Space dynamics is a very general title that can accommodate a long list of activities. This kind of research started with the study of the motion of the stars and the planets back to the origin of astronomy, and nowadays it has a large list of topics. It is possible to make a division in two main categories: astronomy and astrodynamics. By astronomy, we can relate topics that deal with the motion of the planets, natural satellites, comets, and so forth. Many important topics of research nowadays are related to those subjects. By astrodynamics, we mean topics related to spaceflight dynamics.
It means topics where a satellite, a rocket, or any kind of man-made object is travelling in space governed by the gravitational forces of celestial bodies and/or forces generated by propulsion systems that are available in those objects. Many topics are related to orbit determination, propagation, and orbital maneuvers related to those spacecrafts. Several other topics that are related to this subject are numerical methods, nonlinear dynamics, chaos, and control.
The main objective of this Special Issue is to publish topics that are under study in one of those lines. The idea is to get the most recent researches and published them in a very short time, so we can give a step in order to help scientists and engineers that work in this field to be aware of actual research. All the published papers have to be peer reviewed, but in a fast and accurate way so that the topics are not outdated by the large speed that the information flows nowadays.
Before submission authors should carefully read over the journal's Author Guidelines, which are located at http://www .hindawi.com/journals/mpe/guidelines.html. Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking System at http://mts.hindawi.com/ according to the following timetable: