Journal of Volcanology and Geothermal Research 92 Ž1999. 413–429 www.elsevier.comrlocaterjvolgeores

Practical evaluation of steady heat discharge from dormant active volcanoes: case study of Vulcarolo fissure žMount Etna, Italy / Maurice Aubert

)

UniÕersite´ Blaise Pascal et CNRS, Laboratoire de Geologie, 5 rue Kessler, 63038 Clermont cedex, France ´ Received 28 October 1998; accepted 12 May 1999

Abstract From a theoretical approach combined with experimental study, I propose a method to compute both sensible and latent heat transfer in a thermally active volcanic area. This method is applied to a fumarolic thermal fissure located on the south flank on Mount Etna ŽItaly.. The method involves four steps: Ž1. mapping of Self-Potential ŽSP. to quickly identify thermal anomalies; Ž2. mapping of ground temperature ŽT . at constant depth in the area defined by SP mapping; Ž3. measurement of vertical T profiles at selected representative stations including the coldest and hottest; and Ž4. direct measurements of vapor flow where the convective heat flux reaches the surface. When heat flux values vary in the range 10–300 W my2 , T measurements at a depth equal to about 1 m are sufficient to estimate the absolute flux value to within an absolute error less than 30%, and 10% for the relative value. Ground temperature measurement is an easy and accurate method for estimating heat discharge on a limited surface to define the fissure system, monitor time variations of heat flux, or to calibrate remote sensing results on the ground. q 1999 Elsevier Science B.V. All rights reserved. Keywords: heat transfer; heat flux; Mount Etna

1. Introduction One of the more obvious signals of thermal activity at a dormant volcano is heat and gas emission across fumarolic zones ŽWilliams and McBirney, 1979.. Thermal convection above magma bodies is treated by Hardee Ž1982. as a one-dimensional bottom-heated convection process in a permeable medium. Fig. 1 shows the different zones above magma bodies, from the molten magma up to the

)

Tel.: q33-4-73-34-67-49; fax: q33-4-73-34-67-44; E-mail: [email protected]

surface, through a heat exchanger. Because of the assumed presence of this exchanger, variation of flux at the surface may not involve variation of magmaderived flux: the exchanger may absorb or return a part of heat received from the magma. Near the surface, heat flows upward by conduction and twophase permeable convection, the condensation of water being completed in the atmosphere Žfumarolic zone sensu-stricto. or in the soil Žsub-fumarolic zone.. The latent heat of cooling water creates a ground temperature ŽT . anomaly in the sub-surface which is measurable ŽAubert et al., 1984.. In most cases, the top of the convective cell does not reach the surface. This causes sub-fumarolic zones as defined by Aubert et al. Ž1984.. Such sub-fumarolic zones can be de-

0377-0273r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 7 3 Ž 9 9 . 0 0 0 8 8 - 8

414

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

Fig. 1. Thermal zones above magma bodies. ŽA. molten magma body; ŽB. crystalline and progressively crystallising intrusion; F 1 heat flux from magma body in transition zone where conductive flux dominates; ŽC. heat exchanger where two-phase permeable convection flux dominates; F 2 heat flux toward atmosphere. Adapted from Hardee Ž1982., and Harris and Stevenson Ž1997a,b..

tected either by a relative weak T anomaly when the top of the convective cell is near the surface or, if the cell is deeper, they can be detected by geophysical methods such as SP measurements or by gas analysis ŽAubert and Baubron, 1988.. Thermal mapping and monitoring of thermal anomalies at the surface Ždue to active lava bodies such as flows and domes. can be carried out by means of satellite thermal infrared ŽIR. data ŽYuhara et al., 1981; Oppenheimer et al., 1993; Gaonac’h et al., 1994; Harris and Stevenson, 1997a,b; Wooster and Rothery, 1997.. Sekioka and Yuhara Ž1974. proposed a heat balance method to estimate the heat discharge using surface temperatures and some meteorological parameters. Kagiyama Ž1981. uses this method to calculate the thermal energy released by non-eruptive mechanisms associated with active volcanoes in Japan. From this study, it appears that non-eruptive heat discharge is of the same order as that caused by eruptions. From continuous ground T measurements along vertical profiles, Tabbagh and Trezeguet Ž1987. proposed a method for calculation of the heat transfer by both conduction and convection, distinguishing the steady heat flux from the unsteady flux. But heat

flux in areas of predominant convective transfer essentially depends on the liquid vapor content of the convective fluid. These authors do not give flux values, but only propose two cases of liquid vapor content, for which the flux value varies over more than two orders of magnitude. To my knowledge, these remotely sensed thermal data have not been compared with accurate values of thermal flux directly measured from ground temperatures. Therefore I propose a practical evaluation of steady heat discharge by both conduction and convection from a simplified convection model, and using only ground T measurements. I demonstrate the validity of this model, using a well-documented study of a sub-fumarolic zone of Mount Etna volcano ŽSicily.. Values are compared with those obtained by field data including liquid vapor content of the fluid.

2. Theoretical approach Hardee Ž1982. described a complete cycle of permeable convection, from magma bodies to the surface. Here, I only consider hydrothermal convection in the last meter below the surface where ground temperature T can normaly be measured. I first make some simplifying assumptions. These are examined further when comparing the field-collected data and the theoretical model. The assumptions are: ŽA. the medium is porous, with homogeneous physical and thermal parameters; ŽB. at each station, convection corresponds to a one-dimensional distribution ŽT and Q, the heat quantity produced by volume unit, are functions of the vertical axis Z alone, positive downward.; ŽC. the volcano is in thermal steady-state; ŽD. the vertical convective flow Dz goes to zero at the surface. The equation to be solved is: d2 TrdZ 2 q Qrl q Ž dTrdZ . C v r Ž Vrl . s 0

Ž 1.

where T is the ground temperature, Q the heat quantity produced per unit volume in the medium, l the thermal conductivity, r the density of vapor, C v

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

