Geothermal evolution of an intruded dike in the rift zone of Kilauea volcano, Hawaii from VLF and self-potential measurements

Journal of Volcanology and Geothermal Research 302 (2015) 64–80

Contents lists available at ScienceDirect

Journal of Volcanology and Geothermal Research journal homepage: www.elsevier.com/locate/jvolgeores

Geothermal evolution of an intruded dike in the rift zone of Kilauea volcano, Hawaii from VLF and self-potential measurements Paul M. Davis ⁎ Department of Earth, Planetary and Space Sciences, UCLA, United States

a r t i c l e

i n f o

Article history: Received 19 March 2015 Accepted 15 June 2015 Available online 19 June 2015 Keywords: Dike Hydrothermal Electromagnetic

a b s t r a c t Self-potential (SP) and VLF measurements were made in 1973, 1975, 1995, 1997 and 2012 across a basaltic dike that intruded into the Koae fault zone of Kilauea volcano, Hawaii in May 1973. The SP anomaly remained strong throughout. In 2012 it was at about 60% of the strength it had in 1973. In contrast, the VLF anomaly, though diminished, was still observable in 1995/1997, but by 2012 it had disappeared. A hydrothermal dike model, with parameters calibrated by modeling the solidification of Kilauea Iki lava lake, is used to calculate temperatures and conductivity variation. Following Jaeger's (1957) method, we find that the time in years for a dike of width W (m) to solidify is 0.0075W2. Thus, a 1 m dike solidifies within the first few days, and after 39 years is only tens of degrees above ambient. Given the orders of magnitude difference between the conductivities of wet and dry basalt, we infer, that after solidification, the VLF anomalies were caused by induction in a localized veil of wet, hot basalt enveloping the dike, that was generated initially by condensation of steam, and subsequently by condensation of evaporated water as temperatures reduced. The conductivity anomaly persisted until the mid-nineties. By 2012, temperatures and condensation were too small for a VLF signal. The persistent SP anomaly is attributed to localized fluid disruption, with evaporation mainly at the water table and in the vadose zone. Streaming potentials are associated with evaporative circulation in the vadose zone. Next to the dike a positive potential is generated by upward flow of moisture-laden air, with a smaller negative potential on its flanks from downward infiltrating rainwater. The analysis indicates that the combination of SP and VLF measurements can characterize the evolving geothermal regime of intrusions above the water table. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Geothermal reservoirs are found above intrusions of magma such as dikes or dike swarms, which set up hydrothermal circulation generating hot water and steam from which energy can be tapped (Furumoto, 1975; Kauahikaua et al., 1980; Kauahikaua and Mattice, 1981; Jackson and Sako, 1982; Kauahikaua, 1993). The hydrothermal system has associated electromagnetic effects that can be used for exploration as well as predicting the lifetime of a reservoir (Banwell, 1970; Anderson and Johnson, 1976; Combs and Wilt, 1976; Corwin, 1976; Corwin and Hoover, 1979; Fitterman and Corwin, 1982; Kauahikaua et al., 1986; Ishido and Pritchett, 1996, 1999, 2000; Revil et al., 1999). In volcanic zones, intrusions generate electromagnetic changes that are used to monitor volcanic activity (Zablocki, 1976, 1978a,1978b; Jackson et al., 1985; Jackson and Kauahikaua, 1987; Zlotnicki and Le Mouel, 1990; Hashimoto and Tanaka, 1995; Johnston, 1997; Sasai et al., 1997) and may be indicative of stress (Zlotnicki and Le Mouel, 1990). However

⁎ Tel.: +1 310 825 1343; fax: +1 310 825 2779. E-mail address: [email protected].

http://dx.doi.org/10.1016/j.jvolgeores.2015.06.007 0377-0273/© 2015 Elsevier B.V. All rights reserved.

the origins of these signals has been controversial; in particular questions have been raised as to whether SP voltages are generated by hydrothermal circulation beneath the water table or in the vadose zone; or whether high conductivities are due to molten basalt or saline pore fluids (Jackson and Keller, 1972; Kauahikaua, 1993). Our control experiment on the life cycle of a Hawaiian dike seeks to address these questions on a structure with a well-constrained geometry and initial thermal conditions. We report very low frequency (VLF) electromagnetic sounding and self-potential (SP) measurements over a period of 39 years that span two cooling regimes associated with a dike intrusion into Kilauea volcano in 1973. The first is rapid and involves efficient heat transfer through degassing and internal convection of the magma until the dike solidified. The second is orders of magnitude slower involving conduction and evaporation in the vadose zone, and hydrothermal circulation below it in the water table. We used a combination of thermal and electromagnetic modeling to invert the data over the geothermal lifetime of this intrusive body. The results lead us to support the suggestion by Zablocki (1978b) that the SP signal is generated in the vadose zone, rather than beneath the water table. We also find that after solidification, the VLF signals are generated in the vadose zone by concentration of condensed water about the hot dike.

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

VLF and SP signals have long been measured above volcanic intrusions, with the former often attributed to the high conductivity of molten basalt and the latter to streaming potential caused by electrokinetic (EK) effects (Zablocki, 1976, 1978a). For lava tubes VLF anomalies are clearly related to the high conductivity of molten basalt (σ = 1 S/m) since it can be viewed through skylights along its path (Kauahikaua et al., 1998). Thus it has been natural to interpret dike anomalies as arising from molten basalt as well (Zablocki, 1976, 1978a). Geodetic inversion of surface deformation associated with intruding dikes, and associated seismicity, indicate that intruding dikes may extend from near the surface to several km deep in the rift zones in Hawaii (Dvorak et al., 1986; Yang et al., 1992; Cervelli, et al., 2002). Fits of half-space dike models to VLF signals that assume the dike is represented as a sheet of magma in a low conductivity medium (σ = 10−4 S/m) give short dikes, tens of meters long (Zablocki, 1976, 1978a), which suggests that the induction effects concentrate near the surface. On average, exhumed dikes in Hawaii are about 1 m thick (Walker, 1986, 1987, 1988; Spengler and Walker, 2012), a value that corresponds with those of the geodetic inversions. Dike widths show a linear increase with distance from the eruption center obtaining values up to 3 m (Spengler and Walker, 2012). The duration of most flank eruptions on Kilauea has been shown to correspond to estimated time to solidification (Spengler and Walker, 2012) with average values of about 5 days. SP anomalies in geothermal regions have most often been attributed to EK effects associated with hydrothermal circulation (Zohdy et al., 1973; Anderson and Johnson, 1976; Corwin, 1976; Corwin and Hoover, 1979; Fittermann, 1979; Ishido, 1981, 1989; Ishido and Mizutani, 1981; Ishido et al., 1983; Zlotnicki and Le Mouel, 1990; Hashimoto and Tanaka, 1995; Ishido and Pritchett, 1996, 1999; Sasai et al., 1997; Revil et al., 1999). Chemo-electric and thermoelectric mechanisms have been proposed but are thought to be smaller. Laboratory experiments find that fluids streaming in rock generate significant self-potentials (Ishido and Mizutani, 1981; Marsden and Tyran, 1986; Morgan et al., 1989; Antraygues and Aubert, 1993) with positive potentials in the direction of flow. Because negative ions attach to wall rock, positive ions flowing in the water are displaced in the direction of the flow until an equilibrium is reached between the electrostatic restoring force and the viscous drag. Positive self-potential anomalies, with values of 1 V or more, are ubiquitous in geothermal zones. While streaming potentials from EK effects are thought to be the most likely mechanism, it has not been clear what the relative contributions of hydrothermal circulation beneath the water table, and flow in the sub-saturated vadose zone are. Various hydrothermal flows have been modeled (Yasukawa et al., 1993; Ishido and Pritchett, 1996, 1999, 2000). Flowing water containing ions such as dissolved salts is an effective generator in a hydrothermal circulation (Ishido and Mizutani, 1981). However, flowing steam and humid air also generate significant coupling coefficients (Antraygues and Aubert, 1993; Johnston et al., 2001). Dry steam is not effective, but wet steam is (Antraygues and Aubert, 1993). EK effects in gases depend on water content and mode of evaporation. Johnston et al. (2001) conducted experiments in which slugs of water were evaporated in 200 °C rock. They showed that rapid evaporation of water caused positive charges to separate into the evaporated droplets carried by the flow, leaving negative charges behind, generating tens of millivolts positive potentials in the direction of flow. They refer to this mechanism as rapid fluid disruption (RFD) and suggested it could explain the positive polarities of potentials measured on numerous volcanoes (Zablocki, 1976; Hashimoto and Tanaka, 1995; Sasai et al., 1997). They also note that charge separation can occur in low temperature situations where flowing gas evaporates water from water saturated rocks. Evaporation of pore water also decreases conductivity, a necessary condition for the long-term persistence of electric fields, which in high conductivity regions would otherwise short out. Once RFD has taken place, the dipolar separation of charge by flow exhibits an EK effect in which an equilibrium between mechanical and electrostatic forces gives rise to electric polarization.

65