the specific heat of vapor, V the vapor speed. Derivatives are total and not partial in relation to assumption ŽB.. The first term of Eq. Ž1. corresponds to conductive transfer, the second term to heat production in the medium, and the third term to convective heat transfer. The third term is the most difficult to estimate, the main difficulty being the uncertainty in the V and C v values, the latter parameter depending on the liquid vapor state in the fluid. C v values can range from 0.8 = 10 3 J my3 Ky1 for dry gas to 0.4 = 10 6 J my3 Ky1 for moist gas with 10% water content ŽTabbagh and Trezeguet, 1987.. To obtain V, it is necessary to: Ži. measure the vapor flow, Žii. estimate the effective porosity, and Žiii. estimate the percentage of non-condensable gas in the flow Žessentially atmospheric gases.. In fact, when the assumption ŽD. is fully satisfied, and when the object is to estimate the heat transfer from the ground to the atmosphere, one may consider that convective heat creates heat in the medium, and is thus entirely included in Q. Q includes three components: Ž1. heat released by condensation Ž2350 = 10 3 J kgy1 .; Ž2. heat release by cooling of the ascendant vapor flow Ž1500 J kgy1 Ky1 .; and Ž3. heat expended by heating of the descending condensed vapor Ž4000 J kgy1 Ky1 .. The first term dominates, especially because the two others terms partially counterbalance each other. Eq. Ž1. can thus be simplified to: d 2 TrdZ 2 q Qrl s 0

Ž 2.

Fig. 2 shows the typical aspect of a T Ž Z . profile in the first few meters below the surface. This profile shows three different parts, according to the value of Q. Ž1. From the deepest zone to depth Z2 where Q ( 0, the water vapor flow corresponds to the ascent of vapor created by vaporization of meteoric water in deep layers, where the temperature reaches and exceeds the boiling temperature of water. By means of fissures, which create high permeability in the medium, the gas phase Žair, water-vapor, trace gases. is able to rise, the rate of rise depending on the pressure difference of saturated vapor Ž F . along the gas column. Within the vapor flow Ž F ., the partial pressure pe is very high; if pe ) F, then

415

Fig. 2. Variation of the ground temperature T as a function of depth Z, at a station where the heat transfer is both by conduction and convection Žarbitrary units.. The sensible heat transfer is predominant and the dT r dZ ratio constant from the surface to Z1 as well as AT depths greater than Z2 . On the contrary, the latent heat transfer is dominant from Z1 to Z2 .

partial condensation occurs. The condensation rate depends on Ž pe yF., and this rate has an effect on the upstream flow. The heat release Ž Q . balances losses of heat toward the surface and possibly toward the fissure walls by conduction. Ž2. From Z2 to Z1 where Q ) 0, T and vapor flow is much reduced by condensation which is assumed to be nearly complete at Z1 Žflow rate is equal to zero at Z1 .. The heat released by condensation balances the conductive heat flux towards the cold atmosphere. Condensation creates an increase of water content in the ground. Ž3. From Z1 to the surface, Q is equal to zero and the heat transfer is solely conductive. The heat discharge towards the atmosphere is given by the equation:

F s ldTrdZ

Ž 3.

Computation of Q can be made from the T profile ŽFig. 2., using a Taylor approximation formula, whereby Q zrl at depth Z is given by: Q zrl s y Ž TZqd Z q TZyd Z y 2TZ . rd Z 2 where d Z is the differential of Z.

Ž 4.

416

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

Fig. 3. Three aspects of vertical temperature profiles on fumarolic and sub-fumarolic zones. Ža. that is the same case as in Fig. 2, i.e., the water steam flow condensation begins at depth Z2 and is entirely ended at depth Z1 beneath the surface Žsub-fumarolic zone.; Žb. level Z1 is just at the surface; Žc. condensation ended in atmosphere Žfumarolic zone sensu-stricto.. Arrows indicate steam flux.

In the case where Q is assumed constant from Z1 to Z2 , and dTrdZ s 0 from Z2 towards greater depth, integration of Eq. Ž1. yields: Qrl s 2 Ž T2 y T1 . r Ž Z2 y Z1 .

2

Ž 5.

in which T1 and T2 are the temperatures at depths Z1 and Z2 , respectively. In field examples, three cases may also be defined, depending on the position of Z1 with regard to the surface ŽFig. 3.. Ž1. The first case Ža. occurs when Z1 is below the surface and the three sections in the profile of Fig. 2 are well identified. Heat transfer to the atmosphere is calculated from Eq. Ž3.. This heat is also equal to the sum of Q values from Z2 to Z1. Ž2. The second case Žb. occurs when Z1 is at the surface, but condition ŽD. is still satisfied. The heat transfer may be estimated from the sum Q. Ž3. The final case Žc. occurs when Z1 is above the surface and then condition ŽD. is not satisfied. Q is produced partly in the ground, partly in the atmosphere Žfumarolic field strictly speaking.. Only the heat produced in the ground can be measured from ground-temperature.

defined by Kieffer Ž1983., at an average altitude 2900 m. Discrete fissures cut this part of the volcano, and evidence of convective phenomena exists at Vulcarolo, 600 m from the volcano center. At this site ŽVulcarolo, Fig. 4., a condenser was used to provide water to the old observatory destroyed by eruptive activity in 1971. A volume of 300 l dayy1 was condensed, which theoretically converts to a

3. Measurements 3.1. Location of the study area The study area ŽFig. 4. is located at the base of a scoria cone on the south Etna rift zone ŽSicily., as

Fig. 4. Location of the Vulcarolo fissure ŽMount Etna, Italy., between the summit cone and Cisternazza Ždashed line, after Kieffer, 1983.. V ŽVoragine., NE Žnortheast crater., SE Žsoutheast crater., BN ŽBocca Nova.. TDF: Torre del Filosofo. The square marked as the studied area indicates the outline of Fig. 5.

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

heat flux of the order of 1600 W my2 Žfollowing Ponte, 1927.. All the summit craters were in a steady state of magmatic activity at the time of measurements ŽJune 1988.. A previous geophysical study by Aubert and Baubron Ž1988. showed the existence of an active but hidden convective fissure in the study area, only detected by Self-Potential measurements ŽSP. and gas analysis. Tabbagh and Trezeguet Ž1987.

417

made measurements of l and a Žthermal diffusivity. at a station located at 70 m from C ŽFig. 5. in the direction 558N, called station T further below. 3.2. SP and T mapping Spatial correlation between SP and T maps have been established from many studies ŽZohdy et al.,

Fig. 5. SP and T maps Žlocated on Fig. 4.. The contours Žinterval 20 mV for SP, 58C for T . are drawn from data obtained on a grid, 5 m in the short side for both SP and T ; in the long side Ždirection N 1208. 2 m for SP and 4 m for T. Side scale bar is in meters. Detailed measurements on profile ABC are shown Fig. 6.

418

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

1973; Corwin and Hoover, 1979; Aubert and Baubron, 1988; Aubert et al., 1984.. These authors explain this positive correlation because SP anomalies can be produced by a streaming potential, as a consequence of fluid pressure gradients; by a thermoelectric potential, as a consequence of a temperature gradient; by a chemical potential, as a consequence of a chemical gradient; or by these effects possibly acting together. Hydrothermal activity Žtwo-phase water and steam flow. may create ground temperature anomalies, as it is demonstrated in this paper, and then SP anomalies by both an electroelectric effect and a thermoelectric effect. Because of this correlation, an SP map is used to select the locations for T measurements. This allows efficient data collection, because the acquisition rate for SP measurements is faster when compared with the T measurements Žin the order of 10 s for each SP measure, 5 min for each T measure.. SP measurements were made at shallow depth Ž- 5 cm. with a high-impedance voltmeter Ž10 10 V m. and all values were referenced to the same station. The dimensions of the grid were 2 = 5 m, the local mean noise being estimated at 5 mV Žmean value of the difference between measured value at one station and the value used to establish the SP map ŽFig. 5.. T measurements were made at a depth of 30 cm with thermistors Žresolution 0.18C, precision - 0.58C.. The dimensions of the grid were 4 = 5 m. The overall precision of T measurements depends on a number of factors. These essentially are: precision of the depth when the vertical gradient is high, and time interval between placing of the sensor and T acquisition, especially at high temperatures because the steam flux path may be modified by the digging of the placing hole. This time interval was at least 20 min in this study. The error of each measurement is of the order of 0.28C, but the error resulting from an error of position may reach 0.48C. SP and T maps are presented in Fig. 5. The contours are drawn from data obtained on a rectangular grid: for SP, 2 m in the long side direction Ž1208N. and 5 m in the short side; for T, 4 m = 5 m. In agreement with previous quoted studies, there is a good qualitative correlation between the SP and T maps, particularly on an elongated anomaly, oriented 3458N, which extends for a distance of about 40 m. This direction corresponds to the mean strike

of the Vulcarolo fissure. The length of this fissure allows me to consider it as a 2D structure in the Ž X, Z . plane, with X perpendicular to the lengthening direction, and Z along the vertical axis. To refine the correlation between SP and T values, high density measurements were made along the ABC profile, with a sampling interval of 2 m for SP and T ŽFig. 6.. In order to evaluate the validity of the 1-D assumption for the T field at station B, T measurements at depth of 30 cm were made around B using a 0.5 = 0.5 m grid ŽFig. 7.. This assumption is discussed further below. 3.3. Diurnal and seasonal Õariations Theoretically, the mapping of a steady-state-T field in the XY plane or in the XZ section requires in principle continuous T recordings to adjust for the diurnal and seasonal variations. I estimate day-variations at a depth of 0.3 m from a model where the thermal diffusivity a is assumed equal to 6 = 10y7 m2 sy1 from measurements made by Tabbagh and Trezeguet Ž1987.. For measurements made between 9 and 13 h, Table 1a indicates that the half-range of variation is less than 0.28C, i.e., less than the experimental errors. Thus, continuous T recordings were not used because we made measurements for the T

Fig. 6. Detailed profile ABC Žlocation given on Fig. 5.. Closed circles: SP measurements made at surface. Triangles: T measurements made at depth of 0.3 m. Measurement interval is 2 m.

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

419

Table 2 Field of vertical temperature T Žin 8C. field on profile AB, as a function of X Žin m. and Z Žin cm.

Fig. 7. Detailed T map around station B ŽFigs. 5 and 6.. T Žin 8C. is measured at a constant depth of 0.3 m; contour interval is 18C; measurement grid is 50=50 cm.

mapping only in this time interval. Table 1b indicates theoretical seasonal variations of d Trd Z. These variations should be considered when dTrdZ is monitored. 3.4. T-Õertical field One the same profile ABC but limited to section AB, we measured eight vertical T profiles at a

X Žm.

0.0

8.2

12.2

16.5

19.0

21.4

24.0

26.2

Z Žcm. 10 20 30 40 60 80 100

7.4 8.6 11.0 12.2 14.5 16.9 22.9

8.1 10.6 14.0 16.2 22.2 26.7 30.1

12.4 14.6 17.6 22.8 29.7 35.0 37.7

14.5 18.3 24.0 28.9 40.4 45.4 48.7

16.5 19.0 26.0 31.5 44.1 47.0 48.6

21.0 26.2 33.1 38.5 43.4 46.9 50.6

21.3 28.9 32.9 37.3 41.5 44.8 51.0

30.7 40.9 50.8 53.7 56.9 56.5 57.4

maximum depth of 1 m in order to obtain the T field in the XZ plane ŽTable 2.. The error is estimated to be 0.28C, as for the T mapping, and to be 1 cm for thermistor position. This error of position produces an error in T which depends on the vertical gradient value. Except near station B, the mean gradient value is less than 508C my1 . An error of 1 cm thus produces an error of 0.58C. From these measurements, we can obtain the T field in the vertical plane ŽFig. 8.. Vertical T measurements were also made at the station Ref. 1, located 60 m north of A ŽFig. 4., in order to obtain a cold reference. Fig. 9 indicates T profiles in relation to depth for 4 characteristic stations of profile AB Žcold, intermediate and hot. and station Ref. 1.

Table 1 Theoretical model of ground-temperature variations adapted to the case of Etna Žaltitude 2900 m. Assumed thermal diffusivity value: 6 = 10y7 m2 sy1 . Ža. Day-variation in June. Air-temperatures vary between 48C and 208C Žday-mean temperature 128C.. ŽT0.3 y TM .: deviation Žin 8C. between the instantaneous value of T at depth 0.3 m and the day-mean value of T at the same depth; dTrdZ: ŽT0.3 y T0.1 .rŽ0.3 y 0.1. in 8C my1 . Inside the interval time 8.00–10.00 h, dTrdZ values does not overrun 58C my1. Žb. Annual variation. Air-temperatures vary between y58C and 158C Žannual-mean temperature 58C.. dTrdZ: ŽT0.8 y T0.3 .rŽ0.8 y 0.3. in 8C my1 , calculated for the mid-month time. Ža. Local time Žh.

7

8

9

10

13

16

19

T0.3 y T M Žin 8C. dTrdZ Žin8C my1 .