Infiltration and evaporation of rainfall can also give rise to SP potentials that change from negative upwards to positive upwards as evaporation takes over from infiltration (Thony et al., 1997; Darnet and Marquis, 2004). Near the summit of Kilauea volcano the water table was found to be 488 m deep at the location of a research borehole (Keller et al., 1979). It is estimated to be at about 230 m near the location of the intruded Hi'iaka dike (Jackson and Kauahikaua, 1987). The annual rainfall in Hawaii is about 200 cm. Water percolating downwards into the vadose zone is thought to generate a negative SP anomaly due to the EK effect. However in thermal areas relatively positive anomalies are observed where it is inferred that upwelling currents negate the effects of downward percolating rainwater and generate a positive voltage (Zablocki, 1978b). The SP measured on the surface depends on the geometry of the flow and the variation in material properties. Nourbehecht (1963) and Fitterman (1978) showed that a pressure source in a homogeneous half-space generates no SP anomaly on the surface. This occurs because the SP is proportional to pressure and, since the pressure is zero on the surface, the surface is an equipotential. Also hydrothermal circulation in a uniform porous half space gives zero SP on the surface. The current source of the potential is the divergence of the pressure gradient (Sill, 1983). Since velocity is proportional to pressure gradient through Darcy's law, for an incompressible fluid the divergence of velocity is zero, and so then is the divergence of pressure gradient, resulting in a zero source. The zero external fields are effectively a result of pole cancelation. To generate a non-zero external field requires some mechanism to give rise to uncompensated poles. The symmetries associated with the pole cancelation can be broken by non-uniform coupling coefficient, such as might happen across a boundary separating different materials (Fitterman, 1978) or its dependence on temperature or PH (Ishido and Mizutani, 1981). Anisotropic permeability can also break the symmetry. For example, the fissure zone above an intruded dike will channel the flows upwards resulting in a positive dipole. Evaporation and flow of water vapor in a porous medium will also give rise to uncompensated poles, with negative poles at the evaporation surface and positive poles where the vapor either disperses or condenses. 1.1. Vadose zone models Zablocki (1978b) was the first to argue that the SP anomalies in the volcanic zone in Hawaii are generated mainly in the vadose zone. He sought one mechanism to explain the long wavelength negative anomalies that exhibit an inverse correlation with topography, and the shorter wavelength relatively positive anomalies above subsurface localizations of heat. In a homogeneous region, depth to vadose zone correlates with topography. The electrical polarization associated with downward flow gives rise to a negative voltage on the surface, and positive at depth proportional to the distance traveled to the water table. For example, Zablocki (1978b) reported the difference in selfpotential measured between the water table (488 m deep) and the surface in the bore hole drilled near Kilauea's summit was − 655 mv (−1.35 mv/m). In contrast, in thermal zones the downward meteoric flow would tend to be blocked by upward flows of buoyant heated gas containing water vapor, resulting in local positive anomalies compared with the surroundings. This explained the inferred shallow depths (tens of m) of the SP anomalies in regions where the water table is deep but magmatic intrusion shallow. The vadose zone model was (‘partially’) supported by Jackson and Kauahikaua (1987) who compared depth to the water table inferred from vertical electromagnetic soundings (VES) with depths from SP and concluded SP was better correlated with vadose zone thickness than topography. The dikes of Kilauea's rift zones tend to dam the water table, resulting in heterogeneous variation in vadose zone depth compared with topography. This allowed the authors to differentiate between SP proportional to topography or proportional to vadose

66

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

zone depth. However, different constants of proportionality (up to several mv/m) and depths, suggested other factors must be taken into account, such as lateral variation in streaming potential, or contributions from beneath the water table. VLF has been used in Hawaii to detect shallow subsurface intrusion of magma (Zablocki, 1976, 1978a; Kauahikaua et al., 1998) associated with intruded dikes or lava tubes. The frequency of VLF waves (~20,000 Hz) gives rise to a skin depth of about in the highly resistive rock (conductivity σ = 10−4 S/m) near the surface. Magma is orders of magnitude more conductive, such that induced eddy currents generate significant secondary fields at the surface. Laboratory measurements of the conductivity of molten basalt give values around 1 S/m, but on solidification the conductivity of dry basalt reduces by many orders of magnitude (e.g., to 10−10 S/m to 10−12 S/m at room temperature, Presnall et al., 1972; Schloessin, 1976). The higher values of regional conductivity near the surface (∼ 10−4 S/m) are probably due to negligible pore fluids in the vadose zone. In addition, both wet and dry basalt exhibit a significant temperature variation described by an Arhennius relation (Presnall et al., 1972; Schloessin, 1976; Olhoeft, 1977, 1981; Scarlato et al., 2004). The VLF signal is roughly proportional the product of conductivity contrast and thickness (for very high conductivity this effect saturates). Given that dikes are typically 1 m thick, while molten they should give significant anomalies. But on cooling the molten dike effect will rapidly disappear and the dike even becomes a conductivity hole, albeit a thermal radiator. Thus an initially 1200 °C dike, that intersects the water table will cool rapidly, but while temperatures are above the boiling point of pore water (BP), it will generate steam, which, as it rises to the surface condenses to fill a localized region of the vadose zone with water, with a shape similar to the dike itself. As the heat spreads out, the hot, wetted region below BP can then have a significant thickness– conductivity-contrast relative to the cooler surroundings. For example, we estimate the half-width of the heat anomaly after 22 years to be 62 m. A dike-zone with 1/62 the conductivity of a 1 m molten dike will give similar VLF effects. In contrast, an interior zone, including the dike, will be dry enough for SP effects to generate. External to this, the region where water condenses and trickles down will be conductive. We suggest here that this veil-like configuration (see summary Fig. 16) is the origin the VLF and SP signals we see. 2. Hi'iaka eruption On May 5, 1973 a dike intruded into the upper crust of Kilauea volcano associated with the eruption of Hi'Iaka and Pauahi craters along Kilauea's east rift zone (Klein et al., 1987; Tilling et al., 1987). It generated a surface fissure 100 m long that erupted magma WSW of Hi'iaka crater. Geophysical measurements suggested it continued underground in a WSW direction for a further 1.5 km (Fig. 1). The geophysical measurements included SP and VLF electromagnetic measurements made ten days after the eruption by Charles Zablocki (1976, 1978a), a geophysicist at the USGS Hawaiian Volcano Observatory (HVO). He measured a series of 6 SP-VLF profiles that defined an elongate anomaly of about 1 V amplitude, and 40% tilt, respectively, which suggested that the dike had intruded into the Koae fault system for a further 1.5 km without breaking the surface (Fig. 1). We refer to this intrusion as the Hi'iaka dike. The VLF measurements were made using the transmitter NLK (18.6 kHz) in the state of Washington. He also measured the temporal change in VLF across the Hi'iaka dike (digitized versions of his plots are shown with measurements in May 15 and November 11, 1973 and February 20 in 1975.) A further reading in January 1976 was similar in shape to the 1975 measurement, but used the VLF transmitter NPM (23.4 KHz) at Lualualei on the island of Oahu, Hawaii. It had a slightly larger amplitude, probably due to the different frequency and different angle of the VLF field relative to the dike. Our subsequent VLF readings (1995,1997, 2012) were made with NPM. Thermal modeling shows that the dike should have solidified and

cooled in the first few days. The large VLF anomalies over the subsequent years required an explanation that does not involve molten basalt. Zablocki (1976) also measured changes in SP across another intrusion in 1974 showing that as it cooled the height of the anomaly decreased, but the width remained the same. The half-width of geophysical anomalies usually provide an estimate of the depth of their sources; in that case the depth to the top of the hot rock associated with the intrusions. But its constancy in time as the intrusion cooled, led him to infer that half-width was more affected by the properties of the hydrothermal system, such as depth to water table and permeability. For a dike cooling from the top, the half-width would be expected to increase as the top cooled off. 3. Experiment In this work we extended the original measurements to track the evolution of both VLF and SP along one of the profiles (profile 3, Fig. 9, Zablocki, 1978a) measured by Zablocki across the Hi'iaka dike (Fig. 1). In order to compare our measurements with his, we digitized his curves. For the starting point of our surveys we chose a point 0.72 miles along the Hi'iaka road W of Chain of Craters road. GPS was available only for the last survey, and the end points were [start (19o22′8.5 ", 155o14′35.6 " end 19o22′31.1 " 155o14′47.9 ")]. For the earlier surveys we used a compass and tape measure, which will not be as accurate. The measurements included: 1. Self-Potential, SP, which involves inserting two 1-inch wide nonpolarizing Cu/CuSO4 electrodes in the ground to a depth of 6 in., separated by about 20 m and measuring the voltage between them. A survey about 1 km long was performed. A high impedance (10 MOhm) Fluke voltmeter was used to measure voltage. The points were leap-frogged after each reading in which the rearmost electrode was carried forward to become the foremost electrode. In this way any differential contact potentials should cancel out. 2. Very Low Frequency VLF electromagnetic sounding involves using a portable VLF (GEM) meter to measure the tilt and phase of the polarization ellipse of radio waves from the VLF transmitter on Oahu. Both are affected by underground structures containing hot rock or magma. The meter is carried at eye level and the data are recorded along with GPS coordinates (or tape measure). Because the geological structure has an intrinsic inductance depending on its shape and conductivity, the secondary and primary fields are non-parallel and out of phase, and combine to give a tilted elliptical signal. A VLF meter uses a phase shifter on the outputs from two orthogonal detection coils to eliminate the ellipticity, which allows the tilt to be determined by a nulling method, in which the observer rotates the instrument, until the long axis of the detector coil is at right angles to the linearized field. The observer then records both tilt and ellipticity. The tilt is given as 100 times the tangent of the tilt angle from the horizontal, expressed as percent. It is approximately the percentage of 45°. 4. Results Following on from Zablocki's measurements in the period 1973– 1975, our first survey was made in May 3 1995. The approximate location of his 1973 line 3 was obtained from his published maps (Fig. 1) and an SP and VLF profile was measured (Figs. 2,3) with the surprising result that the SP field was nearly as strong as it was in 1973, while the VLF field though decayed, was still observable at about 30% of the original values. There was no evidence of steam or open fissures in the region. The region near the peak of the anomaly was open, having little vegetation. The surveys were repeated in October 21 1997 and most recently in December 18 2012. The data show decay with tme in the VLF field (Fig. 3) such that the anomaly is unrecognizable in 2012. However,

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