y0.5 q10.5

y0.6 q5.0

y0.7 q1.0

y0.8 y4.0

y0.5 y15.5

0 y18.5

q0.5 y10.5

Žb. Mid-month

January

March

May

June

July

September

November

dTrdZ Žin8C my1 .

q2.4

y1.3

y4.4

y4.4

y2.6

q1.8

q4.2

420

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

Fig. 8. Vertical temperature T field on profile AB. The filled circles indicate measurement points. T isoline interval is 48C. Vertical exageration= 10.

3.5. Water Õapor measurements The instrument ŽFig. 10. is devised to measure vapor flow at temperature - 1008C and without non-condensible gases. It is made of stainless steel and has a cylindrical shape. It consists of: a cylindrical base ŽA, Fig. 10. through which the steam is conveyed into the upper conical section, inside the condensation cone. There, steam condenses on the internal wall that is kept cold by ice cubes contained in the tank B. This tank is isolated from air and

Fig. 9. Representative vertical T profiles at four stations, from the coldest ŽRef. 1. to the hottest Žstation B..

surrounding ground by a empty tank C Žlow pressure.. The condensed water flows down well polished internal walls of the cone, and is drained through a gutter D where it is collected and measured. The efficiency of the instrument is 100%, so long as tank B contains ice cubes and the steam temperature is Ts ) 458C. It diminishes when Ts

Fig. 10. Instrumental section of the device used for the fumarolic vapor steam condensation. ŽA. condensation cone; ŽB. cooling tank, volume 2 l; ŽC. empty isolating tank; ŽD. gutter to collecting condensate, volume 0.1 l. The section through which the steam flows has a surface area of 200 cm2 .

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429 Table 3 Water vapor flow measurements at station B and supplementary data Depth Condensate T Žm. Ž8C. 10y6 kg my2 sy1 0.3 0.8

85 130

Water Water Water PH percent percent percent Žweight. Žsatur.. Žrelative.

52.9 7.9 57.5 7.9

10.3 14.0

76 56

7.1 6.8

diminishes. Italiano et al. Ž1998. describe an instrument like this but used with other conditions Žtemperature ) 1008C and presence of uncondensible gases.. We made two measurements of vapor flow-rate at depths of 30 and 80 cm, one at point B, where the efficiency of the instrument is equal to 100%, and the second at the cold reference station ŽRef. 1., where efficiency is not known but is probably ) 20%. Table 3 reports the data obtained at station B Žat Ref.1, the condensate values are equal to zero., and the percentage water content in the ground. Note that the ground porosity is not exactly the same at depths of 30 and 80 cm because of the different percentage values of water content at saturation. The flow measured at a depth of 80 cm corresponds to a water velocity of 130 = 10y6 mm sy1 , which also represents a water height equivalent to 4000 mm yeary1 . Rainfall is estimated at 500 mm yeary1 at this site ŽCosentino, 1974.. Hence, we assume that the ascending water flow can originate in two parts: first, downward condensed water flow; secondly, rain-flow infiltrated and vaporized at depth in the soil surrounding the fissure which acts as a drain. 4. Data analysis and estimation of the total heat flux 4.1. SP and T maps The present study is not concerned with the quantitative analysis of the relationship between T and SP. Nonetheless, we observe a very good correlation between these two parameters in the XY plane ŽFig. 5., especially when considering only the maximum values of SP and T which both correspond to the direction 3458. This result further supports the hypothesis that SP measurements are an effective method for mapping convective areas ŽCorwin and

421

Hoover, 1979; Aubert et al., 1984.. If we consider SP and T measurements along the ABC profile perpendicular to the direction 3458 ŽFig. 6., the linear correlation coefficient is 0.9, and the ratio between average values of SP and T is about 4.5 mV 8Cy1 . In my opinion, this indicates that the SP anomaly is not produced by only a thermo-electric effect, producing values of about 0.2 mV 8Cy1 ŽCorwin and Hoover, 1979., but by both thermoelectric and electro-kinetic effects. A previous experimental study ŽAntraygues and Aubert, 1993. showed that both the SP and T anomalies are linked to the vapor flow. 4.2. Calculation of Q r l Õalues field on profile AB From Eq. Ž2., Qrl values are calculated along the ABC profile using the Laplacian function of the T field. Computation is carried out in a two step process: Ži. smoothing using a B-spline function ŽFortran, 1985. on a single variable Z for each X value ŽT is assumed to be only a function of Z ., then Žii. smoothing of the Laplacian so obtained using a B-spline function on two variables X and Z over the entire domain. To maintain the same weighting on all measurements, the computation does not take account of the supposed symmetrical plane passing through B and the limit conditions above and below the studied field. Table 4 gives all derived Qrl values from a depth of 5 to 105 cm Žusing slab thicknesses of 10 cm., and Qrl values for the whole 15–95 cm Table 4 Vertical Qr l values Žin 8C my2 ., as a function of X Žin m. and Z Žin cm. The values are computed assuming a slab thickness of 10 cm. X Žm.:

0

8.2

12.2 16.5 19.0 21.4 24.0 26.2

Interval Žcm. 5–15 q5 y2 y2 y9 y6 y73 15–25 y3 y17 y4 y14 y36 y39 25–35 q3 q6 0 y8 y15 y11 35–45 q4 q8 q3 y2 q1 q10 45–55 q3 q8 q5 q3 q13 q24 55–65 0 q6 q5 q7 q20 q30 65–75 y3 q3 q6 q11 q22 q28 75–85 y8 y3 q5 q13 q19 q17 85–95 y14 y11 q3 q14 q9 y4 95–105 y22 y2 q1 q14 y6 y34 15–95 y18 0 23 24 33 55

q24 q21 q19 q15 q12 q9 q5 q1 y3 y8 79

q93 q76 q60 q45 q31 q17 q5 y7 y18 y27 209

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

422

Fig. 11. Vertical heat field of Qrl values Žin 8C my2 . along profile AB, computed from T values ŽTable 2., and limited to the 30–80 cm depth interval. Full circle: measurement points.

interval, excluding the first Ž5–15 cm. and the last Ž95–105 cm. lines where the accuracy is decreasing Žsee above.. Fig. 11 shows the Qrl field for the 30 to 80 cm interval where the accuracy is assumed to be homogeneous. According to the previous theoretical model ŽFig. 2., the condensation rate is highest where Qrl is at a maximum, i.e., at the roof of a convective cell. Figs. 8 and 11 suggest that this cell may consist of two encased structures: a broad structure Žhalf-breadth about 10 m, roof located at a depth of 0.5 m., and a more narrow Žstations x s 24.0 and 26.2. and more active structure, with a roof located at the surface. Table 5 gives supplementary values for stations 24.0 m and 26.2 m where latent heat is produced as far as the surface Žinterval 0–5 cm.. These values obtained by extrapolation make it possible to calculate the sum Qrl in the 0–95 cm interval. These

Table 5 Supplementary computed Qr l values at stations x s 24.0 and x s 26.2 m, for three slabs: 0–5, 0–95 and 30–80 cm Interval Žcm.:

0–5

0–95

30–80

x Žm. 24.0 26.2

14 55

117 357

51 127

values are then used to estimate the flux Žsee Table 7.. Table 5 also indicates values obtained for the 30–80 cm interval, to be compared at station 26.2 m with direct measured values of Q Žsee Table 6.. 4.3. Particular analysis of data at station B (x s 26.2 m) At station B, the flow-rate measurements allow us to evaluate the range of the convective term and the validity of the computation of Qrl. 4.3.1. Validity of 1D approximation The 1D assumption around station B needs to be checked because the presence of one fissure is in Table 6 Latent heat Q produced at station B within three intervals: 0.0–0.3, 0.3–0.8 and 0.0–0.0.8 m D: water vapor flow measured by condensation; Ža. at depth of 0 m, D is assumed to be equal to zero; Q1: heat computed from D measurements; Q2 : heat computed from Qr l values, with l s 0.8 W my1 8Cy1 . Depth Žm.: y6

y2

y1 .

D Ž10 kg m s Q1 from D ŽW my2 . Q2 from Qr l, l s 0.8 ŽW my2 .

0.0

0.3

0.8

0 Ža. ²200:²105:²305: ²203:²100:²303:

85

130

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

principal inconsistent with this assumption. From map Fig. 7, I calculate the second derivative of T along the two horizontal axes Ž708N and 1608N. at a depth of 0.3 m, yielding values of y8 and 08C my2 corresponding to 708N and 1608N, respectively. These values should be compared with the value obtained along the vertical axis Žy608C my2 given in Table 4 for the slab thickness of 0.1 m at a depth of 0.3 m.. Thus, the 1D approximation at station B is valid as regards the Qrl computation, being all the more applicable for the other stations where the T and Qrl fields are more regular. 4.3.2. EÕaluation of conÕectiÕe flux In an assumption of 1D heat flow, the convective flux is equal to wŽ dTrdZ .C v r V x. The gas is a mixture of saturated water vapor and air in unknown proportions, with other gases assumed to be present in trace amounts. We look for a rough estimate to show that the convective flux is normally smaller than the latent heat flux. Then we estimate each term of this flux: dTrdZ - 1008C my1 , C v can be assumed to be 1500 J kgy1 Ky1 and r is taken as 1 kg my3 Ž1.293 for air, 0.9 for water.. However V is more difficult to evaluate. The water flow rate measured at 80 cm Ž130 = 10y6 kg my2 sy1 . corresponds to a velocity of 130 = 10y6 mm sy1 in free space, a value that must be multiplied by a porosity factor which is taken as 10. At a temperature of 508C, the partial pressure of water is equal to 95 mm Hg. By assuming that the total vapor pressure is equal to atmospheric pressure Ž526 mm Hg at altitude 3000 m., the air pressure can be calculated as six times higher than the water vapor pressure. At most, we can assume that a maximum gas flow velocity is 60 times higher than the initial water flow-rate, that is to say, a velocity of approximately 10y5 m s y 1. Under these conditions, the maximum convective flux may be estimated as 1.5 W my2 . 4.3.3. Validity test of Q r l computation Following the assumptions given previously, the Q value derived from Qrl calculated in the 30–80 cm interval Ž1278C my2 , from Table 5. should correspond to the Q value obtained from flow difference measured at depths 80 and 30 cm. Because of the very low flow velocity, we assume that this flow is

423

mono-phasic, and that the flow difference Ž45 = 10y6 kg my2 sy1 . corresponds to the latent heat, 105 W my3 . The l value obtained is l s 105r127s 0.83 W my1 Ky1 . This value is higher than the value Ž0.53 W my1 y1 . K obtained by Tabbagh and Trezeguet Ž1987. at station T Žsee over. located on the same scoria cone, outside the fissure. At this station, the flux is only conductive and equal to 0.75 W my2 Two reasons explain a high value of l in this context ŽWoodside and Messmer, 1961; Baladi et al., 1981.: increasing of water content, and context of moisture migration phenomena. Without l and Q measurements at all the stations, we assume a value of l rounded to 0.8 W my1 Ky1 . Further investigations would be necessary to define l with more precision in the context of unconsolidated media with three phases, solid, liquid, gaseous. Table 6 compares the values of Q produced in the 0–30 cm and 30–80 cm intervals computed by the two methods: Q1 derived from flow measurements, and Q 2 derived from Qrl values. For the interval 0.3–0.8 m, the difference is due to the assumed rounded value of l. For the interval 0.0–0.3 m, the good agreement between these two values is reasonable grounds for validating the assumption that flow becomes equal to zero at the surface and the method for Qrl computation. In consequence of this result, it is probable that the F value computed by the Qrl method is underestimated when the F value exceeds 286 W my2 , because the part of Q produced in the atmosphere is not taken into account. 4.3.4. Other Õalidation tests Two other tests are used to check the Qrl computation. Firstly, the Qrl values derived by the method described above are shown on Fig. 12, along with the results of a non-smoothed computation of d 2 Trd Z 2 from measured T values using Taylor approximation, for a slab thickness of 20 cm. The second method clearly leads to a greater variability of the values than the smoothed method, but the sum of Qrl values between 30 and 80 cm is nearly the same in both methods Ž1328C my2 for the direct method, 1278C my2 for the smoothed one.. Secondly, the extrapolated value toward the surface obtained by the smoothed method yields a surface temperature of 11.28C, while the extrapolated surface

424

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

Fig. 12. Qr l variation at station B versus depth, computed from d2 T r dZ 2 using two different methods: Ž1. continuous line with smoothed values by spline function; Ž2. vertical segments computed for each slab Žslab thickness is 20 cm. by Taylor approximation.