67

19o 22' 31.1"N 155o 14' 47.9"W

19o 22' 8.5''N

155o 14' 35.6"W

Fig. 1. Location of VLF and self-potential measurements (dashed line) relative to the inferred location of the Hi'iaka dike (solid line). (modified after http ://hvo. wr. usgs. gov/gallery/kilauea/ erz/hiiaka−erupt. html; Last accessed June 2015).

the SP anomaly remained strong, but had broadened, and decreased by about 40% (Fig. 2). In addition, the region had become dense vegetation. The time histories of the peak-to-peak VLF and maximum SP values are shown in Fig. 4. 5. Analysis 5.1. Thermal modeling Our objective in the following analysis is to construct a hydrothermal model to explain the temporal behavior seen in the SP and VLF measurements. While the dike is molten the VLF field is generated by the conductivity of the magma which at about σ = 1 S/m (e.g., measurements by Presnall et al., 1972). But a 1 m wide dike will

solidify in days. Thus, the VLF anomalies, seen at 10, 190, 638 days after the eruption and up to 22–24 years later, must be a sub-solidus effect. As the dike cools and temperatures become lower than BP, the relevant conductivity is that due to hot wet basalt and that can be many orders of magnitude greater than that of dry basalt. The temperature variation of a one-dimensional cooling dike was examined by Jaeger (1957). The effect of latent heat of solidification was taken into account by modifying the specific heat for parts of the dike above the solidus. The heat equation was then solved numerically with contrasting diffusivities for molten and solid rock. We recalculated the Jaeger curves to confirm that our numerical scheme gave his results. We then determined the parameters of his model that fit the time history of the solidifying crust of Kilauea Iki lava lake (Richter and Moore's, 1966; Fig. 8).

Fig. 2. Self-potential measured across the Hi'iaka dike as it cooled over a period of 39 years. The inferred hydrothermal circulation broadens with time and decreases in amplitude by about 40% in contrast with the VLF anomaly that decreased by over an order of magnitude.

68

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

Fig. 3. VLF tilt measurements above the Hi'iaka dike between 1973 and 2012 showing the effect of cooling.

Jaeger approximated the effect of latent heat, L by modifying the specific heat, Cp by Cp ¼ Cp þ

L T 2 −T 1

ð1Þ

where T2 is the intrusion temperature and T1 the solidus. Then the 1-D diffusion equation is solved numerically with corresponding different diffusivities used for molten and solid material. In an eruption in 1959 Kilauea Iki lava lake was filled with 135 m of magma (Richter and Moore, 1966; Barth et al., 1994). Over the next three years the thickness of the solid crust, which reached 14 m, was measured seven times by drilling. Cooling was thought to be due to a combination of thermal diffusion and evaporation of rainfall. Following Jaeger, we also used a 1-D model, since the dimensions of the cooling crust are small relative to those of the lake. We fit the model to the time variations of crustal thickness (Richter and Moore, 1966) assuming the initial temperature of the magma was Tm = 1150oC the solidus Tm = 950 °C with the surface boundary condition T = 25 °C (Fig. 5). Thermal diffusivities for molten and solid material which fit the data are D = 1.0 × 10−6m2/s and D = 0.39 × 10−6m2/s respectively. However an average of D = 0.5 × 10−6m2/s gives a comparable fit. Geophysical measurements showed that the crust immediately above the magma had very low conductivity (Barth et al., 1994); consistent with it being dehydrated by the heat from the magma, but closer to the surface the conductivity was higher. We turn this model on its side to approximate a dike in the vadose zone. Jaeger noted that the cooling curves are self-similar for dikes of any width. He gave the time in years for solidification of a molten half-space to distance W as ts = (W/2)2/(4λ2/D1), D1 is the diffusivity. For Kilauea Iki we used his Eqs. (3-5) to find λ = 0.37 assuming L = 80 cal/g. We obtain ts = 0.015 W2 for Kilauea Iki, which is similar to his results (i.e., after 3 years a thickness of 14 m). For the dike case, since heat leaves in both directions, we divide the half-space equation by 2. Then for a dike time to solidification is t s ¼ 0:0075 W2 where W is the width in m and ts is time in years.

ð2Þ

As has been mentioned, exhumed dikes in Hawaii are mostly W = 1 m, and on rare occasion have W = 3 m been seen. Thus the Hi'iaka dike will have solidified in 2-18 days after emplacement. The VLF measurements after 10 days, and certainly after 190 and 638 days and later were most likely above a solidified dike, and as we shall see at 10 and 190 days dike temperatures are above 100 °C, and so it would also be dry. The temperature variation at the times of each survey of a 1-D cooling 1 m dike model is shown in Fig. 6. For the sub-solidus variation the analytic solution of Carslaw and Jaeger (1986) was used with the aforementioned average diffusivity of k = 0.5 × 10−6m2/s.

T ðx; t Þ ¼

T0 2

( er f

! !) ðD=2−xÞ ðD=2 þ xÞ pffiffiffiffiffiffiffiffiffiffiffi þ er f pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðk1 t Þ ð2 ðk1  t Þ

ð3Þ

where T0 is the intrusion temperature and the surroundings are initially at T = 0 °C. A background ambient temperature T = 25 °C is superposed on the resulting distribution. The temperature variation (Fig. 6) shows that after 10 days the field had diffused 1 m away from the dike with a maximum (sub-solidus) temperature of 500 °C. By 1995 it had dropped to 42 °C spread over 62 m. If the dike was 2 m thick rather than 1 m the temperatures can be multiplied by 2. However as well as diffusion a full description of the thermal field about a cooling dike needs to take hydrothermal circulation into account, which we examine in the next section. 5.1.1. Hydrothermal modeling The results of the Jaeger model of the last section show that the therpffiffiffiffiffiffiffiffiffiffi mal field diffuses at about 1.5 m/ year away from the dike. Thus after one year its width is 3 m. Since heat is conserved, the temperature of a 1 m dike is proportionally reduced by diffusion by a factor of 6 to be about 200 °C, i.e., above BP (which in a hydrothermal system depends on depth, ranging from 100 °C at the surface to 220 °C at the water table, in our case). Then the VLF signal is likely due to wet basalt as steam condenses in the surrounding rock. The region right next to the dike above BP will be dry and have a very low conductivity, but the

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

69

Fig. 4. Double exponential decay fitted to the differential VLF tilts from Fig. 3 and linear decrease in the maximum of the self-potential anomaly (Fig. 2). The initial rapid decay is related to the solidification and water saturation of the vadose zone near the dike as steam condensed, whereas the final decay is due to conductive cooling of the solid body. The initial decay is governed by internal convection in the molten dike whereas the latter is by thermal conduction resulting in an orders of magnitude larger decay constant.

return flow will generate pore water from condensing steam in the exterior region. That can give rise to a concentrated, dike-like, high conductivity zone, which produces a VLF signal. We examine this behavior using a porous hydrothermal model in this section.

Inversion of the 1973A VLF data showed the top of the dike was at about 20 m depth and it dipped to the south at 80–90° (Zablocki, 1976; Davis, 2014). The total circulation involves hydrothermal circulation beneath the water table and circulation of steam and water vapor in

Kilauea Iki Solidification 14

12

Thickness (m)

10

D1 = 1e−06 m2/s

8

D2 = 3.9375e−07 m2/s 6

4

2

0 0

0.5

1

1.5

2

2.5

3

Time (years) Fig. 5. Time variation of the depth of the crust that formed above Kilauea Iki Lava lake fit to Jaeger's 1-D numerical model with contrasting diffusivities that model the effect of latent heat of solidification. An average D = 0.5 × 10−6m2/s provides a similar fit.

70

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

T 1973

T 1995 50

Temperature (°C)

Temperature (°C)

600 500 HH= −1 (m)

400 300 200 100 −4

−2

0

2

45 40 35 30 25 −200

4

HH= −31 (m)

−100

Distance (m) T 1997

100

200

50

45

Temperature (°C)

Temperature (°C)

50 HH= −33 (m)

40 35 30 25 −200

0

Distance (m) T 2012

−100

0

100

200

45

HH= −42 (m)

40 35 30 25 −200

−100

Distance (m)

0

100

200

Distance (m)

Fig. 6. Time variation of the temperatures of a cooling 1 m dike. The half-width (HH) were used to set the half width of the SP generation zones in Figs. 14, 15.

the vadose zone. The coupled problem is beyond the scope of what is attempted here. Instead, we examine two cases separately (1) that of hydrothermal circulation beneath the water table and (2) convective circulation of air in the vadose zone. After non-dimensionalizing the problem (see Supplementary Appendix A) the convective heat flow equation (Nield and Bejan, 1992; Turcotte and Schubert, 2002) in a porous medium is given by   ^ ^ T^ ^ þ ψ ^ ^ T^ ^ þ ∇2 T T^^t ¼ Ra −ψ y x x y

ð4Þ

^ is the normalized where Ra is the porous medium Rayleigh number, ψ stream function. The first term on the right is the advection and the second the (normalized) diffusion. The relative strengths of each depend on the Rayleigh number (See supplementary material Appendix A for more details). Our model consists of a 1001 × 1001 m 2-D box with a dipping 1 m wide dike at its center beginning at a depth of 25 m from the top and extending to the bottom. The intrusion temperature was set to T0 = 1200 °C, while the rest of the medium was at T = 0. The surface was held a constant temperature T = 0 °C and the heat flux was set to zero at the sides and bottom. Eq. (4) was solved at each time step using finite differences. The result showed that the dike channels an upward flow in its vicinity that circulates outward at the top to complete the circuit (Fig. 7). The Rayleigh number for the hydrothermal case 103 is nearly an order of magnitude larger than for the air case. In both cases the diffusion term dominates. The main effect of advection was to redistribute the thermal energy making the top of the dike a little hotter and the bottom cooler than in the diffusion-only case. (An example for a vertical dike in a hydrothermal circulation is shown in the supplementary material, Figure A1). For the air circulation case the difference between an advecting model and a diffusion-only model is not significant enough to warrant the extra computing. So for the subsequent analyses we used the 2-D dipping dike diffusion model.

5.2. Basalt conductivity In order to model VLF variation we construct the conductivity crosssection as a function of time. This involves converting temperatures to conductivity along with making estimates of water content. The conductivity of basalt decreases by 2.5 orders of magnitude going from the liquidus to solidus (Presnall et al., 1972; Schloessin, 1976; Scarlato et al., 2004) and then decreases exponentially as it cools. Since these samples started as molten we assume that on solidification they are relatively dry, giving estimates of σ (T)dry. For wet basalt we use data for brine-saturated basalt of 14% porosity containing (0.1 m) aqueous solution of NaCl, (from Olhoeft, 1977) to obtain σ (T)wet. The laboratory data are plotted in Fig. (8). We see in this figure that at temperatures below melting, water-saturated rock has a much higher conductivity than dry rock at the same temperature. The dependence of conductivity on temperature T is given by σ ðT Þ ¼ a1 eð−a2 =T Þ

ð5Þ

where the temperature is in o K. We fit the a1 and a2 to the published sub-solidus data. For the wet case we obtain, σ ðT Þwet ¼ 331:5eð−3:139

Þ

1000 T

ð6Þ

whereas for dry basalt we obtain σ (T)dry. 1000 σ ðT Þdry ¼ 145eð−8:365 T ÞÞ

ð7Þ

The fits of these equations are shown in Fig. 8. At room temperature the dry case has a conductivity of 10−10S/m. Schloessin (1976) obtains 10−12S/m for his ultra-dry samples in a vacuum.

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

71

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

Fig. 7. Convective circulation about a 1 m dike. (Units are normalized fractions of 1000 m.).

Conductivity Wet Dry Basalt 2

o 727 C

227 oC

60 oC

1

0

logσ (S/m)

0 −0.5

Dry

Wet

−1

−2

−3 1000 K

−4 0.5

1

500 K

1.5

2

333 K

2.5

3

3.5

1000/TK Fig. 8. Conductivity of wet (red) and dry (blue) basalt after Olhoeft (1977) and Presnall et al. (1972) fit to σ ðTÞ ¼ a1 eð−a2 =TÞ .

0.5

72

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

5.3. VLF modeling of dike injection Initially we inverted the VLF field taken in 1973 using a half-space model (Davis, 2014) to obtain approximate dimensions of the dike. Generally much of the induced field occurs within a few skin depths of the surface. The skin depth for VLF waves ( f = 20, 000 Hz conductivity σ1 = 1 × 10−4 S/m) in Hawaii is about 360 m. Dikes in Hawaii typically have dimensions greater than several skin depths. Thus we tried a representative model of an inclined dike 1000 m long (Davis, 2014). However a much better fit was obtained with a dike 30 m long, starting at a depth of 20 m, which was a puzzle because long dikes are thought to extend km into the Earth. After numerical tests we have concluded this apparent shortness may be due to an over-simplified model, and that one has to take into account the increasing conductivity of the surrounding medium with depth. The background σ1 = 1 × 10−4 S/m applies to near-surface moderately dry basalt. The brackish water below the water table has a conductivity orders of magnitude higher. Jackson and Keller (1972) have noted that in some regions with sufficient salinity conductivity may even be similar to molten basalt at 1 − 2 S/m. Since the exponential decay associated with a 300 m skin depth has significant energy at the depth of the water table (here 230 m) its effect needs to be included. Using a finite difference (FD) program we tested the effect of the water table conductive layer on the VLF tilt field generated by a deeply penetrating dike. The finite difference program was benchmarked against the boundary element program (Hohmann, 1971; Davis, 2014) for the case of a dike in a uniform half space. The boundary conditions of the E field must match the exponentially decaying incident field at each side of the model. The Helmholtz equation ∇E = k2E is solved where k(x, z) is zero in the air layer above the surface and has values in the half space and dike corresponding to their respective conductivities and the frequency of radiation; the permeability, μ 0 is assumed to be that of a vacuum. Adequately close agreement between the halfspace models and the finite difference calculation was obtained for a 25 m dike at a depth of 25 m if the model space is 2001 × 2001 m in extent, in which the upper half of the model space is air (Fig. 9) and the lower half includes the dike. Too small a model space results in

Air σ=0

layer −4 σ=10 S/m

Dike L

layer −2 σ=10 S/m

Fig. 9. Finite difference model for the effect of a deep dike that crosses a conductive layer associated with the water table. The dark shaded region represents a dike that extends to depth. The white shaded region is the apparent size of the dike if a homogeneous half space model is used to invert the data. The effect of the conductive layer is to attenuate the incident field resulting in a tilt signal on the surface that corresponds to a significantly truncated dike. This explains why short dikes (tens of m) fit half-space VLF models.

interaction of the perturbation field with the boundary conditions, which are only accurate at infinity. For the boundary conditions for the two-layered case, corresponding to vadose zone and water table, we used Cagniard's (1953) equations for the two-layer MT problem interfaced to the air layer by matching the boundary conditions of tangential E and H fields at the surfaces. The incident E field, which is linear in the air layer continues its almost linear variation in the upper layer, and then decays rapidly in the highconductivity lower layer. It has a much lower amplitude and faster decay than the half-space exponential. The saturated zone below the water table is given a higher conductivity, (σ = 10−2 S/m), below a vadose layer of conductivity, (σ = 10− 4 S/m), with the dike having a conductivity, (σ = 1 S/m). The half-space VLF tilt curve for a short-dike 25 m deep of 25 m length is closely reproduced by a long dike in the two-layer model when the upper depth is 12 m and lower depth 1000 m, and the water table layer is at 230 m. Thus the effect of the water table is to attenuate the incident field resulting in a tilt signal on the surface that appears to correspond to a significantly truncated dike. This explains why models of short dikes (tens of m) in a homogenous half space fit VLF data, whereas geological considerations expect them to be long. Therefore, for the subsequent analyses we used the finite difference calculations. The half-space model indicated that the dike had a depth of about 20 m and a dip of 80 − 90° (Zablocki, 1978a; Davis, 2014). Given the short distance from the eruptive center at Hi'iaka crater, we fixed the width to 1 m based on the exhumed dike data (Spengler and Walker, 2012). Also the diffusivity was set at D = 0.5 × 10− 6m2/s consistent with the Kilauea Iki modeling (Fig. 5). Thus the geometry was fixed for the fit to each of the 5 sets of measurements (1973A, 1973B, 1975, 1995 and 1997). The thermal field after intrusion of a T = 1200 °C dike was calculated, and the E field was used to calculate tilt to compare with the observations. We calculate the VLF effect in the two regimes, dry, wet, by solving the thermal equations and then converting to conductivity using Eqs. 6, and 7. We found the best fit to the data had the second layer at a depth of 150 m, which, while short of the estimate of the water table at 230 m, may be the result of using just one layer to model a steady increase in conductivity. Alternatively, this could be the true depth since the VES location was some distance away. For temperatures below BP we used the Olhoeft relation (Eq. 6), while for those above BP the Presnall relation (Eq. 7). For a 1 m dike at 10 days it has solidified. Therefore the signal is mainly due to the high conductivity veil comprised of wet basalt. Temperatures below BP are shown in Fig. 10. In detail, the modeling comprised the following steps: (1) We calculated the thermal field numerically at each time step on a 1001 m × 1001 m cross-section, for a 25 m deep, dipping dike. (2) We constructed conductivity models by taking the upper 500 × 1001 cells of the temperature model and adding an air model above of the same size. We assigned background conductivities of σ = 0 in the air layer, σ = 10−4S/m in the vadose zone and σ = 10−2 S/m beneath the water table at 230 m below the surface. (3) Because the temperature field diffuses slowly many of the values are close to ambient. So much of the VLF model has background values of conductivity. For temperatures higher than ambient we used Eq. (7, 8) to convert temperature to conductivity depending on whether T was above (dry Eq.6) or below (wet Eq. 7) BP. (4) We then calculate k(x, y) everywhere and solved the Helmholtz equation for E(x,y) (Fig. 11), subject to the boundary conditions given by the Cagniard solution. To obtain the fits in Figs. (12, 13) four parameters were varied at each time step using least squares (1) the effective saturation porosity ϕ relative to the reference data (Olhoeft, 1977, Eq. 7) of 14% porosity (2) the vertical, (3) horizontal offsets and (4) regional slope. The latter three parameters take into account background regional effects including sub-water table ones. We concluded that most of the VLF signal is due to the wet hot basalt region either side of the dike. While we do not expect to reproduce the real situation with lateral heterogeneities, and vertical variation of porosity, and water content, this simple

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

=573 oC

1973A T

=137 oC

1973B T

max

73

max

250 200 150 100 50

50 40 30 20 10 20

40

60

80

100

100

o

1975 T

=86.2 C

200

1995 T

max

300

400

500

o

=42.2 C

max

250 200 150 100 50

250 200 150 100 50 100

200

300

400

500

100

=41.3 oC

1997 T

200

2012 T

max

300

400

500

=37.8 oC

max

250 200 150 100 50

250 200 150 100 50 100

200 300 400 Distance (m)

500

100

200 300 400 Distance (m)

500

Fig. 10. Temperatures below the boiling point of water and more than ten degrees above ambient (yellow). Finite difference model for temperatures at the times of the VLF and SP measurements in Hawaii. Note the scale is expanded for the first measurement in 1973A, to show the region where condensed water may accumulate. For both 1973A and 1973B dike temperatures are above BP such that the condensation zone forms a veil above the dike.

E Field 0.4

1000

900

0.35

800

distance (m)

0.3 700 0.25

600

0.2

500

400

0.15

300 0.1 200 0.05

100

100

200

300

400

500

600

700

800

distance (m) Fig. 11. Finite difference model E field (volts/m) for a dipping dike.

900

1000

74

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

1973A 50

40

30

Tilt (%)

20

10

0

−10

−20

−30

−40 −400

−300

−200

−100

0

100

200

300

400

Distance (m) Fig. 12. Fit of FD mode to 1973 15 May data.

1973B 40

20

20

20

0

Tilt (%)

40

0

0

−20

−20

−20

−40 −400−200 0 200 400

−40 −400−200 0 200 400

−40 −400−200 0 200 400

Distance (m) 1995

Distance (m) 1997

Distance (m) 2012

40

40

20

20

20

0

Tilt (%)

40

Tilt (%)

Tilt (%)

1975

40

Tilt (%)

Tilt (%)

1973A

0

0

−20

−20

−20

−40 −400−200 0 200 400

−40 −400−200 0 200 400

−40 −400−200 0 200 400

Distance (m)

Distance (m)

Distance (m)

Fig. 13. VLF tilt measurements above the Hi'iaka dike between 1973 and 2012 and finite difference model fits. Calculated temperatures from a cooling dike were converted to conductivity depending on whether they were above or below the boiling point of water. The VLF equations were solved using finite differences. The dots are digitized data from Zablocki (1976) for 1973A, 1973B and 1995, and our measurements in 1995 1997 and 2012. The lines are the output of the model. The fitting parameters were vertical and horizontal offset, regional slope and effective water saturation. The fitting parameters are given in Table 1.

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

∇:~J cond

Table 1 Inverted parameters of VLF anomalies. Year

ϕ (%)

Hoffset (m)

Voffset (%)

Slope (%/m)

1973A 1973B 1975 1995 1997 2012

34 ± 0.39 15 ± 0.31 11 ± 0.17 2.8 ± 0.17 2.7 ± 0.13 −0.59 ± 0.19

12 ± 0.35 13 ± 0.75 28 ± 0.57 26 ± 3.8 42 ± 3.1 49 ± 22

5.5 ± 0.26 4.4 ± 0.38 4.3 ± 0.22 4.5 ± 0.27 1.5 ± 0.21 6.9 ± 0.29

0.016 ± 0.0021 0.011 ± 0.0033 0.0061 ± 0.002 0.015 ± 0.0017 0.024 ± 0.0015 −0.017 ± 0.0022

model fits the data remarkably well (Figs. 12,13). The values of the fit parameters and their uncertainties are given in Table 1. The values for the saturation porosity depend on how saline the water is, but are within the range seen in water content values from the logs of the summit drill hole in Kilauea (Keller et al., 1979). 5.4. Hydrothermal SP Given the VLF model of the previous section we propose that the SP effect is due to the upward flow of vapor-dominated air near the dike. The hydrothermal calculations indicated it takes the shape of an upward current parallel to the dike of breadth proportional to the thermal regime (e.g., as seen in the analytic calculations of Cheng and Minkowycz, 1977). We thus seek here a simple model to invert the data. We follow the approach of Sill (1983) and show how the source can be approximated by a 2-D electrically polarized trapezoid about the dike with uncompensated positive and negative poles on its upper and lower surfaces, and a width equal to that of the half width of the thermal field. The self-potential, ΔV, that develops due to fluid flow in a porous medium is given by the Helmholtz-Smoluchowski equation. ΔV ¼

εζΔ P μσ

ð8Þ

where ε is the dielectric constant, ζ is the streaming potential, ΔP the pressure difference, μ the dynamic viscosity and σ the conductivity. and the reThe total current is that generated by SP convection, ~J conv

turn current in the medium, ~J cond . The current ~J conv due to SP effects is modeled to be proportional to the gradient of pressure P, which through Darcy's law is proportional to velocity, i.e., ~¼− u

κ∇P μ

Lij are the coupling coefficients(L21 = εζ/μ and L11 = μ/κ), Eq. (9) becomes ~J ¼ − εζ ∇P ¼ − εζ u ~: μ κ

ð11Þ

~ ~J ¼ ~J T cond þ J conv :

ð12Þ

Since ∇:~J T ¼ 0 ∇:~J cond ¼ −∇:~J conv

V¼

I 2πσ r

V SP

1 ¼ 2πσ

ð15Þ   ∇:P~ re0

Z

ð13Þ

j~r −~r0 j

V0

dV 0 :

ð16Þ

For our purposes given we do not know the details of streaming potential, porosity, permeability, pressure, that make up (εζ/κũ) we approximate it with an average value P~ 0 . For example we consider a model of upwards streaming moist air comprised of water evaporated from the water table and vadose zone and lost to the atmosphere. We assume P~ goes from 0 to P~ 0 at the point of evaporation and remains constant and aligned with the dike until the point where moist air escapes to the atmosphere (or condenses) where it drops to zero again. Thus ∇ ~ 0 is current. The integral in Eq. (16) :P~ 0 is non-zero only at the ends. P~ 0 :dS can be converted to a surface integral given by V SP ¼

1 2πσ

Z 0

S

P~ 0 ~ 0 :dS : ~r −~r 0

ð17Þ

In two dimensions (x′, y′), after integrating in the z′ direction, this equation becomes, 1 πσ

Z s0

~ 0 logðj~r −~r 0 jÞP~ 1 :ds

ð18Þ

~ 0 is an elementary lengthwhere now P~ 1 is the polarization per length, ds segment in the x − y plane. Thus the problem becomes one of summing the effects of the uncompensated poles on the upper and lower boundaries. If we consider a sheet of uncompensated poles (as in Supplementary Appendix B Figures B 1,2,3) we can integrate to obtain the potential at a point on the surface (x, y) due to an elementary line of poles at x′ i.e., of current P1 dx′ as, V SP ðx; yÞ ¼ −

With no external sources the total current JT is

ð14Þ

Then the self-potential is given by integrating the current density using the Green's function for a point current in a half space,

V SP ¼ −

ð10Þ

75

∇. (εζ ũ/κ) are the current sources (per unit volume) that give rise to the external potential. In a uniform hydrothermal circulation ∇. ũ = 0. If εζ/κ is constant, hydrothermal circulation in a medium of constant streaming potential generates zero SP voltage variation on the surface (Fitterman, 1978, 1979, 1983; Sill, 1983). Generation of an external field requires the divergence to be non-zero. Thus SP fields rely on the divergence of the pressure gradient, or contrasts in the streaming potential (Fitterman, 1978). One such example is flow across boundaries separating media of contrasting properties. This occurs if any of the ε, ζ, κ, ũ have contrasting properties. Then the divergence becomes non-zero (provided they are not mutually canceling). The physical phenomenon is the electric equivalent of pole cancelation in magnetostatics. Uncompensated poles give rise to external ~. ∇:P~ is the current density. fields. Let the electric polarization, P~ ¼ εζ =κ u

ð9Þ

where κ is permeability. Thus, if as in Sill (1983) we write ~J ¼ −L21 ∇P ¼ −L21 L11 u ~:

  εζ ~ : ¼ ∇:ðL21 ∇P Þ ¼ ∇: u κ

V SP ðx; yÞ ¼ −

1 πσ

ZA P 1 log

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx−x0 Þ2 þ y2 dx0

ð19Þ

−A

P1 ½ logðr 1 r 2 Þ þ x logðr 2 =r 1 Þ−2A þ yθ  πσ

ð20Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x−A where θ ¼ atanðxþA r 1 ¼ ððx−AÞ2 þ y2 Þ , r2 y Þ−atanð y Þ , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ¼ ððx þ AÞ þ y2 Þ and y is the depth to the surface. For the lower

76

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

boundary y = W the depth of the water table and for a tilted trapezoid at angle ϕ the offset will be x0 = (W − y)/tan(ϕ). After adding the contribution from the water table at the bottom, we obtain for the potential at a point x on the surface a distance y above the upper surface of the trapezoid, V SP ðx; yÞ ¼ −

      P1 r1 r2 r2 r3 þ x log þ x0 logðr4 =r3 Þ þ yθ1 −Wθ2 : A log πσ r3 r4 r1 r4

ð21Þ Given that the dipolar sources for the SP variation are located in the vadose zone and that the signal has persisted over 39 years, a period when the dike would be expected to have cooled below the boiling point of water, we suggest that the SP signal is due to the contrast between evaporation of meteoric water in the deep vadose zone and the influx of meteoric water. Recalling that significant SP ζ potentials have been observed associated with evaporation of rain water (Thony et al., 1997; Darnet and Marquis, 2004), a dike several tens of °C above ambient will cause differential evaporation in its vicinity. Contrasts in the flow occur where the water evaporates and when the moisture-laden water either condenses or reaches the low porosity zone near the surface and ultimately the surface itself. Down-flowing meteoric water will be diverted by the up flow but as distance increases it will give rise to a negative potential e.g., as observed in the drill hole of 1.35 mv/m (Keller et al., 1979). 5.5. SP inversion For a start we inverted the SP data taken each year using a simple di~ made up of pospole model. As in magnetism an SP dipole of strength D ~ itive and negative currents separated by a distance dr has a voltage given by V ¼−

1 e0 1 D :∇ 2πσ r

ð22Þ

~ ¼ P 1  area, in two diLet the dipole strength per unit distance be D mension Eq. (22) becomes V¼

1 ~ D:∇logðr Þ πσ

ð23Þ

We invert the data using Eq.(23) including a regional slope and offsets (Table 2). The fits are satisfactory and show that the depth of the sources increase monotonically with time in the range 50–190 m, that is, in the vadose zone. However in all such inversions there is a tradeoff between depth and finite size of the body (Parker, 1975). Hydrothermal modeling indicated the upflow extends either side of the dike over a distance equal to the half-width of the temperature field. We then approximated a finite source using a model comprised of a 2-D dipping, trapezoidal, cross-section of constant electric polarization parallel to the long axis in the vertical plane. The voltage then is due to uncompensated poles on its top and its bottom. We fixed the bottom at 230 m, the inferred depth of the water table. The tilt of the trapezoid was fixed at the value of the tilt of the dike obtained from the VLF inversion. The cross-sectional widths were fixed equal to the half-widths of the

temperature field calculated for each survey time (Fig. 6). The sides are taken to be parallel to the dike. We inverted for the five parameters: polarization strength, the depth of the upper surface, the vertical and horizontal offsets and the regional slope. Fig. 14 shows that the trapezoidal model fits the SP data well, but with so many parameters this may not be too surprising. Fig. 15 shows the inverted models that produced these fits. Table 3 lists the assumed widths from the thermal modeling, and inverted parameter values and their standard deviations. All of the parameters are well constrained. Note the dipole strength is calculated as P1 × area to compare with dipole strength in Table 2. To convert to volts/m we divide by the area to obtain 978, 20, 21, and 78 mv/m for the surveys in 1973, 1995, 1997, and 2012, respectively. The first of these values is high, but that may be the result of too small a width (2 m) of the upflow zone at a time of intense degassing of the dike. Values of 25 mv/m have been measured for evaporating rainfall in soils (Thony et al., 1997). Thus the remaining values are plausibly generated by evaporation caused by heating by a dike. This EK field is superposed on the background field of negative SP (1.35 mv/m) from influx of meteoric water. The SP source widens, and deepens with time, as might be expected. We also performed inversions with the depth of the water table set to 150 m, that is, the value that fit the VLF data. The fits are similar. The data alone cannot distinguish between the two cases. We note in 2012 the deep thermal anomaly may explain why there was no VLF signal at that time. At a depth of N 150 m, and a smaller temperature contrast, the anomaly is both too weak and situated too low to generate a measurable VLF field. 6. Discussion and conclusions Self-Potential and VLF measurements over 39 years above an intruded dike in Kilauea's rift zone show distinctly different behavior. The VLF signal decayed rapidly while the SP remained at 60% of it original value. The more rapid decay of the VLF signal compared with that of the SP anomaly is explained by the generation and decay of a hot saturated vadose zone about the dike having higher conductivity than the surroundings. While, in contrast, the SP anomaly is due to evaporation that depends on stored heat which has spread out from the dike over some 100 m owing to the slow rate of diffusion with only a small amount lost from the surface. Although effects from below the water table may occur, all the VLF and SP signals could be modeled by temperature and water circulation in the vadose zone. Hydrothermal circulation beneath the water table is not expected to generate a large SP anomaly on the surface, unless there are marked contrasts in properties. VLF anomalies are not expected from there either, because the inducing field has decayed to near zero owing to high conductivity. However the long-wavelength regional trends required to fit the data, though smaller than vadose zone effects, may originate from below the water table. The thermal evolution about the intruded dike was modeled based on diffusivity values that fit the thickening crust of Kilauea Iki lava lake. Hydrothermal modeling suggested that in the vadose zone the dominant redistribution of heat is by diffusion rather than advection. The diffusion anomaly from a dipping dike 25 m below the surface shows a progressive expansion centered about the dike widening to 84 m after 39 years.

Table 2 Inverted parameters of self-potential anomalies assuming a dipole source. Year

D (Vm)

Depth (m)

Hoff (m)

Voff (mv)

Slope (mv/m)

1973 1995 1997 2012

290 ± 17.6 213 ± 15.6 243 ± 9.21 416 ± 43.9

84 ± 4.4 80 ± 5.6 1.2e + 02 ± 3.6 2e + 02 ± 14

−7.9 ± 1.6 −38 ± 2.6 −10 ± 1.5 −46 ± 5.2

−0.11 ± 0.024 0.042 ± 0.019 0.062 ± 0.0087 −0.0083 ± 0.029

−0.76 ± 0.069 −0.31 ± 0.049 0.28 ± 0.018 0.28 ± 0.034

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

Self Potential 1973

Self Potential 1995 1

SP (mVolts)

SP (mVolts)

1

0.5

0 −500

0

500

−500

0

500

Distance (m)

Distance (m)

Self Potential 1997

Self Potential 2012 1

SP (mVolts)

SP (mVolts)

0.5

0

1

0.5

0 −500

77

0.5

0 0

500

−500

Distance (m)

0

500

Distance (m)

Fig. 14. Self-potential measurements above the Hi'iaka dike between 1973 and 2012 showing the effect of cooling. The fit of a simple trapezoidal model is shown by the lines. The parameters are listed in Table 2.

A finite difference VLF model was constructed in which the dike intrudes a two-layered medium consisting of water table and vadose zone. For the first few years the VLF and SP signals remained at constant half-width. This is expected from the thermal modeling, since the

diffusion of heat is so slow. However after 22 and 39 years appreciable widening is seen, consistent with a diffusing thermal field. A thermal diffusion model of a dipping dike in which temperature above and below the boiling point of water were converted to conductivity based

Dike Model 1995 0

−50

−50

Depth (m)

Depth (m)

Dike Model 1973 0

−100 −150 −200 −500

−100 −150 −200

0

500

−500

Distance (m)

500

Dike Model 2012 0

−50

−50

Depth (m)

Depth (m)

Dike Model 1997 0

−100 −150 −200 −500

0

Distance (m)

−100 −150 −200

0

Distance (m)

500

−500

0

500

Distance (m)

Fig. 15. Inverted models for region developing SP polarization. We propose that the widths are due to thermal diffusion as the dike cools and define the region where incipient heat evaporates ground water. The rising ground water develops SP polarization. The sources are assumed to develop between the water table (230 m) and the top of the still-hot dike increases in depth with time as it loses heat to the surface. The absence of a VLF signal in 2012 is explained as the cooled dike is deeper than the zone of VLF sensitivity.

78

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

Table 3 Inverted parameters of trapezoidal self-potential source assuming base at a depth of 230 m. Year

W

D (Vm)

Depth (m)

Hoff (m)

Voff (v)

slope (mv/m)

1973 1995 1997 2012

2 62 66 84

385 ± 13 264 ± 7.5 257 ± 6.3 396 ± 41

33 ± 3.2 17 ± 2.6 48 ± 3.3 1.7e + 02 ± 25

−57 ± 1.6 −94 ± 1.5 −54 ± 1.6 −57 ± 6.9

−0.22 ± 0.02 −0.02 ± 0.01 0.047 ± 0.007 −0.0072 ± 0.03

−0.7 ± 0.06 −0.23 ± 0.03 0.32 ± 0.02 0.28 ± 0.03

on laboratory results, provides an excellent fit to the VLF data. The parameter, ϕ describing effective water saturation decays monotonically over time (Table 1). The effective depth of the water table was found to be 150 m, rather than the 230 m inferred from VES in the vicinity. We are unsure if the difference is due to inadequacy of our layered VLF model or if it is real. Thus our final model (Fig. 16) is an initially hot dike (1200 °C) 1 m wide, dipping at 72°, which intruded the crust for a distance of 2 km with a lower depth of over 1 km and upper depth of 25 m. After 10 days it had cooled to (500 °C) when the first measurements were made. The dike would then be degassing steam and boiling the water table and vadose zone, but at a distance away from it, where temperatures were below the boiling point, condensation would occur. There the vadose zone would become saturated with hot water, presenting a conductivity anomaly. Closer to the dike the unsaturated region of upwelling steam generates the self-potential.

(a)

As the temperature decreases below the boiling point, the recharge of the saturated vadose zone ceases and it drains away. But evaporation continues, albeit spread over a larger region. By 2012, the temperature had dropped to 38 °C (100 °F), still capable of significant evaporation. The inferred streaming potentials of tens of mv/m are consistent with values seen in evaporation of rainfall in soils. Rapid fluid disruption of water molecules (Johnston et al., 2001) provides the best explanation of the SP anomalies in the early stages when temperatures are above boiling point. Boiling of the water table and fluid in the vadose zone would form an upward flow of steam next to the dike, which either escapes to the atmosphere, or condenses away from it to join the downward flow of meteoric water. Thus the electric polarization of the upwards flow is truncated at each end, leaving uncompensated (current source) poles separated by a low conductivity medium, conditions necessary for generation of SP on the surface.

Dike water table

Sea level

Brackish water Sea floor Salt water

(b)

+

Rainfall SP

Condensed Water

+

upflow moist air or steam Water table

-

-

+ Evaporation

Fig. 16. Model of generation of SP anomaly in the Vadose zone by the contrast between downward flow of meteoric water to the water table and upward flow of evaporation caused by the thermal anomaly of a cooling dike. (a) Schematically shows the water table beneath the island and how condensation can effectively uplift it as a veil over the dike. (b) Shows an expanded version of this effect where the inner zone of low-conductivity rising moist air generates the SP signals.

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

As temperatures cool well below boiling, such as in the mid nineties or 2012, evaporation, fluid disruption, and up-flow, will still occur next to a hot dike, but distributed over a much wider area, as the thermal field has diffused outwards. Thus the continued strength of the SP field and its broadening with time are likely the result of fluid disruption as the surrounding vadose zone and underlying water table is dried by evaporation. The time-varying regional field (Tables 1–3), that was fit with a sloping line, may reflect longer wavelength effects generated beneath the water table, but they are small compared with the vadose zone affects modeled here. We conclude that the combination of SP and VLF measurements can characterize the thermal regime of dikelike intrusions above the water table that could quantify the associated geothermal resource. Acknowledgments Staff at the Hawaiian Volcano Observatory who hosted the field trips are thanked, especially James Kauahikaua who hosted our visits, and supported the experiments and analysis. Members of field classes from UCLA, Professors Robert McPherron and Paul Tackley, and Cecily Davis are thanked for assistance in the field work. Comments by an anonymous reviewer helped improve the manuscript. This work was funded in part by the Office of Instructional Development, UCLA. The U.S. National Park Service is thanked for permission to carry out the study. Supplementary Data. Appendices A and B Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.jvolgeores.2015.06.007. References Anderson, L.A., Johnson, G.R., 1976. Application of self-potential method to geothermal exploration in Long Valley, California. J. Geophys. Res. 81 (8), 1527–1532. Antraygues, P., Aubert, M., 1993. Self-Potential generated by two-phase flow in a porous medium: experimental study and volcanological applications. J. Geophys. Res. 98 (B12), 22,273–22,281. Banwell, G.J., 1970. Geophysical techniques in geothermal exploration. Proceedings of the U.N. Symposium of the Development and Utilization of Geothermal Resources. vol. 1, p. 32. Barth, G.A., Kleinrock, M.C., Helz, R.T., 1994. The magma body at Kilauea Iki lava lake: potential insights into mid-ocean ridge magma chambers. J. Geophys. Res. 99 (B4), 7199–7217. Cagniard, L., 1953. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics 18 (3), 605–635. http://dx.doi.org/10.1190/1.1437915. Carslaw, H.S., Jaeger, J.C., 1986. Heat Conduction in Solids. Oxford University Press (520 pp.). Cervelli, P., Segall, P., Amelung, F., Garbeil, H., Meertens, C., Owen, S., Miklius, A., Lisowski, M., 2002. The 12 September 1999 upper East Rift Zone dike intrusion at Kilauea volcano, Hawaii. J. Geophys. Res. 107 (B7). http://dx.doi.org/10.1029/2001JB000602. Cheng, P., Minkowycz, W.J., 1977. Free convection about a vertical flat plate embedded in a porous medium with application to heat transfer from a dike. J. Geophys. Res. 82 (14), 2040–2044. Combs, J., Wilt, M.J., 1976. Telluric mapping,telluric profiling, and self-potential surveys of the Dunes geothermal anomaly, Imperial Valley,California. Proc. 2nd U.N. Symposium on the Development and Use of Geothermal Resources. vol. 2, pp. 937–945. Corwin, R.F., 1976. Self-Potential exploration for geothermal reservoirs. Proceedings of the Second U.N. Symposium on the Development and Use of Geothermal Resources. vol. 2, pp. 937–945. Corwin, R.F., Hoover, D.B., 1979. The self-potential method in geothermal exploration. Geophysics 44, 222–226. Darnet, M., Marquis, G., 2004. Modelling streaming potential (SP) signals induced by water movement in the vadose zone. J. Hydrol. 285, 114–124. Davis, P.M., 2014. Electromagnetic induction in a conductive strip in a medium of contrasting conductivity: application to VLF and MT above molten dykes. Geophys. J. Int. 199 (2), 1176–1186. Dvorak, J.J., Okamura, A.T., English, T.T., Koyanagi, R.Y., Nakata, J.S., Sako, M.K., Tanigawa, W.T., Yamashita, K.M., 1986. Mechanical response of the south flank of Kilauea Volcano, Hawaii, to intrusive events along the rift systems. Tectonophysics 124, 193–209. Fitterman, D.V., 1978. Electrokinetic and magnetic anomalies associated with dilatant regions in a layered earth. J. Geophys. Res. 83, 5923–5928. Fitterman, D.V., 1983. Modeling of self-potential anomalies near vertical dykes. Geophysics 4 (8), 171–180. Fitterman, D.V., Corwin, R.F., 1982. Inversion of self-potential data from the Cerro Prieto geothermal field. Mexico. Geophysics 47, 938–945.

79

Fittermann, D.V., 1979. Theory of the electrokinetic effect in a faulted half space. J. Geophys. Res. 84 (6-31-6041). Furumoto, A.S., 1975. A systematic program for geothermal exploration onthe island of Hawaii. paper presented at the 45th Annual International Meeting. Soc. of Explor. Geophys, Denver, Colo (Oct. 12–16). Hashimoto, Y., Tanaka, Y., 1995. A large self-potential anomaly on Unzen volcano, Shimbara peninsula, Kyushu island, Japan. Geophys. Res. Lett. 22 (3), 191–194. Hohmann, G.W., 1971. Electromagnetic scattering by conductors in the Earth near a line source of current. Geophysics 36 (1), 101–131. Ishido, T., 1981. Streaming potential associated with hydrothermal convection in the crust: a possible mechanism of self-potential anomalies in geothermal areas. J. Geophys. Res. Soc. Jpn. 3, 87–100 (in Japanese with English abstr.). Ishido, T., 1989. Self-potential generation by subsurface water flow through electrokinetic coupling. In: Merkler, G.P., et al. (Eds.), Detection of Subsurface Flow Phenomena, Lecture Notes in Earth Sciences vol. 27. Springer-Verlag, pp. 121–131. Ishido, T., Mizutani, H., 1981. Experimental and theoretical basis of electrokinetic phenomena in rock–water systems and its applications to geophysics. J. Geophys. Res. vol. 86 (B3), 1763–1775. Ishido, T., Pritchett, J., 1996. Numerical simulation of electrokinetic potentials associated with subsurface fluid flow. Proc. Twenty-First Workshop on Geothermal Reservoir Engineering. Stanford University, Stanford, California, pp. 143–149 (Jan. 22-24). Ishido, T., Pritchett, J., 1999. Numerical simulations of electrokinetic potentials associated with subsurface fluid flow. J. Geophys. Res. 104, 15,247–15,259. Ishido, T., Pritchett, J., 2000. Using numerical simulation of electrokinetic potentials in geothermal reservoir management. Proceedings World Geothermal Congress 2000 Kyushu — Tohoku, Japan, pp. 2629–2634 (May 28–June 10). Ishido, T., Mizutani, H., Baba, K., 1983. Streaming potential observations, using geothennal wells and in situ electrokinetic coupling coefficients under high temperature. Tectonophysics 91, 89–104. Jackson, D.B., Kauahikaua, J.P., 1987. Regional self-potential anomalies at Kilauea volcano. U. S. Geol. Surv. Prof. Pap. 1350 (40), 947–959. Jackson, D.B., Keller, G.V., 1972. An electromagnetic sounding survey of the summit of Kilauea volcano, Hawaii. J. Geophys. Res. 77 (26), 49574965. Jackson, D.B., Sako, M.K., 1982. Self-potential surveys related to probable geothermal anomalies, Hualalai volcano, Hawaii. U.S. Geol. Surv. Open File Rep. 82–127. Jackson, D.B., Kauahikaua, J.P., Zablocki, C.J., 1985. Resistivity monitoring of an active volcano using the controlled-source electromagnetic technique: Kilauea, Hawai'i. J. Geophys. Res. 90, 12,545–12,555. Jaeger, J.C., 1957. The temperature in the neighborhood of a cooling intrusive sheet, Amer. J. Mar. Sci. 255, 306–318. Johnston, M.J.S., 1997. Review of electric and magnetic fields accompanying seismic and volcanic activity. Surv. Geophys. 18, 441–475. Johnston, M.J.S., Byerelee, J.D., Lockner, D., 2001. Rapid fluid disruption: a source for selfpotential anomalies on volcanoes. J. Geophys. Res. 106 (3), 4327–4335. Kauahikaua, J.P., 1993. Geophysical characteristics of the hydrothermal systems of Kilauea volcano, Hawai'i. Geothermics 22 (4), 271–299. Kauahikaua, J.P., Mattice, M.D., 1981. Geophysical reconnaissance of prospective geothermal areas on the island of Hawai'i using electrical methods. U.S. Geol. Surv. Open File Rep. 81–1044 (50 pp.). Kauahikaua, J.P., Mattice, M.D., Jackson, D.B., 1980. Misc-a-la-Masse mapping of the HGP—a geothermal reservoir. Geotherm. Res. Count. Trans. 4, 65–68. Kauahikaua, J.P., Jackson, D.B., Zablocki, C.J., 1986. Resistivity structure to a depth of 5 km beneath Kilauea volcano, Hawai'i from large-loop-source electromagnetic measurements (0.04–8 Hz). J. Geophys. Res. 91, 8267–8283. Kauahikaua, J.P., Cashman, K.V., Mattox, T.N., Heliker, C.C., Hon, K.A., Mangan, M.T., Thornber, C.R., 1998. Observations on basaltic lava streams in tubes from Kilauea Volcano, island of Hawai'i. J. Geophys. Res. 103 (B11), 27303–27323. Keller, G.V., Trowbridge Grose, L., Murray, J.C., Skokan, C.K., 1979. Results of an experimental drill hole at the summit of Kilauea volcano, Hawaii. J. Volcanol. Geotherm. Res. 5, 345–385. Klein, F.W., Koyanagi, R.Y., Nakata, J.S., Tanigawa, W.R., 1987. The seismicity of Kilauea's magma system. In: Decker, R.W., Wright, T.L., Stauffer, P.H. (Eds.), Volcanism in Hawaii: U.S. Geological Survey Professional Paper. 1350, pp. 1019–1185. Marsden Jr., S.S., Tyran, C.K., 1986. The potential generated by flow of wet steam in capillary tubes. Proceedings, Eleventh Workshop on Geothermal Reservoir Engineering Stanford University, Stanford, California, January. 21–.23, 1986 SGP-TR-93, pp. 225–229. Morgan, F.D., Williams, E.R., Madden, T.R., 1989. Streaming potential properties of Westerly granite with applications. J. Geophys. Res. 94, 12,449–12,461. Nield, D.A., Bejan, A., 1992. Convection in Porous Media. Springer-Verlag, NY (408 pp.). Nourbehecht, B., 1963. Irreversible Thermodynamic Effects in Inhomogeneous Media and their Applications in Certain Geoelectric Problems, Ph.D. thesis, M. I.T. Olhoeft, G.R., 1977. Electrical properties of water saturated basalt, Preliminary results to 506 K (233 °C). U.S. Geol. Surv. Open File Rep. D 77–688. Olhoeft, G.R., 1981. Electrical properties of granite with implications for the lower crust. J. Geophys. Res. 86, 931–936. Parker, R.L., 1975. The theory of ideal bodies for gravity interpretation. Geophys. J. R. Astron. Soc. 42, 315–334. Presnall, D.C., Simmons, C.L., Porath, H., 1972. Changes in electrical conductivity of a synthetic basalt during melting. J. Geophys. Res. 77 (29), 5665–5672. Revil, A., Schwaeger, H., Cathles, L.M., Manhardt, P., 1999. Streaming potentials in porous media, 2, Theory and application to geothermal systems. J. Geophys. Res. 104, 20,033–20,048. Richter, D.H., Moore, J.G., 1966. Petrology of the Kilauea Iki lava lake Hawaii. U. S. Geol. Surv. Prof. Pap. 537B, B1–B26. Sasai, Y., Zlotnicki, J., Nishida, Y., Yvetot, P., Morat, P., Murakami, H., Tanaka, Y., Ishikawa, Y., Koyama, S., Sekiguchi, W., 1997. Electromagnetic monitoring of

80

P.M. Davis / Journal of Volcanology and Geothermal Research 302 (2015) 64–80

Miyake-jima Volcano, Izu-Bonn Arc, Japan: a preliminary report. J. Geomagn. Geoelectr. 49, 1293–1316. Scarlato, P., Poe, B.T., Freda, C., 2004. High-pressure and high-temperature measurements of electrical conductivity in basaltic rocks from Mount Etna, Sicily, Italy. J. Geophys. Res. 109 (B02210), 1–11. http://dx.doi.org/10.1029/2003JB002666. Schloessin, H.H., 1976. A comparison of the effects of pressure temperature and H2O adsorption and desorption on the electrical conductivity of basalts. J. Geophys. Res. 81 (23), 4077–4086. Sill, W.R., 1983. Self-potential modeling from primary flows. Geophysics 48, 76–86. Spengler, S.R., Walker, G.P.L., 2012. Dike Width Characteristics Observed in Intra-plate Volcanoes, Element Environmental, Hawaii (Website, http://www.e2hi.com/wpcontent/uploads/2012/10/11–Dike-Width-Paper.pdf). Thony, J.L., Morat, P., Vachaud, G., Le Moue¨l, J.L., 1997. Field characterization of the relationship between electrical potential gradients and soil water flux. CR Acad. Sci. Paris 325, 317â€'321. Tilling, R.I., Christiansen, R.L., Duffield, W.A., Endo, E.T., Holcomb, R.T., Koyanagi, R.Y., Peterson, D.W., Unger, J.D., 1987. The 1972–1974 Mauna Ulu eruption, Kilauea Volcano: an example of quasi-steady-state magma transfer. In: Decker, R.W., Wright, T.L., Stauffer, P.H. (Eds.), Volcanism in Hawaii: U.S. Geological Survey Professional Paper. 1350, pp. 405–469. Turcotte, D.L., Schubert, G., 2002. Geodynamics. University press, Cambridge (466 pp.). Walker, G.P.L., 1986. Koolau dike complex, Oahu: intensity and origin of a sheeted-dike complex high in a Hawaiian volcanic edifice. Geology 14, 310–313.

Walker, G.P.L., 1987. The dike complex of Koolau volcano, Oahu: internal structure of a Hawaiian rift zone. Volcanism in Hawaii. U. S. Geol. Surv. Prof. Pap. 1350, 961–993. Walker, G.P.L., 1988. Three Hawaiian calderas: an origin through loading by shallow intrusions? J. Geophys. Res. 93 (B12), 14773–14784. Yang, X., Davis, P.M., Delaney, P.T., Okamura, A.T., 1992. Geodetic analysis of dike intrusion and motion of the magma reservoir beneath the summit of Kilauea volcano, Hawaii: 1970–1985. J. Geophys. Res. 97 (B3), 3305–3324. Yasukawa, K., Bodvarsson, G.S., Wilt, M.J., 1993. A coupled self-potential and mass-heat flow code for geothermal applications. GRC Transactions. vol. l7, pp. 203–207. Zablocki, C.J., 1976. Mapping thermal anomalies on an active volcano by the self-potential method, Kilauea, Hawaii. Proc., 2nd U.N. Symposium on the Development and Use of Geothermal Resources, San Francisco, Calif., May 1975. 2, pp. 1299–1309. Zablocki, C.J., 1978a. Applications of the VLF induction method for studying some volcanic processes of Kilauea volcano, Hawaii. J. Volcanol. Geotherm. Res. 3, 155–195. Zablocki, C.J., 1978b. Streaming potentials resulting from the descent of meteoric water—a possible source mechanism for Kilauean self-potential anomalies. Geotherm. Res. Counc. Trans. 747–748 (Vol., 2 July 1978). Zlotnicki, J., Le Mouel, J.-L., 1990. Possible electrokinetic origin of large magnetic variations at La Fournaise, Volcano. Nature 343, 633–636. Zohdy, A.A.R., Anderson, L.A., Muffler, L.J.P., 1973. Resistivity, self-potential, and induced polarization surveys of a vapor-dominated geothermal system. Geophys. J. R. Astron. Soc. 38, 1130–1144.