temperature at Ref. 1 is equal to 88C. This discrepancy of 3.28C seems to be realistic, according to the maximum of condensation rate observed at the surface, but without water flowing across the groundrair surface. 4.4. Computation of heat discharge on profile AB Table 4 clearly shows two features of the Qrl field, according to the location of the stations along the profile. For the two hottest stations Ž x s 24.0 m and 26.2 m., Qrl increases with decreasing depth, and heat discharge can be only estimated from heat produced in the medium. For all the other stations, the data of Table 4 indicate that the top of the convective cell does not reach the surface, thus conductive transfer can then be computed. Table 7 and Fig. 13 present heat discharge values for all stations. These come from the data in Table 4 obtained as follows: Ž1. from 0 to 21.4 m by conductive transfer, and Ž2. from 0 to 26.2 m by heat produced in the medium. From station 21.4 to station 0, the discrepancy increases between the two methods of calculating heat discharge. To explain this discrepancy, we assume that Ž1. a part of the latent heat production may be located below a depth of 1

m, and Ž2. conductive heat transfer takes place from the deep hot zone towards the surface. Table 1 indicates the theoretical value of dTrdZ computed from T values in June in the interval 9 h–13 h where measurements were made. The apparent computed fluxes with a l value of 0.8 W my1 Ky1 vary from y2.7 to y6.7 W my2 . These rough diurnal and seasonal corrections are not made in the present case. Nevertheless, it indicates that the minimum computed value of F at station x s 0 m Ž6 W my2 . is really a minimum value, higher than the value 0.72 W my2 at station T computed by Tabbagh and Trezeguet Ž1987.. Fig. 14 indicates heat flux values F in relation to T at a depth of 0.3 m. This relation is clearly non-linear in the interval 0–100 W my2 ; it seems to be linear in the interval 100–300 W my2 and surely becomes non-linear to higher F values, because T values are limited by the water boiling temperature at atmospheric pressure Ž918C at altitude 2900 m.. 4.5. Total heat flux Õalue on the area corresponding to T map (Fig. 5) We assume that the relationship between T and F indicated in Fig. 14 is available on the whole surface on this map. The total flux value is the sum of F values corresponding to each surface limited between two iso-T lines Žinterval 58C., where the mean F value is estimated from data indicated in Fig. 14 Žfor example, the mean flux value in the interval 25–308C is 56 W my2 .. Then the total heat flux is equal to about 210 kW, on a whole surface equal to 2100 m2 . Table 7 Heat discharge along profile AB, computed by two different methods By conductive transfer dT r dZ from x s 0 to x s 21.4, and by latent heat Q for stations x s 24.0 and 26.2. l is taken as equal to 0.8 W my1 8Cy1 for the whole profile. X Žm.:

0

8.2 12.2 16.5 19.0 21.4 24.0 26.2

dT r dZ 12 30 40 Ž8C my1 . Qr l Ž0–95 cm. y18 0 23 Ž8C my2 . F Ž l s 0.8. 10 24 32 ŽW my2 .

54

60

60

–

–

24

33

55

117 357

43

48

48

94 286

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

425

Fig. 13. Heat flux F Žstar. computed on profile AB, with ground temperature Žtriangle. at depths of 0.1 m ŽT10 . and 0.3 m ŽT30 ..

4.6. Error eÕaluation made on flux Õalue 4.6.1. Error on absolute flux Õalue The error on the absolute flux value depends on two terms: first the error in dTrdZ or in Q, secondly

the error in l. Let us consider separately the fluxes obtained by the dTrdZ and by the Qrl computation. Ža. Error in dTrdZ. Each T and Z measurement is made with an absolute error less than 0.28C and 2

Fig. 14. Heat flux F computed on profile AB, versus temperature T at depth of 0.3 m. Extrapolation to higher T is non-linear because T is limited at 918C Žboiling T at altitude 2900 m..

426

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

cm, respectively. T values obtained near the ground surface are more inaccurate because of the diurnal variation. This effect can be corrected, either by recording the diurnal variation, or assuming this variation from a theoretical model. When the flux is computed from at least four values, as on profile AB, we estimate that the relative error in the flux is less than 10%, except at station x s 0 where the error appears to be greater Ž20%. because of the low T variation versus Z. Žb. Error in Qrl. For Qrl, such a method of error estimation is too approximate considering the number of processing steps involved. When Qrl can be checked by direct flow measurement, and when this test confirms the calculated Qrl values, we consider that this method gives a Qrl value with an error less than 20%. Žc. Error in l. l may be estimated from direct measurement ŽTabbagh, 1985; Tabbagh and Trezeguet, 1987. when steam flux is insignificant, or from vapor steam flow measurement combined with Qrl computation as is demonstrated when F ) 100 W my2 ŽTable 7.. In the two cases, we assume that the error in l is about 10%. Finally, we estimate the error in the absolute value of F between 20 and 30%. 4.6.2. Error in relatiÕe Õalue of f The relative error should be considered when the objective is thermal monitoring, i.e., measurement of time variations of F Žmeasurement taken at the same place with the same experimental device.. In this case, only T measurements should be considered, and the error becomes less than 10%.

5. Discussion 5.1. Practical method for heat flux measurement On the basis of the above results, we can propose a practical method for estimating the heat flux on an active hydrothermal convective cell: Ž1. SP mapping should be first carried out at the surface in order to delimit broadly the active zone. Using the SP map as a guide, T mapping at a constant depth, which is a compromise between shallow sampling depth to facilitate easy measurements

and sufficient depth to minimize diurnal variation, should then be carried out. From the T field described in this study, it appears that the suitable sampling depth is located between 20 and 30 cm. Ž2. At a selection of stations, ranging between the hottest and the coldest, vertical T profiles should be measured in order to calculate conductive and convective heat fluxes at these stations. When F values vary in the range 10–300 W my2 , it seems that measurement depth equal to about 1 m is sufficient to estimate the flux with an error less than 30%. Flux may be calculated by dTrdZ measurement in the range 10–100 W my2 , and by Qrl in the range 100–300 W my2 . A borehole to 1 m depth requires granular ground throughout. If it is not the case, dTrdZ measurement is possible with a borehole depth of 0.3 m, and a F value may be calculated only in the range 10–100 W my2 . Supplementary data Ž l and steam flow measurements. are necessary to calculate the absolute value, but not the relative value. Ž3. When flux values exceed 300 W my2 , it is probably no longer valid to assume non-convective flow across the surface. In such cases flux values may be underestimated, and direct flow rate measurements are required to obtain a correct computation. Ž4. When the flux is less than 10 W my2 , the accuracy of the evaluation decreases, firstly because the T difference decreases between depths Z s 10 and 100 cm, and secondly because this difference approaches that due to diurnal and seasonal variations. Thus, it becomes necessary to record diurnal variations at one station, or to set an accurate theoretical model to correct for these effects. Ž5. To calculate the Qrl field with accuracy, a mathematical approach is required for the computation of heat fluxes. For a preliminary evaluation in the Qrl field, Eq. Ž5. gives a rough value. At station B, in the 10–60 cm interval, this value is equal to: 2

Qrl s 2 Ž T0.6 y T0.1 . r Ž 0.6 y 0.1 . s 210 W my3 , This rough value is underestimated with an error about 20% Žthe accurate value found above is 267 W my3 ., and the sense of the error is always the same.

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

5.2. Time-Õariations of heat discharge It is quite clear that heat fluxes computed from vertical T profiles from 5 cm below the surface to a depth of 1 m may be attributed to two effects: Ži. heat from deep zones, and Žii. heat exchange between atmosphere and ground which is influenced by, for example, diurnal and seasonal T variations Žsee Table 1., and precipitation. The relationship between these two effects has been recently described ŽYasuaki and Hurst, 1998.. Only a comparison between stations, located inside and outside the hydrothermal area but subject to the same microclimate, allows us to distinguish heat variations originating from external and internal effects. 5.3. Calibration of remote sensing measurements by ground-T measurements For flux values between 10 and 300 W my2 , it is possible to measure heat flux on a surface covering a relatively large pixel, for example 100 m = 100 m. Subtle thermal anomalies at fumarolic and subfumarolic areas can also be detected and analysed using remotely sensed data Žsee for example Gaonac’h et al., 1994; Harris and Stevenson, 1997a,b.. It would be interesting to compare the results obtained by the two methods on the same site and approximately at the same time, to determine how these two approachs to thermal phenomena could be complementary for heat flux mapping and monitoring. 5.4. Importance of hydrothermal transfer Although this part is speculative, it is interesting to compare heat flux measured at the surface and heat flux provided from different magma sources. We calculate Ži. the volume of dyke beneath the surface of the studied fissure necessary to provide a magmatic heat flux equal to the measured heat flux at the surface, and Žii. the cumulate hydrothermal flux corresponding to heat flux produced by known lava flows. 5.4.1. Volume of dyke beneath the fissure On the half-profile ABC Žsection AB., the mean flux is equal to 73 W my2 , corresponding at an energy per year equal to 2.3 = 10 9 J my2 . The heat

427

energy produced by a dyke model ŽHarris and Stevenson, 1997b. characterized by a crystallised fraction equal to 45% and by a cooling equal to 3008C is equal to 1.2 = 10 9 J my3 . Then one may imagine a body beneath the Vulcarolo fissure, having the same length and breadth Ž52 m. and a vertical thickness which is cooling and crystallising each year equal to 2 m Ž2.3r1.2.; or a more realistic body having a breadth equal to 10 m, as the active part of the fissure. Measurements were made in 1988. If the heat flux at the surface is supposed constant since the re-opening of the Vulcarolo fissure Ž1971 after Kieffer, 1983., the total cooled thickness of this supposed body would be equal in 1988 to about 180 m Ž17 = 2 = 5.2.. Kieffer Ž1983. indicated that this fissure, as the other south-east fissures of this zone, participates to the distension tectonic activity which affects the central parts of Etna. 5.4.2. Cumulate hydrothermal flux corresponding to heat flux produced by known laÕa flows Barberi et al. Ž1993. indicate that the total volume of lavas from 1971 to 1989 was equal to 435.7 = 10 6 m3 , that corresponds to a mean volume of 23 = 10 6 m3 per year. The heat produced by crystallisation and by cooling equal to 8508C is equal to 3 = 10 9 J my3 . Then the heat provided by these lava flows is equal to 7 = 10 16 J per year, equivalent to a flux equal to 2000 MW. Three sorts of convective areas can provide hydrothermal flux over Etna. First, on the summit cone Žsurface area 2 km2 , from altitude 2900 m to the summit., Aubert and Dana Ž1994. measured on a profile located on the south side of the cone a mean ground temperature at a depth of 0.3 m equal to 528C, that corresponds to a mean flux of the order of at least 280 W my2 . If we assume that this profile is representative of the whole cone, the total thermal power would be equal to 560 MW. Secondly, the system of fissures on Etna ŽKieffer, 1983. is well developed around the summit at a radius of about 6 km. This system could provide a total length of fissures similar to the studied section of the Vulcarolo fissure equal to at least 10 km, that represents a total flux of about 40 MW. Third, the third part could be represented by a large zone around the cone Žradius 3.5 km, surface

428

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429

area 40 km2 . where the mean flux value could be equal to about 0.7 W my2 , the value measured by Tabbagh and Trezeguet at station T where the flux is only conductive. The total flux provided from this surface would be equal to about 30 MW. Finally, the total flux value provided from these three sorts of hydrothermal zones would be equal to 600 MW, i.e., 30% of the flux provided by lava flows. This heat flux needs a cooling magma volume of 0.2 m3 sy1 , a value on the same order that the values obtained by Harris and Stevenson Ž1997b. in the cases of Vulcano and Stromboli. Allard Ž1997. gives a time average Ž1975–1995. figure of 4.5–8 m3 sy1 of degassed magma to account for degassing at the Central Craters of Etna. These estimations are largely approximate but realistic. A direct consequence of them is that non-erupted magma bodies could have a similar volumetric importance as the erupted magma. Kagiyama Ž1981. also obtains the same conclusion for the Japanese Island Arc. 6. Conclusions Ground-temperature measurements near the surface may provide three sorts of results: Ža. identification of hidden thermal fissures which arise from the internal structure of the volcano. The breadth of these fissures can be very narrow, only a few meters ŽAubert et al., 1984.. These results can only be obtained with difficulty by other methods; Žb. monitoring of heat flux, for a volcano monitoring method; Žc. thermal calibration of remote sensing methods by mapping an area corresponding to a pixel and the surrounding zone, with a error less than 30% for an absolute value, 10% for a relative value. The method is easy to use when the superficial ground is granular, such as on a scoria cone, to a depth of 0.3 m for measuring heat flux in the range 10–100 W my2 , or a depth of 1 m for the range 100–300 W my2 . Outside this range 10–300 W my2 , the precision diminishes and the calculation needs supplementary data. Acknowledgements I am grateful for constructive reviews by L. Wilson, A.J.L. Harris and F. Italiano that greatly im-

proved the final manuscript version. I also thank Philippe Antraygues for his help in the field, Gerard Touchard for computation process, Aleth Prud’Homoz for the drawings, and M.S.N. Carpenter for correcting the English.

References Allard, P., 1997. Endogeneous magma degassing and storage at Mount Etna. Geophys. Res. Lett. 24 Ž17., 2219–2222. Antraygues, Ph., Aubert, M., 1993. Self-Potential generated by two-phase flow in a porous medium: experimental study and volcanological applications. J. Geophys. Res. 98 ŽB12., 22273–22281. Aubert, M., Auby, R., Bourlet, Y., Bourlet, F., 1984. Contribution a` la surveillance de l’activite´ de l’Etna a` partir de l’etude des ´ zones fumerolliennes. Bull. Volcanol. 47-4 Ž2., 1039–1050. Aubert, M., Baubron, J.C., 1988. Identification of a hidden thermal fissure in volcanic terrain, using a combination of hydrothermal convection indicators and soil–atmosphere analysis. J. Volcanol. Geotherm. Res. 35, 217–225. Aubert, M., Dana, I.N., 1994. Interpretation des profils de polari´ sation spontanee ´ ŽPS. en volcanologie. Possibilites ´ d’application de la methode PS a` la surveillance des volcans actifs. ´ Bull. Soc. Geol. ´ France 165 Ž2., 113–122. Barberi, F., Carapezza, M.L., Valenza, M., Villari, L., 1993. The control of lava flow during the 1991–1992 eruption of Mt Etna. J. Volcanol. Geotherm. Res. 56, 1–34. Baladi, J.Y., Ayers, D.L., Schoenhals, R.J., 1981. Transient heat and mass transfer in soils. Int. J. Heat Mass Transfer 24, 449–458. Corwin, R.F., Hoover, D.B., 1979. The self-potential method in geothermal exploration. Geophysics 44 Ž2., 226–245. Cosentino, M., 1974. Il regime pluviometrico sull’ Etna negli anni 1921–1971. Atti della Accademia Gioenia di Scienze naturali in Catania 6, 37–52. Fortran scientific subroutine library 2, 1985. Gaonac’h, H., Vandemeulebrouck, J., Stix, J., Halbwachs, M., 1994. Thermal infrared satellite measurements of volcanic activity at Stromboli and Vulcano. J. Geophys. Res. 99 ŽB5., 9477–9485. Hardee, H.C., 1982. Permeable convection above magma bodies. Tectonophysics 84, 179–195. Harris, A.J.L., Stevenson, D.S., 1997a. Thermal observations of degassing open conduits and fumaroles at Stromboli and Vulcano using remotely sensed data. J. Volcanol. Geotherm. Res. 76, 175–198. Harris, A.J.L., Stevenson, D.S., 1997b. Magma budgets and steady-state activity of Vulcano and Stromboli. Geophys. Res. Lett. 24 Ž9., 1043–1046. Italiano, F., Pecoraino, G., Nuccio, P.M., 1998. Steam output from fumaroles of an active volcano: tectonic and magmatic–

M. Aubertr Journal of Volcanology and Geothermal Research 92 (1999) 413–429 hydrothermal controls on the degassing system at Vulcano ŽAelian arc.. J. Geophys. Res. 103 ŽB12., 29829–29842. Kagiyama, T., 1981. Evaluation methods of heat discharge and their applications to the major active volcanoes in Japan. J. Volcanol. Geotherm. Res. 9, 87–97. Kieffer, G., 1983. Deplacement progressif vers l’Est des voies ´ d’alimentation magmatique de l’Etna ŽSicile.. Bull. Soc. Geol. Fr. 7, 25, Ž3., 369–378. Oppenheimer, C., Francis, P.W., Rothery, D.A., Carlton, R.W., Glaze, L.S., 1993. Infrared image analysis of volcanic thermal features: Lascar volcano, Chile, 1984–1992. J. Geophys. Res. 98 ŽB13., 4269–4286. Ponte, G., 1927. Il vulcarolo sull’Etna e la utilizzatzione del suo vapore aequeo. Bull. Acc. Gioenia Sc. Nat. 57 Ž1–2., 14–17. Sekioka, M., Yuhara, K., 1974. Heat flux estimation in geothermal areas based on the heat balance of the ground surface. J. Volcanol. Geotherm. Res. 79, 2054–2058. Tabbagh, A., 1985. A new apparatus for measuring thermal properties of soils and rocks in situ. IEEE Trans. Geosci. Remote Sens. GE 23-26, 896–900. Tabbagh, A., Trezeguet, D., 1987. Determination of sensible heat flux in volcanic areas from ground temperature measurements

429

along vertical profiles: the case study of Mount Etna ŽSicily, Italy.. J. Geophys. Res. 92 ŽB5., 3635–3644. Williams, H., McBirney, A.R., 1979. Volcanology. Freeman, Cooper. Woodside, W., Messmer, J.H., 1961. Thermal conductivity of porous media: I. Unconsolidated sands. J. Appl. Phys. 32 Ž9., 1688–1699. Wooster, M.J., Rothery, D.A., 1997. Thermal monitoring of Lascar volcano, Chile, using infrared data from the along-track scanning radiometer: a 1992–1995 time series. Bull. Volcanol. 58, 566–579. Yasuaki, S., Hurst, A.W., 1998. Temperature changes at depths to 150 metres near the active crater of Aso Volcano: preliminary analysis of seasonal and volcanic effects. J. Volcanol. Geotherm. Res. 81, 159–172. Yuhara, K., Ehara, S., Tagomori, K., 1981. Estimation of heat discharge rates using infrared measurements by a helicopterborne thermocamera over the geothermal areas of Unzen volcano, Japan. J. Volcanol. Geotherm. Res. 9, 99–109. Zohdy, A.A.R., Anderson, L.A., Muffler, L.J.P., 1973. Resistivity, self-Potential, and induced-polarization surveys of a vapordominated geothermal system. Geophysics 38 Ž6., 1130–1144.