Diffuse degassing and the hydrothermal system at masaya volcano, nicaragua

Diffuse degassing and the hydrothermal system at masaya volcano, nicaragua

Material Information

Diffuse degassing and the hydrothermal system at masaya volcano, nicaragua
Pearson, Sophie
Place of Publication:
[Tampa, Fla]
University of South Florida
Publication Date:


Subjects / Keywords:
Geophysical survey
Gas flux
Groundwater convection
Fracture flow
Dissertations, Academic -- Geology -- Masters -- USF ( lcsh )
non-fiction ( marcgt )


ABSTRACT: Hydrothermal systems change in response to volcanic activity, and in turn may be sensitive indicators of volcanic activity. Fumaroles are a surface manifestation of this interaction. We use time series of soil temperature data and numerical models of the hydrothermal system to investigate volcanic, hydrologic and geologic controls on this diffuse degassing. Soil temperatures were measured in a low-temperature fumarole field located 3.5 km from the summit of Masaya volcano, Nicaragua. They respond rapidly, on a time scale of minutes, to changes in volcanic activity also manifested at the summit vent. The soil temperature response is repetitive and complex, and is characterized by a broad frequency signal allowing it to be distinguished from meteorologic trends. Geophysical data reveal subsurface faults that affect the transport of fumarole gases. Numerical modeling shows that these relatively impermeable faults enhance flow through the footwall. On a larger scale, modeling suggests that uniform injection of fluid at depth causes groundwater convection in a permeable 3-4 km radial fracture zone transecting the entire flank of the volcano. This focuses heat and fluid flux and can explain the three distinct fumarole zones located along the fracture. We hypothesize that the rapid response of fumarole temperature to volcanic activity is due to increased flow of gas through the vadose zone, possibly caused by changes in the subsurface pressure distribution. Numerical models show that an abrupt injection of hot gas, at approximately 100 times background rates, can cause the rapid increase in temperature observed at the fumaroles during volcanic activity. A decrease in hot fluid injection rate can explain the gradual decrease in temperature afterwards. Mixing with surrounding vadose-zone fluids can result in the consistent and abrupt decreases in temperature to background level following hot gas injection. Fumaroles result from complex interaction of the volcanic-hydrologic-geologic systems, and can therefore provide insight into these systems. Increases in fumarole temperature correspond to increased gas flux related to changes in volcanic activity, suggesting that monitoring of distal fumaroles has potential as a volcano monitoring tool, and that fumarole temperatures can provide insight into the response of shallow gas systems to volcanic activity.
Dissertation (PHD)--University of South Florida, 2010.
Includes bibliographical references.
System Details:
Mode of access: World Wide Web.
System Details:
System requirements: World Wide Web browser and PDF reader.
General Note:
Title from PDF of title page.
General Note:
Document formatted into pages; contains X pages.
Statement of Responsibility:
by Sophie Pearson.

Record Information

Source Institution:
University of South Florida Library
Holding Location:
University of South Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
E14-SFE0004512 ( USFLDC DOI )
e14.4512 ( USFLDC Handle )

Postcard Information



This item has the following downloads:

Full Text
xml version 1.0 encoding UTF-8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd
leader nam 22 Ka 4500
controlfield tag 007 cr-bnu---uuuuu
008 s2010 flu s 000 0 eng d
datafield ind1 8 ind2 024
subfield code a E14-SFE0004512
XX9999 (Online)
1 100
Pearson, Sophie.
0 245
Diffuse degassing and the hydrothermal system at masaya volcano, nicaragua
h [electronic resource] /
by Sophie Pearson.
[Tampa, Fla] :
b University of South Florida,
Title from PDF of title page.
Document formatted into pages; contains X pages.
Dissertation (PHD)--University of South Florida, 2010.
Includes bibliographical references.
Text (Electronic thesis) in PDF format.
Mode of access: World Wide Web.
System requirements: World Wide Web browser and PDF reader.
3 520
ABSTRACT: Hydrothermal systems change in response to volcanic activity, and in turn may be sensitive indicators of volcanic activity. Fumaroles are a surface manifestation of this interaction. We use time series of soil temperature data and numerical models of the hydrothermal system to investigate volcanic, hydrologic and geologic controls on this diffuse degassing. Soil temperatures were measured in a low-temperature fumarole field located 3.5 km from the summit of Masaya volcano, Nicaragua. They respond rapidly, on a time scale of minutes, to changes in volcanic activity also manifested at the summit vent. The soil temperature response is repetitive and complex, and is characterized by a broad frequency signal allowing it to be distinguished from meteorologic trends. Geophysical data reveal subsurface faults that affect the transport of fumarole gases. Numerical modeling shows that these relatively impermeable faults enhance flow through the footwall. On a larger scale, modeling suggests that uniform injection of fluid at depth causes groundwater convection in a permeable 3-4 km radial fracture zone transecting the entire flank of the volcano. This focuses heat and fluid flux and can explain the three distinct fumarole zones located along the fracture. We hypothesize that the rapid response of fumarole temperature to volcanic activity is due to increased flow of gas through the vadose zone, possibly caused by changes in the subsurface pressure distribution. Numerical models show that an abrupt injection of hot gas, at approximately 100 times background rates, can cause the rapid increase in temperature observed at the fumaroles during volcanic activity. A decrease in hot fluid injection rate can explain the gradual decrease in temperature afterwards. Mixing with surrounding vadose-zone fluids can result in the consistent and abrupt decreases in temperature to background level following hot gas injection. Fumaroles result from complex interaction of the volcanic-hydrologic-geologic systems, and can therefore provide insight into these systems. Increases in fumarole temperature correspond to increased gas flux related to changes in volcanic activity, suggesting that monitoring of distal fumaroles has potential as a volcano monitoring tool, and that fumarole temperatures can provide insight into the response of shallow gas systems to volcanic activity.
Advisor: Charles Connor, Ph.D.
Geophysical survey
Gas flux
Groundwater convection
Fracture flow
Dissertations, Academic
x Geology
t USF Electronic Theses and Dissertations.
4 856
u http://digital.lib.usf.edu/?e14.4512


Diffuse Degassing and the Hydrothermal System at Masaya volcano, Nicaragua by Sophie C. P. Pearson A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Geology College of Arts and Sciences University of South Florida Major Professor: Charles B. Connor, Ph.D. Sarah Kruse, Ph.D. D iana C. Roman, Ph.D. Ward E. Sanford, Ph.D. Mark Stewart, Ph.D. Date of Approval: April 29 2010 Keywords: monitoring, modeling, TOUGH 2, geophysical survey fumarole temperature, gas flux groundwater convection, fracture flow Copyright 2010, Sophie C. P. Pearson


DEDICATION To my parents, John and Su, and my sisters, Alexa and Cassie In memory of my grandfath er, Wilfred who would have loved to see this


ACKNOWLEDGEMENTS I would like to thank my advisor, Dr. Charles Connor, for giving me the opportunity to work on this exciting project, to explore the world of scientific research, and to visit field sites and conferences around the world His guidance was invaluable. I am gr ateful to my committee members, Dr. Ward Sanford, Dr. Diana Roman, Dr. Mark Steward and Dr. Sarah Kruse for their support, insight and advice. I am also grateful to the Instituto Nicaraguense de Estudios Territoriales, particularly Allan Mo rales, Wilfried Strauch, Virginia T e norio, Martha Herrera and Antonio and to Heather Lehto and the University of South Florida 2007 Volcanology class for their support during fieldwork. Thank you to Jona than Wynn for an interesting discussion on the possible importance of geochemi stry in this study Thank you to Mikel, Alain, Armando and Koji for being great office mates and field companions. Thank you to Brad, John P., Wayne a nd John O. for sharing ideas and support I would especially like to thank Jenn for all of our discussions, both professional and personal. Her friendship through exams, conferences, fieldtrips, classes, dissertation writing, and life in Tampa made all the difference. Many than ks also to Mandie and Leah, for being such great colleagues and friends. I gratefully acknowledge f inancial support from the Volcano Hazards Program of the U.S. Geological Survey I particularly wish to thank the University of South Florida Presidential Fellowship program for the financial backing and freedom th at it provided throughout my PhD.


TABLE OF CONTENTS LIST OF TABLES ................................ ................................ ................................ ............. iii LIST OF FIGURES ................................ ................................ ................................ ........... iv ABSTRACT ................................ ................................ ................................ ...................... vii CHAPTER 1 : INTRODUCTION ................................ ................................ ....................... 1 1.1 Previous fumarole studies ................................ ................................ ................. 3 1.2 The hydrologic system at Masaya ................................ ................................ ... 15 1.3 Volcanic activity at Masaya volcano (adapted from Tenorio et al., 2010) ..... 19 1.4 Numerical modeling ................................ ................................ ........................ 27 1. 4 .1 Governing equations (adapted from Pruess et al., 1999; Pruess, 2004) ................................ ................................ .................. 29 1. 4 .2 Solving the equations ................................ ................................ ........ 32 1. 4 .3 TOUGH2 validation (adapted from Oldenburg and Pruess, 1994) ................................ ................................ .............................. 35 CHAPTER 2 : RAPID RESPONSE OF A HYDROLOGIC SYSTEM TO VOLCANIC ACTIVITY ................................ ................................ ....................... 37 2.1 Introduction ................................ ................................ ................................ ..... 37 2.2 Masaya Volcano ................................ ................................ .............................. 38 2.3 Method ................................ ................................ ................................ ............ 39 2 .4 Results ................................ ................................ ................................ ............. 39 2.5 Discussion and Conclusion ................................ ................................ ............. 45 CHAPTER 3 : INTEGRATED GEOPHYS ICAL AND HYDROTHERMA L MODELS OF FLANK DEGA SSING AT MAS AYA VOLCANO, N ICARAGUA ................................ ................................ ................................ ....... 48 3.1 Introduction ................................ ................................ ................................ ..... 48 3.2 Methods for Delineating Flow Paths ................................ .............................. 49 3.3 Geologic Setting ................................ ................................ .............................. 51 3.4 Data Collection and Results ................................ ................................ ............ 54 3.5 Modeling and Interpretation ................................ ................................ ........... 59 3.5.1 Magnetic data ................................ ................................ .................... 59 3.5.2 Hydrothermal m odels ................................ ................................ ........ 60 Vadose zone ................................ ................................ ....... 60 3.5.2 .2 Saturated zone ................................ ................................ .... 67


ii 3.6 Di scussion ................................ ................................ ................................ ....... 72 3.7 Conclusions ................................ ................................ ................................ ..... 76 CHAPTER 4 : NUMERICAL MODELING OF FLANK FUMAROLE TEMPERATURE VARIATIO NS RELATED TO VOLCAN IC ACTIVITY ................................ ................................ ................................ ............ 78 4.1 Introduction ................................ ................................ ................................ ..... 78 4.2 Geologic s etting ................................ ................................ .............................. 80 4.3 Fumarole t emper ature d ata ................................ ................................ ............. 83 4.4 Numerical modeling ................................ ................................ ........................ 86 4.5 Results ................................ ................................ ................................ ............. 91 4.6 Discussi on ................................ ................................ ................................ ....... 97 4.7 Conclusion ................................ ................................ ................................ .... 103 LIST OF REFERENCES ................................ ................................ ................................ 105 APPENDIX A: Finite difference model of 2D heat transfer .... .................................. 1 1 7 APPENDIX B: Fumarole temperature data plotted by year ABOUT THE AUTHOR END PAGE


iii LIST OF TABLES Table 3.1. Statistical results from magnetic, CO 2 and SP surveys ................................ ................. 57 Table 3.2. Parameters used in numerical models ................................ ................................ ........... 63 Table 3.3. Parameters used i n Rayleigh number calculations. ................................ ....................... 71 Table 3.4. Parameters used to determine whether convection will occur along the fracture at Masaya vol cano ................................ ................................ ................................ ............ 72 Table 4 1 Parameters used in TOUGH2 model ................................ ................................ ............ 91


iv LIST OF FIGURES Figure 1.1. Fumarole temperature measurements recorded at a) Telica volcano and b) Cerro Negro volcano, Nicaragua in 2007 ................................ ................................ ........... 9 Figure 1.2 CO 2 flux measured at Comalito cinder cone during this study ................................ .... 12 Figure 1.3. a) Topographic map of Masaya volc ano. b) Elevation of the water table as estimated from TEM soundings ................................ ................................ ........................ 17 Figure 1. 4 Groundwater distribution at Masaya in meters above sea level, according to MODFLOW2000 results. ................................ ................................ ................................ 18 Figure 1.5. Estimated water table depths. ................................ ................................ ...................... 19 Figure 1.6 Timeline of crater activity at Masaya volcano during the study period ...................... 20 Figure 1.7. Real time Seismic Amplitude (RSAM) and total number of earthquakes recorded per month during the study period. ................................ ................................ .... 21 Figure 1.8. Fumarole temperature recorded about once a month by INETER sci entists ............... 21 Figure 1.9. Sulfur dioxide flux from the active crater during the study period.. ........................... 22 Figure 1.10. Images of the new lava lake that developed in the crater in October 2006. .............. 23 Figure 1.11. On average temperatures decreased slowly at the 2006 vent, even after a rockfall blocked it in early 2008. ................................ ................................ ...................... 23 Figure 1.12. Degassing from the vent formed in 2006 ................................ ................................ .. 24 Figure 1.13. New fumaroles formed to the east of the active crater in October 2007. .................. 25 Figure 1.1 4 Strong, consistent degassing in September 2008 ................................ ....................... 26 Figure 1.1 5 Degassing levels dropped in late 2009 ................................ ................................ ...... 26 Figure 1.16. By the end of 2009 activity had returned to similar levels to early 2006. ................. 26


v Figure 1.17. Simple model of advective diffusive heat flow across a vertical fracture at Comalito cinder cone. ................................ ................................ ................................ ....... 28 Figure 1.18. Model domain in the analytical solution. ................................ ................................ .. 30 Figure 2.1. Aerial photograph showing location of studied fumaroles adjacen t to Comalito cinder cone on flanks of Masaya volcano ................................ ......................... 38 Figure 2.2. (a) Temperatures recorded at different depths in a fumarole zone.. ........................... 41 Figure 2.3. (a) Spectrograms showing dominant frequencies in the fumarole soil temperature time series. ................................ ................................ ................................ .... 43 Figure 2.4. Time series of monitored data, with anomalous episodes highlig hted in gray ........... 44 Figure 3.1. (a) Topographic image of Masaya volcano, showing the active crater, Santiago ....... 52 Figure 3.2. (a) Topographic map of the central fumarole zone, showing offset extending to the NE, and locations of GPS and magnetic measurements (gray dots) and SP and CO 2 flux measurements (dashed lines 0 2). ................................ ............................... 56 Figure 3 3 Magnetic and SP profiles ................................ ................................ ............................. 58 Figure 3.4. Profiles of SP (solid green lines) and CO 2 flux (blue dashes) ................................ ..... 59 Figure 3.5 Example time series of fumarole temperature recorded at Comalito cinder cone ................................ ................................ ................................ ................................ ... 61 Figure 3.6. (a) SP measurements recorded along profile C .................. 62 Figure 3 7 Near surface temperatures as a function of permeability ................................ ............ 65 Figure 3.8. Fluid flow and heat transfer when there is a strong permeability contrast between the scoria and lava layers ................................ ................................ .................... 67 Figure 3.9. (a) Geometry used for numerical model, representing a cross section of the flank fracture at Masaya (Figure 4.1b) ................................ ................................ ............. 69 Figure 3.1 0. Conceptual model of the hydrothermal system at Masaya volcano, as determined by geophysical observations and numerical models ................................ ...... 73 Figure 4.1. a) Topographic image of Masaya volcano. ................................ ................................ .. 81 Figure 4.2. Flank fumarole temperature measurements ................................ ................................ 84 Figure 4.3. Details of fumarole temperature measurements ................................ .......................... 86


vi Figure 4 4 TOUGH2 model geometry ................................ ................................ .......................... 8 9 Figure 4 5 Changes in surface fumarole temperature are modeled with abrupt change s in mass flux and fluid enthalpy ................................ ................................ ............................. 9 2 Figure 4.6. Simulation of near surface temperature, modeled with abrupt changes in gas enthalpy and injection rate ................................ ................................ ................................ 94 Figure 4 7 Correlation between fumarole temperature and permeability ................................ ..... 9 6 Figure 4.8. Relationship between surface temperature and depth to the gas source ...................... 97 Figure 4 9 Heat of wetting curve according to Prunty and Bell (2005) ................................ ...... 100 Figure A 1. Early numerical simulations, of magma conduit and wall rock next to it. ............... 122 Figure B 1. Fumarole temperature data in 2006 ................................ ................................ .......... 124 Figure B 2 Fumarole temperature data in 2007 ................................ ................................ .......... 125 Figure B 3 Fumarole temperature data in 200 8 ................................ ................................ .......... 126 Figure B 4 Fumarole temperature data in 200 9 ................................ ................................ .......... 127 Figure C 1. Top: Spectrogram from thermocouple at 33 cm depth ................................ ............. 129 Figure C 2 Top: Spectrogram from thermocouple at 6 3 cm dep th ................................ ............. 1 30 Figure C 3 Top: Spectrogram from thermocouple at 95 cm depth ................................ ............. 131 Figure C 4 Top: Spectrogram from thermocouple at 150 cm depth ................................ ........... 1 3 2 Figure C 5 Top: Spectrogram from thermocouple 7 m away ................................ ..................... 133 Figure C 6 Top: Spectrogram of air temperature time series ................................ ..................... 129


vii Diffuse Degassing and the Hydrothermal System at Masaya volcano, Nicaragua Sophie C. P. Pearson ABSTRACT Hydrothermal systems change in response to volcanic activity, and in turn may be sensitiv e indicators of volcanic activity. Fumaroles are a surface manifestation of this interaction. We use time series of soil temperature data and numerical models of the hydrothermal system to investigate volcanic, hydrologic and geologic controls on this diff use degassing. Soil temperatures were measured in a low temperature fumarole field located 3.5 km from the summit of Masaya volcano, Nicaragua. They respond rapidly, on a time scale of minutes, to changes in volcanic activity also manifested at the summit vent. The soil temperature response is repetitive and complex, and is characterized by a broad frequency signal allowing it to be distinguished from meteorologic trends. Geophysical data reveal subsurface faults that affect the transport of fumarole gases. Numerical modeling shows that these relatively impermeable faults enhance flow through the footwall. On a larger scale, modeling suggests that uniform injection of fluid at depth causes groundwater convection in a permeable 3 4 km radial fracture zone tra nsecting the entire flank of the volcano. This focuses heat and fluid flux and can explain the three distinct fumarole zones located along the fracture.


viii We hypothesize that the rapid response of fumarole temperature to volcanic activity is due to increased flow of gas through the vadose zone, possibly caused by changes in the subsurface pressure distribution. Numerical models show that an abrupt injection of hot gas, at approximately 100 times background rates, can cause the rapid increase in temperature ob served at the fumaroles during volcanic activity. A decrease in hot fluid injection rate can explain the gradual decrease in temperature afterwards. Mixing with surrounding vadose zone fluids can result in the consistent and abrupt decreases in temperature to background level following hot gas injection. Fumaroles result from complex interaction of the volcanic hydrologic geologic systems, and can therefore provide insight into these systems. Increases in fumarole temperature correspond to increased gas fl ux related to changes in volcanic activity, suggesting that monitoring of distal fumaroles has potential as a volcano monitoring tool, and that fumarole temperatures can provide insight into the response of shallow gas systems to volcanic activity .


1 CHAPTER 1 INTRODUCTION A comprehensive understanding of volcanic systems is critical for reliable forecasting of eruptive activity, a major goal of volcanology (Sparks, 2003). Development of this understanding depends in part on understandi ng the relationship between groundwater and magma at active volcanoes in terms of heat and mass transfer. This understanding is needed because groundwater interaction on short timescales can increase the explosivity of eruptions (e.g., Mastin, 1991). On longer timescales, hydrothermal systems may cause significant alteration of volcan ic edifices, leading to weakening and collapse (e.g., Reid, M. E., 2004; Scott et al., 2005 ). Conversely, changes in the groundwater system may, in some circumstances, re sult from changes in volcanic activity ( Tanguy, 1994; Newhall et al., 2001; Shibata and Akita, 2001 ). At many volcanoes, groundwater magma interactio n is manifested at the surface by low temperature fumaroles and diffuse degassing, which occur where magm atic gases and water vapor reach the surface by porous flow over a relatively broad area. The spatial and temporal variations of low temperature fumaroles and diffuse degassing on volcanoes reflect the interaction between magmatic, hydrothermal, and meteor ological systems. These interactions, however, are complex. Through a series of works I address the physical controls on diffuse degassing on the flanks of Masaya volcano, an active volcano located in Nicaragu a. I explore the timescales of variation in flank fumarole


2 temperature and model these variations with the aid of geophysical measurements and spatially and temporally varying simulations of the hydrothermal system. Cumulatively, thes e studies provide insight into the causes of temperature variations in the flank fumaroles and processes that give rise to diffuse degassing at Masaya. Most importantly, these models seek to explain how fumaroles located some distan ce from the summit c rater at Masaya can respond rapidly, apparently on a time scale of minutes, to changes in the magmatic system. An extensive literature describes previous studies of fumaroles at active volcanoes and what they have sug gested about the primary controls on diffuse degassing (see Section 1.1 ). To be able to interpret fumarole temperature variations in light of volcanic activity the hydrologic system must be constrained, which a number of studies have attempted to do at Masaya volcano (see Section 1 .2). The volcanic activity must also be well describ ed (see Section 1.3) Numerical modeling of the hydrothermal system can then be attempted using the multicomponent, multiphase fluid flow simulation program TOU G H2 (Transport of unsaturated groundwater and heat) a code based on the physics of fluid flow as described in Section 1 .3. The relationship between measured fumarole temperatures and atmospheric and volcanic effects is explore d at Masaya volcano in Chapter 2 I find a correlation between fumarole temperatures and volcanic activity in the crater 3.5 km a way primarily identified from analysis in the frequency domain. However, rainfall also appears to have a sig nificant effect. I propose that fumarole temperatures increase due to elevated levels of gas flux from a pressure wave originating in the volcanic system. The rainfall may increase pore pressure to also enhance this effect.


3 A good conceptual model of the h ydro volcanic system is vital t o interpret fumarole temperature measurements I use geophysical constraints to create a hydrothermal model of the subsurface, including local and regional geolog y (see Chapter 3 ). These models show that fumaroles can result from steady convection of groundwater within the saturated zone, and that relatively impermeable faults can redirect gas and heat flow locally. During volcanic activity the fumarole temper ature signals gain distinct characteristics, which can be partially explained with numerical models. Geophysical constraints combined with information from a previous groundwater model allow a temporally varying model of the local hydrothermal system to be created (see Chapter 4 ). Model results shows that the fumarole temperatures c an be described by variations in gas flux, but could not be explained by any other single known mechanism This highlight s the u sefulness of fumarole temperatures as a potential volcano monitoring tool and as a proxy for total gas flux in some circumstances All three of the later chapters are written to stand alone; Chapter 2 was published in Geology in December 2008, Chap ter 3 is under review with Geophysics, Geochemistry, and Geosystems and Chapter 4 will be submitted in the near future. The term 'we' is used to refer to the collaborative nature of this project, particularly in terms of the geophysica l analysis of data in Chapter 3 1 .1 Previous fumarole studies Diffuse degassing has been studied for decades as a potential method for monitoring volcanoes. Fumarole temperatures, in particular, have the advantage of being


4 relative ly inexpensive to measure, and cost only a few thousand dollars to monitor continuously. N ot all active volcanoes have fumaroles (Varley and Armienta, 2001), but most do. Often, fumaroles are located in the crater of a volcano, making them danger ous to monitor I n some cases volcanoes also have fumaroles distal to the active crater (Baubron et al., 1991) Th ese allow easier and safer access particularly during elevated volcanic activity when deducing th e state of a volc ano is vital for minimizing damage to people and property. However, t he interpretation of the signals from these flank fumaroles is often problematic because their temperature and gas flow reflect a combination of both hydrol ogic and magmatic systems and are strongly affected by varia bility in the hydraulic properties of the rocks through which they pass. Previous studies help to understand controls on fumarole temperatures, thereby improving understanding of the m echanics of diffuse degassing and interpretation of its variations. A relatively early study of fumarole temperatures in the literature was done by Barquero (1988). He took regular measurements of fumarole temperature s at Volcan Poas, Costa Rica from 1980 to 1987. He detected increases in temperature prior to phreatic eruptions in both 1981 and 1987, although there were other small phreatic eruptions that occurred during a long term cooling trend. Another early study by Stoiber et al. (1975) detected a long term cooling in fumarole temperatures during waning activity following eruption s of Izalco Volcano, El Sal vador Neumann van Pedang (1963) looked at crater temperatures at a range of volcanoes in Indonesia but could not find any strong precursory signals related to volcanic activity. Sudo and Hurst (1998) found variations in temperature associated with volcanic activity at Aso volcano, Japan, during a study lasting 9 years. They measured soil


5 temperatures in two drill holes 200 m from the active crater, down to 150 and 70 m respectively. These measurements were not in a fumarole area, and temperatures va ried from 6C at the surface to 30C at depth. The large depths and high time resolution (every 10 minutes) allowed detection of atmospheric effects on temperature that had an increasing lag with greater depth, until other signals were detectable at 100 m and the atmospheric effects were minimal at 150 m. At that depth they observed the maximum temperature in November December 1989, just after the most vigorous volcanic activity from September to November 1989. A study by Yamashina and Matsushima ( 1999) at Mt Unzen found similar results to those of Sudo and Hurst (1998). There, they measured ground temperature in a horizontal cave outside of the fumarole areas. They measured temperatures every 5 10 minutes between 1991 and 1997. N o diurnal fluctuations in temperature were observed in the cave although there were seasonal changes. A lava extrusion reached its peak in June December 1991, and ground temperatures peaked in 1992, several months later. Connor et al. (1993a ; 1993b ) continuously m on itor ed high temperature (300 5 50C) fumaroles at Colima volcano, Mexico. They did two studies, first ly measuring fumarole temperatures at four locations at the summit of the volcano intermittently through 1990 (1993a) They observed formation of new fumaroles and daily temperature fluctuations of between 30C and 60C associated with barometric effect s, but otherwise the fumarol e temperatures remained consistently elevated. Modeling showed that the variation in wall rock temperature due to the fumaroles could explain demagnetization patterns observed in the area A negative magnetic anomaly at the same location as the fumaroles became the site of a new lava dome in 1991, which Connor et al. (1993a)


6 suggested was due to a preexisting fracture at depth which had previously been the site of steady state degassing and so formed an easier pathway for magma t o rise Another study by Connor et al. (1993b) m onitored high temperature fumaroles in five locations near the summit of Colima volcano between May 1991 and May 1992. The lack of obvious changes in volcanic activity allowed them to do an in depth study of the correlation between fumarole temperature and atmosphe ric variables, showing that rainfall only had an impact if it directly hit a thermocouple. Atmospheric pressure, however, had a strong correlation with fumarole temperature variations and clear diurnal variations in fumarole temperature were attributed to atmospheric forcing. There was also a cooling of the fumarole temperatures during winter. Connor et al. (1993b) then went on to create a steady state model of flow in a narrow fracture, showing that fumarole temperature will increase with increased mass f low, fracture width, and gas flow, and will decrease with greater distance to the magma body Other, more recent studies have found different correlations. For example, at Merapi volcano, Indonesia, Richter et al. (2004) measured strong dec reases in temperature immediately after rainfall, and found no significant correlation between atmospheric pressure and fumarole temperature. They recorded fumarole temperatures every minute in a high temperature fumarole zone 150 m from the active crater, and found that an increase in temperature of a few degrees often correlated with a certain type of seismicity, termed Seismic Cluster Ultra long Period Seismicity (SCULP). The fumarole temperatures had the same duration but varying magnitude responses to these events. Richter et al. suggested that both seismicity and temperature were due to a sudden release of pressurized gas from a crack 100 m below the active crater.


7 Zimmer and Erzinger (2003) looked in more detail at the hydrologic relationship with fu marole temperature at Merapi volcano. They found that with lower gas rates the ratio of magmatic gas to meteoric water decreased, causing cooler fumaroles. With a higher gas flux there was also a build up in pressure, resulting in relatively less meteoric water input and so greater fumarole temperatures. Geochemical measurements showed that immediately after rainfall the fumaroles had a higher concentration of water vapor, and also suggested that the input of meteoric and magma derived water ma y fluctuate regularly. Tedesco et al. (1991) also looked at hydrothermal inputs at Vulcano, Italy. After some recorded seismicity they noted an increase in fumarole temperature from 200 to 330C, as well as aerial extensions of the fumarole field and new fractures across the crater rim associated with increased gas flow. However, during their continuous monitoring study the seasonal variations dominated and there was no new seismic activity. From chemical variations they were able to deduce that the fumaro les contained a mixture of hot crater fluids and cooler hydrothermal fluids from residual magma degassing near the coast and that variations in fumarole temperature were therefore due to variations in mixing between the two. At Stromboli volcano De Gregorio et al. (2007) continuously monitored soil temperature in a low temperature fumarole on the southern rim of the crater, but only recorded measurements every hour. They measured temperatures anomalies in 70C and 55C fumaro les, mainly involving a decrease in temperature of up to 30C followed by an increase of a few degrees On one occasion however they observed a rapid and long term increase in temperature, and a loss of diurnal variations. This occurred a gain several


8 months later but with a smaller change in temperature They found that these anomalies corresponded well with measured variations in total dissolved gas pressure. As total dissolved gas pressure variations result from variati ons in the carbon dioxide present in the aquifer, De Gregorio et al. attributed the anomalies to changes in the degassing regime. They also noted a cooling effect of rainfall and a drop in temperature associated with volcano tectonic events, which they att ributed to degassing switching to the main conduit. In contrast, Chiodini et al. (2007) noted an increase in soil temperature surrounding a high temperature fumarole following seismic swarms at Vulcano, although this did not occur at the fumarole itself. T hey attributed this to an influx of hot fluids from a deeper part of the hydrothermal system, with a time lag for the fluids to reach the surface. Other studies have found similarly variable correlations between fumarole temperatures volcanic activity, se ismicity and/or surface features. Okada (1990) noted an increase in geothermal and fuming activity before an eruption at Mt Tokachi, Japan. Ingebritsen et al. (2001) found that fumarole temperatures were relatively stable around 90C at a variety of quiesc ent volcanoes in the western USA and were controlled primarily by snow melting. They were not directly related to mass flow, which did show some relationship with volcanic activity My own measurements at Masaya volcano, Nicaragua show a good correlation with volcanic activity (Pearson et al., 2008) but at Telica and Cerro Negro volca noes, Nicaragua the strongest correlation i s with soil moisture content, and therefore by extension with rainfall (Fig ure 1 1). Although there is considerable disparity between these systems, with study areas ranging from high temperature fumaroles in the crater to low temperature fumaroles far


9 from the crater to sites with no diffuse degassing, there is still commonality between the results that suggests that they can be useful for our study. In both high temperature fumaroles and non degassing areas, correlations with volcanic activity were found (Barquero, 1988; Sudo and Hurst, 1998; Yamashina and Matsushima, 1999) In both high and low temperature fumaroles, rainfall, seasonal affects, barometric pressu re and changes in gas flux were found to have an effect on temperature (Tedesco et al., 1991; Connor et al., 1993b; Richter et al., 1994; Zimmer and Erzinger, 2003; De Gregorio et al., 2007). Therefore it is necessary to understand the hydrothermal s ystem, both the thermal and hydrologic components, and to consider atmospheric e ffects on this system, in order to understand the response of fumarole temperatures to volcanic activity Figure 1 1 Fumarole temperature measurements recorded at a) Telica volcano and b) Cerro Negro volcano, Nicaragua in 2007. The strongest c orrelation is with soil water potential (SWP), a measure of soil moisture content. No known volcanic activity occurred at eit her volcano during this measurement period. There were considerable problems with instrumentation, which may explain the difference in magnitude between the SWP measurements at the two volcanoes.


10 There are also several other ways to monitor fumaroles that can provide insight into the hydrothermal system and help us to understand fumarole temperature variations. Self potential (SP) is one such method, as it provides a way to track changes in movement of hydrothermal and meteoric fluid. For example, Hashimoto and Tanaka (1995) noted an increase in SP before a lava e xtrusion at Mt Unzen, and a change in the SP signal afterward which they attributed to changes in hydrothermal convection Friedel et al. (2004) measured noticeable changes in SP that correlated to rainfall at Mt Merapi; there was no volcanic activity duri ng their study. Zlotnicki et al. (2008) measured changes in SP, total magnetic field and ground temperature on the flank of Taal volcano, Philippines, associated with a seismic crisis which prompted self evacuation of hundreds of people. Although this was not followed by a volcanic eruption, it suggested shallow magma and a possible new site of activity along a fissure in the northern part of the volcano. At Masaya volcano, a consistent positive SP anomaly was attributed to upwelling of hydrothermal fluids ( Lewicki et al., 2003; Lehto, 2007). However, there was no change in the anomaly during observed changes in volcanic activity Instead, the SP signal was dominated by atmospheric effects, with diurnal signals that correlated with barometric pressure. Dur ing rainfall the SP decreased and took a minimum of several hours to return to average levels. The lack of a correlation between volcanic activity and SP could be due to the suppressing effect of rainfall, or because any new fluids have the same io nic charge as surrounding fluids (Lehto et al., 2007). SP at Masaya shows that hydrothermal fluids are constantly circulating, but that atmospheric effects must be included when studying diffuse degassing at this volcano.


11 CO 2 flux i s another more expensive, method to track changes in fumaroles, and is often done concurrently with soil temperatures ( e.g. Chiodini et al., 1998; Rogie et al., 2001 ; Aiuppa et al., 200 4 ) For example, at the Phlegrean fields, Italy soil CO 2 and temperature were measured once or twice per month (Granieri et al, 2003). Some correlation was observed between soil CO 2 flux and ground deformation, but environmental effects dominated. Approximately 58% of CO 2 flux variance could be explai ned by soil humidity, and therefore there was also a significant relationship with rainfall. They attributed this to the fact that the monitoring station was on high ground, and so rainfall flowed into the surrounding area and diverted gas flow toward the CO 2 monitoring station. At Masaya volcano previous studies have shown a good spatial correlation between elevated SP, CO 2 and temperature (Lewicki et al., 2003; Lehto, 2007) The subsurface geology is therefore likely to be an important factor in controlli ng subsurface fluid flux Continuous CO 2 flux measurements collected at a fum arole zone near Comalito cinder cone 3.5 km from the active crater varied between 300 and 11000 g/m 2 d (Figure 1.2) This is well within the range measured by Lewicki et al (2003). There were no clear long term trends, other than the increase between 2007 and 2008 which may due to instrument error as the fluxmeter stopped working in August 2007 and/or possibly a change from the rainy season to the dry season There is a diurnal signa l which suggests atmospheric forcing, and the consistently elevated flux show s that the hydrothermal system is constantly active.


12 Figure 1 2 CO 2 flux measured at Comalito cinder cone during this study It is constantly elevated, with predominantly d iurnal va riations The large increase in temperature between 2007 and 2008 may be due to instrument error or the change from the rainy season to the dry season. Todesco et al. (2003) did a study of gas fluxes at the P hlegrean Fields, focusing on numerical simulations. They modeled the shallowest 1500 m of the system using geophysical, geochemical, petrological and well log data as constraints. With a fully saturated model they discovered that injection of hot water and CO 2 could reproduce the distribution of heat obser ved at the surface at the Phlegrean Fields and also the single phase gas region at the crater inferred from geochemistry. The results were dependent on the temperature of injected gas, the rate, its sourc e depth, and host rock properties. High permeability or a deeper source resulted in distal areas remaining cold and no dry gas zone. Increasing permeability at shallow depths caused a decrease in pressure and temperature in the single phase gas region. Enhanced permeability near the single phase region, however, changed the pressure and temperature conditions to enhance boiling and change both gas composition and flux rate. Even with injection of hot fluids from a steady sou rce, the surface discharge composition and rate could be highly variable both in time and space and was more limited by rock properties. Changes in the source fluid took longer to be recognized than changes in permeability (on the order of years to decades )


13 although periods of intense magma degassing could still induce important changes in the reservoir that fed the surface fumaroles (Todesco et al., 2004). To understand the impact of variations in the hydrothermal system, Todesco et al. (2004) modeled defo rmation associated with it. Using the temperature and pressure from the numerical model as input in a FLAC3D model, Todesco et al. found that periods with higher injection rates caused heating and build up of pore pressure that could trigger significant rock deformation. The models suggested that periods of intense degassing resulted in fast uplift followed by slower subsidence and large changes in gas composition. The deformation would be focused at the end of these periods while the gas c omposition would change several months after the injection rate was reduced again, as observed at the Phlegrean Fields. The hydrothermal model of the Phlegrean Fields (Todesco et al., 2003) was also compared to time varying gravity signals to further v alidate the model. The distribution of fluid density is affected by ascent of hot fluids, phase changes and variable gas composition (Todesco and Berrino, 2005). Although the gravity signals associated with this are small compared to injection of new magma they are also very shallow so may still be detected. Residual gravity signals increase when the fluid injection rate increases because mass injection increases, water vapor condenses due to increased pore pressure, and liquid water is displaced upwards. Gravity changes computed as a function of volume, depth, distance from the observation point and average fluid density were summed over the numerical model from Todesco et al. (2003) who showed that by varying the injection rate, composition and durat ion of volcanic crises, both the distinct gravity trend and variable gas composition recorded at the Phlegrean Fields could be


14 recreated. The volcanic crises were attributed to a pulsating source of hot fluids which discharged large amounts of CO 2 rich flu ids, also affecting hydrothermal circulation. However, in 2005 the rate of subsidence declined and new uplift began, with gradual gas enrichment over 5 years and a significant drop in the gravity residual. This could not be recreated with the previous mode l (Todesco et al., 2006). Todesco et al. suggested that the mechanism feeding the fumaroles had probably changed, as had the magma composition and subsurface permeability distribution. A wider source region over the previous few years and more frequent and shorter unrest crises suggested that the rate of fluid injection, timing and duration needed to be reassessed to explain more recent signals, and that variable magmatic degassing and emplacement of deep source fluids was probably occurring (Todes co et al., 2006). Although the system at the Phlegrean Fields is not analogous to the one at Masaya, the numeri cal models are useful for our study. A t Masaya volcano a lava lake occ asionally forms in the crater, showing that the system is open ( Rymer et al., 1998; Williams Jones et al., 200 3; Stix 2007) unlike that at the Phlegrean Fields This makes it unlikely that increased injection would result in the same degree of deformation. The fumaroles a t the Phlegrean Fields are also much hotter, but the dominant controls are still likely to be similar Therefore that we need to particularly consider gas injection rate, rock permeability and to a lesser extent source depth and gas injection temperature in our models Previous studies suggest that the relationship between diffuse degassing and fumarole temperature is complex, but that variations in the temperature are often attributable to changes in gas flux through the system. Despite being unable to ex plain


15 and/or recreate all of the signals observed at fumaroles on active volcanoes, these studies have provided important insight into factors controlling diffuse degassing, and its link with volcanic activity. 1 .2 The h ydrologic system at Masaya The hydrologic system at the Phlegrean Fields is well constrained, and modeling has shown the importance of incorporating it into models of diffuse degassing. However, t he hydrologic system at many active volcanoes is poorly known (Hurwitz et al., 2003) and hydrologic properties are rarely characterized. At Masaya volcano most of the numerous studies have focused on the geological structure (e g. Williams, 1983 a ; McBirney, 1956; Rymer et al., 1998) or the gas composition (e g. Duffell et al., 2003; Lew icki et al., 2003), which are described in the relevant chapters of this dissertation. However, the hydrologic system has been previously assessed and will be described here. In 1993, Instituto Nicaraguense de Acueductos y Alcantarillados an d Japan International Cooperation Agency used well soundings outside of the caldera at Masaya volcano to infer a water table elevation of 190 m eters a bove s ea l evel (masl) to the south of the caldera and 130 masl to the north of the caldera. They ass umed fairly consistent gradients across the caldera. Transient electromagnetic soundings (TEM) showed that this was an oversimplification and that a more detailed study was necessary to approximate groundwater flow within Masaya caldera (MacNeil et al., 20 07). TEM soundings at 30 sites throughout Masaya caldera detected a highly resistive layer overlying one or more conductive layers (MacNeil et al., 2007). This was interpreted as dry basaltic rocks overlying porous, saturated lava flows. Therefore the


16 boun dary between the two layers corresponded to the water table. A simple two layer model suggested that the water table was a subdued reflection of topography (Fig ure 1.3 ), although the near vent area was possibly vapor dominated and was too complex for the caldera scale model to recreate. One important feature of the geology that was included in the hydrologic model is that c aldera bounding faults form walls that are up to 400 m in height and have bounded lava flows in the past (Williams, 198 3 b ). Dikes have possibly been injected along these faults (Williams, 1983b). This has created a strong lithologic and hydrologic boundary that minimizes flow of groundwater between the caldera and the surroundings, causing the caldera to be relatively isol ated hydrologically Therefore evapotranspiration, degassing and rainfall can all be assumed to be contained within the caldera. MacNeil et al. (2007) used the TEM soundings to create a 3D groundwater flow model with MODFLOW2000. The average measured rain fall of 1.6 m/year, degassing rate of 400 kg/s (Burton et al., 2000), and estimated evapotranspiration rate of 1.1 m/yr were used as constraints in the model. They found that, as with the two layer model, the water table was at its shallowest depth close t o Laguna Masaya, which has the lowest elevation in the caldera (Fig ure 1.4 a). Close to the active vent, which has the highest elevation, dramatic gradients occurred and the water table followed topographic gradients The modeling study concluded that a hig h permeability zone beneath the volcanic edifice must exist at depth (Fig ure 1.4 b) that transmit s groundwater toward the active vent.


17 b a Figure 1 3 a) Topographic map of Masaya volcano. b) Elevation of the water table as estimated from TEM soundings. Reproduced from MacNe il et al., 2007 (with permssion from Elsevier).


18 Figure 1 4 Groundwater distribution at Masaya in meters above sea level accor ding to MODFLOW2000 results. Arrows indicate direction of groundwater flow. a) At the surface. b) At depth, where a high permeability zone results in groundwater flowing toward the active crater (Reproduced from MacNeil et al., 2007 with permission from El sevier). A study of self potential wavelet tomography by Mauri et al. (2010) is in good agreement with the model of MacNeil et al. (2007), although their estimated depths to the water table are consistently less (Fig ure 1.5 ). This supports the model of a water table that approximately follows topographical contours, with a lower hydraulic head further


19 from the vent (Figure 1.5a) but with groundwater flowing toward the active crater at depth (Figure 1.5b) Figure 1 5 Estimated water table depths. Brown line represents land surface, points represent water table estimated from SP wavelet tomography, and gr ay profile represents water table estimated by MacNeil et al. (2007) (Reproduced from Mauri et al., 2010 with permission from Elsevier). 1 .3 Volcanic a ctivity at Masaya volcano ( adapted from Tenorio et al., 2010) Although the volcano hydrologic system at Masaya volcano is relatively stable and the volcano has been consistentl y degassing for over a decade both the magnitude and the style of degassing is constantly fluctuating The last major explosion was in 2001, when a number of people were injured by ash and rocks that were ejected f rom a new vent that formed in the crater. From 2001 until 2005 occasional ash explosions and incandescence were observed In 2005 a landslide blocked the active vent and degassing dropped to low levels. From 2006 until 2009, my study period, activity varied between low level degassing, ash explos ions, incandescence and the formation of a new lava lake (Figure 1.6) This provid es an excellent opportunity to observe any changes in fumarole temperature related to changes in volcanic activity


20 Figure 1 6 Timeline of crater activity at Masaya volcano during the study period. At the beginning of the study period, f rom October 2005 until May 2006 gas emanations from the crater were strong and consiste nt. On May 25 2006 at 11:30 am the Volcan Masaya National Park staff reported a rockfall in the crater and on May 30 incandescence was visible in the active vent The park staff some incandescen ce in the old crater ve nt There was no corresponding change in seismicity or fumarole temperature (Figure 1.7 1.8 ). In June 2006 the crater although only at night. The incandescence was visibl e at night for the next year, until May 2007.


21 Figure 1 7 Real time Seismic Amplitude (RSAM) and total number of earthquakes recorded per month during the study period. There was no clear correlation with volcanic activity. Figure 1 8 Fumarole temperature recorded about once a month by INETER scientists. Th ere are no obvious correlations with volcanic activity. Crater degassing remained consistent from June until October 2006, when surface features changed in the crater Mini DOAS measurements and visual observations showed a marked increase in both degassi ng and magmatic activity in October (Figure


22 1. 9 ) A fracture formed in the crater floor during the first two weeks of October, and National Park staff noted u nusual incandescence on October 21. On October 22, National Park staff observed that a new vent h ad formed in the North East part of the crater, with an estimated diameter of more than 5 m. Strong degassing and l ight w ere visible from a new lava lake at t he ven t (Figure 1. 10 ) The vent formed from the larger explosion in 2001 still hosted more degassi ng than this new vent however. Intense rain from 20 22 October hindered observations of the new vent formation, and also resulted in acid rain that damaged vegetation and trees to the east of the crater. Figure 1 9 Sulfur dioxide flux from the active crater during the study period. There was a marked increase corresponding to the formation of the lava lake in October 2006 SO 2 flux (t/d)


23 Figure 1 10 Image s of the new lava lake that developed in the crater in October 2006. For the next 7 months gas emissions from the new crater remained strong and incandescence was still visible at night. Seismicity changed slightly, with an increase in RSAM units (Fig ure 1.7 ) and a change in dominant frequency from 1.5 Hz to 5 Hz in February 2007. In May 2007 emissions decreased slightly, enough that the temperature in the new vent could be measured ( Figure 1.1 1 1.12 ). It was discovered that neither the new mini DOAS installation nor the broadband seismic station were working. Figure 1 11 On average temperatures decreased slowly at the 2006 vent, even after a rockfall blocked it in early 2008


24 Figure 1 12 Degassing from the vent f ormed in 2006. Lava was no longer visible by May 2007 From July 2007 activity changed as the plume became more ash rich and reached higher eleva tions, resulting in a number of reports from the Washington VAAC ( Global Volcanism Program, 2010 ) In July 200 7 seismicity, fumarole temperature and gas emissions remained consistent with the month before (Figures 1.7 1.9) However, on 15 July National Park staff reported ash deposits on both the road and the lookout station to the southeast of the crater. On 16 J uly they observed that the plume was a darker color than the usual white, water vapor dominated plume INETER scientists attributed this to either an increase in ash or sulfur content. A large, ash rich plume was visible from July 2007 until May 2008 drif ting in a variety of directions. The degassing was primarily from the 2001 vent, particularly after a rockfall closed the 2006 vent in January 2008. In September there were some rockfalls in the crater and in October acid rain affected vegetation and new fumarole activity was observed in the east part of the crater (Figure 1.13)


25 Figure 1 13 New fumaroles formed to the east of the active cra ter in October 2007 On 18 June 2008 2 ash explosions were felt seismically, and emitted medium amounts of ash and gas to 1 km. There were no visible rockfalls or ash deposits afterwards. From June 2008 until July 2009 the volcano was the site of some low level tremor and persistent, strong degassing although w ith reduced ash content (Figure 1.14) In August 2009 the crater gas emissions visibly decreased, and activity remained low (Figure 1.15) until late November 2009, when abundant gas emissions were visible from th e SE of the crater (Figure 1.1 6 )


26 Figure 1 14 Strong, consistent degassing in Figure 1 15 Degassing levels dropped in September 2008. late 2009. Figure 1 16 By the end of 2009 activity had return ed to similar levels to early 2006.


27 1.4 Numerical modeling Improvements in computational power have revolutionized volcanology, allowing increasingly sophisticated numerical simulations to approximate a volcanic system or certain aspects of it. We bega n by creating simple models of heat and mass transfer using a finite difference, implicit code we w r ote in MATLAB (see Appendix) and with the COMSOL M ultiphysics program a finite element code to solve partial differential equations. W e then used COMSOL multiphysics to create a 2 dimensional heat flow model of a fracture and its adjacent rock at Comalito cinder cone (Figure 1.6) which was observed to be the site of considerable degassin g. we used the 2D convection and conduction mode of COMSOL The interior of the model was only subject to heat where q is the heat flux, k is the thermal conductivity, and is the local temperature gradient. The top and bottom boundaries were at fixed temperature to represent t he heat source and the surface, respectively. The right boundary was no flo w. The left boundary of the simulation was given by forced convect ion to simulate heat flow through the fracture : where q o is the heat flux across the boundary into the model (in this case 0) h is the heat transfer coefficient and T inf is the temperature outside of the model This resulted in a simple conductive profile, with temperature decreasing moving through the wall rock away from the convective boundary (F igure 1. 17 )


28 Figure 1 17 Simple model of advective diffusive heat flow across a vertical fracture at Comalito cinder cone. The left, convective boundary simulates the fra cture. To create more sophisticated models we used the TOUGH2 numerical code. This is a numerical simulator of multicomponent, multiphase fluids in porous rock (Pruess, 1991). It can simulate high temperatures and pressures, f low of gases and liquids, and flow through porous and fractured media, making it ideal for models of heat and mass transfer around and within a fracture on the flank of an active volcano. Originally developed by Lawrence Berkeley National Laboratory as a r esearch code called MULKOM, TOUGH (Transport Of Unsaturated Groundwater and Heat) was first released to the public in 1987 (Pruess, 2004). It was designed for geothermal reservoir simulation primarily, particularly for the Yucca Mountain project, and took into Fracture


29 account the nonisothermal nature of flow, the importance of both boiling and condensation, and the highly nonlinear nature of two phase, water steam flow (Pruess, 2004). It was gradually expanded to include more applications and modules (TOUGH2), non aqueous phase liquids in contamination problems (T2VOC and TMVOC), coupling with geochemical reactions (TOUGHREACT) and coupling with rock mechanics (TOUGH FLAC ; Pruess, 2004). My study focuses on flow of groundwater and gases through porous rock and their affect on surface temperature and gas distributions, and therefore we worked exclusively with TOUGH2. 1 4 .1 Governing equations (adapted from Pruess et al., 1999; Pruess, 2004) TOUGH2 solves energy and mass balance equations for flow of multicomponen t, multiphase fluids through porous and permeable material. These can be written as: mass accumulation = mass flux + sinks and sources or; [1] where k is the co mponent 1... nk (air, water etc.), V n is an arbitrary subdomain of the flow system, n is the closed surface bounding it, and n is a no rmal vector on surface element n pointing inward into V n (Fig ure 1 18 ) Essentially, the rate of change of fluid mass ( M ) in V n is equal to the net inflow across V n ( F ) plus the net flux from other sources ( q )


30 Figure 1 18 Mod el domain in the analytical solution. Looking at the mass accumulation term first, M k is the total mass of component k obtained by summing over fluid phases : [2] where is the porosity, S is the saturation of the phase (the fraction of the pore volume occupied by that phase), is the density of the phase, and X is the mass fraction of component k presen t in phase Similarly, the heat accumulation term of a multiphase system is given by: [3] where r is the density of the rock, c r is its specific heat capacity, T is temperature and u is the specific internal energy in phase In this case heat accumulation is not dependent on the component, and so there is only one term per cell. Next the advective mass flux term is calculated by summing over the phases: [4]


31 where F is the multiphase version of Darcy's Law: [5] Here, u is the Darcy velocity, is absolute permeability, is relative permeability of the phase, is the visc osity of the phase and g is gravity. P is the fluid pressure, given by the sum of the pressure of the reference phase (generally a gas) and the capillary pressure. The heat flux has a similar form and is calculated to include both conduction and conv ection: [6] where is the thermal conductivity and h is the specific enthalpy in each phase The diffusive fluxes of all phases are also included by summing the fluxes of each component k in phase [7] where is the tortuosity (the former depends on the porous medium, the latter on the phase saturation), d is the diffusion coefficient and X is the mass fraction as above. Other complexities can also be included that are not relevant to this study and so are not described, for example hydrodynamic dispersion and the gas phase permeability increase that occurs at low pressures. To transform these equations into a partial differential equation from which a numerical solution is generally derived, Gauss' theorem can be applied to equation [1]: [8]


32 1 4 .2 Solving the equations To solve the equations above, a numerical approximation is made. The spatially and temporally continuous integrals that make up the mass accumulation and mass flux terms are discretized. This is done by averaging for each cell and then summing all cells, the mass accumulation over volume and the mass flux over area. To discretize the mass flux term, the change in pressure between nodal grid blocks is divided by the dist ance between the blocks to give the P term in Darcy's Law ( E quation 5). The gravity is the net value in the appropriate direction between the two blocks, and all other values are averaged over adjoining cells. This is then substituted back into equati on 1 to give a first order differential equation. The flux and source and sink terms are calculated for each new time step to give a fully implicit, first order backward finite difference scheme. This results in a series of coupled, non linear algebraic eq uations in each cell, with one for each mass component and one extra for the energy balance. These are solved simultaneously and then Equation 1 is rearranged to give a residual on the left and all the original terms on the right. Newton R aphson iteration is used to get the residual below a preset convergence tolerance, with automatic time step adjustment. There is 100% upstream weighting of flux terms at interfaces for unconditional stability, and time steps that are not too smal l Solving these equations is dependent on the fluid properties (pressure, volume, temperature), and the interaction between the fluids and permeable medium. TOUGH2 includes a number of equations of state (EOS) that describe these. We used EOS3, air and wa ter, for all of our modeling as this allows us to model the vadose zone. We therefore assume that air is representative of the volcanic gases, although at Comalito cinder cone


33 the fumarole gases were measured to be around 50% air (Chiodini et al., 2005). H owever, compared to the uncertainties in our models and the relatively small range of temperatures and pressures of the fluid, this is unlikely to be a major source of error. A number of assumptions are included in EOS3, including the fact that air is an i deal gas, air and vapor partial pressures are additive in the gas phase, and that the solubility of air in water is a constant, although the latter's small magnitude and slowly varying dependence on temperature make this unlikely to be a large source of er ror. Complications arise from the multiphase nature of TOUGH2 that have to be addressed. For a single phase model, temperature and pressure are the inputs. However, for a multiphase model pressure is not independent of temperature and so another variable m ust be used in place of one of these. A particular quirk of EOS3 is the fact that air mass fraction is input with pressure and temperature for the single phase case, while gas saturation is the input for the two phase case. As both of these fall between 0 and 1, the input representing the gas saturation is actually S g + 10 so that the single and two phase simulations are numerically distinguishable. To determine if the model changes phase, the pressure is the determining factor. For a single phase liquid m odel, the saturated pressure must be less than the fluid pressure. For a two component (air + water) model, the sum of the saturated vapor pressure and the air pressure must be less than the fluid pressure otherwise a gas phase will evolve. The models can switch between the two, with the input of air mass fraction being replaced by very small gas saturation. To actually implement these models, the permeability, porosity, capillary pressure relations thermophysical properties of the fluid, initial an d boundary conditions of the flow system and sinks and sources must be estimated and assigned These are user inputs,


34 except for most of the thermophysical properties of the fluid, which are calculated within TOUGH2. The model geometry, program o ptions, computation parameters and time stepping information are also specified by the user. In this project we used fairly simple, regular grids and therefore were able to use the Petrasim interface to TOUGH2, a more graphical approach. This allows the us er to specify values in tables and on images of the model domain rather than using the text file format utilized by TOUGH2. It also outputs images of the model results and time cell history plots. Probably the most critical part of TOUGH2 modeling i s determining the boundary conditions. For Neumann conditions the fixed heat and/or mass flux is simply prescribed as a source or sink. This can vary with time. For no flux conditions the connections across the boundary are simply removed. Dirichlet condit ions, with fixed temperature and/or pressure, are more difficult to implement. For a fixed Dirichlet condition which does not vary with time, a cell can be set as inactive in TOUGH2 (fixed state in Petrasim). Inactive cells appear in flow connections and i nitial conditions but are ignored otherwise. Another way to set a Dirichlet condition is to specify very large volumes for the boundary cells so that the thermodynamic conditions do not change from fluid or heat exchange. For time dependent Dirichlet co nditions, a sink or source has to be introduced into the cell at the correct rate to get the desired variations. For a simple temperature boundary the porosity and permeability are set to negligible values to limit mass flow into and out of the cell, and t he heat flow is specified to give the desired temperature. A similar approach is taken for fixed pressure conditions, but with zero thermal conductivity and mass flow into or out of the cell specified instead.


35 For more information on code architecture, im plementation and applications, see Pruess et al. (1999), Pruess (2004) and the Petrasim users manual. 1 4 .3 TOUGH2 validation (adapted from Oldenburg and Pruess, 1994) To assess the validity of the numerical approximations created by TOUGH2, a number of comparisons were made with models established in t he literature. Of particular relevance for the advective diffus ive part of TOUGH2 as utilized in our model s the seawater intrusion problem of Henry (1964) and the free convection problem of Elder (1967) we re modeled and compared with previous results. In the Henry problem, injection of fresh water into the left hand side of a two dimensional model results in constant flow of fresh water to ward the right where the boundary is set to be that of seawater and the pressure is set at hydrostatic pressure This results in density variations of 2.5 %, with seawater along the right boundary and the base of the model near this boundary and fr esh water throughout the rest of the model domain. Comparisons with results from Voss and Souza (1987) which are representative of the literature show good agreement Results are found to be relatively insensitive to grid dimensions. In the Elder problem heat is inject ed from below resulting in free convection. A common adaptation of this problem is for solutal convection, where there is a salt source at the top rather than a heat source at the b ase Oldenburg and Pruess (1994) simulated the solutal c onvection problem with TOUGH2 They found similar results to those calculated by Elder (1967) and Voss and Souza (1987) for similar grid dimensions However, Elder (1967) also did a laboratory experiment with a Hele Shaw cell and


36 observed a number of convection cells forming compared to the two of the numerical sim ulations. The center of the apparatus was also the site of upwelling fluid compared to downwelling fluid in the numerical simulations. By reducing the grid spacing in the TOUGH2 simulations, much better approximations of the laborato ry results co uld be computed. C omparisons of TOUGH2 with previous experimental and numerical simulations show that it can adequately represent the systems we are trying to describe. Grid dimensions are found to be important in the Elder problem and therefore in othe r models involving free convection For each of my models in this study, I tested a var iety of grid resolutions to ensure that the heat and mass flow was not being aliased.


37 CHAPTER 2 RAPID RESPONSE OF A HYDROLOGIC SYSTEM TO VOLCANIC ACTIVITY 2 .1 Introduction Time series of fumarole soil temperature variations are of interest both in monitoring active volcanoes ( De Gregorio et al., 2007; Yamashina and Matsushima, 1999 ), and for studying the nature of heat and mass transfer in volcanic systems ( Tedesco et al., 1991; Connor et al., 1993b; Sudo and Hurst, 1998 ). H igh temperature fumaroles located in the craters of vol canoes have been observed to increase before eruptions (Barquero, 1988 ) T emperature responses of lower temperature flank fumaroles are often less clearly correlated with changes in volcanic activity ( Richter et al., 2004; Sudo and Hurst, 1998 ; Yamashina and Matsushima, 1999) in part due to interaction between volcanic and hydrologic systems. Yet this interaction can play an important role in the nature and timing of volcanic eruptions (e.g., Nakada et al., 2005). Here we present time series showing very rapid and significant temperature changes in low temperature fumaroles in response to volcanic activity during two episodes at Masaya volcano, Nicaragua. Rapid and multi channel data acquisition allows us to identi fy cyclic variations in temperature in both of these volcanic episodes, revealing new details about the response of the hydrologic system to changes in volcanic activity.


38 2 .2 Masaya Volcano Masaya volcano, Nicaragua (Figure 2 .1), is one of the most p ersistently active volcanoes in central America (McBirney, 1956; Williams, 1983 a ; Stoiber et al., 1986). It has been constantly degassing from Santiago crater (Figure 2 .1) since 1993 ( Duffell et al., 2003; Mather et al., 2006) and has one of the largest r eported noneruptive gas fluxes at any volcano, with intermittent lava lake development associated with an extensive magma plumbing system (Stoiber et al., 1986; Rymer et al., 1998; Williams Jones et al., 2003). Low temperature (~65C) fumaroles zones are l ocated along a northeast trending fracture system that extends down the edifice of the volcano for ~4 km (Figure 2 .1) Near Comalito cinder cone, located along the fracture 3.5 km from and 200 m below Santiago crater, some of the highest known diffuse car bon dioxide fluxes from low temperature fumaroles worldwide have been identified and mapped. The presence of magmatic gases in the fracture system ( Lewicki et al., 2003; Shaw et al., 2003 ) increases the likelihood tha t Comalito fumaroles might respond to change s in volcanic activity Figure 2 1 Aerial photograph showin g location of studied fumaroles adjacent to Comalito cinder cone on flanks of Masaya volcano (Instituto Nicaraguense de Estudios Territoriales (INETER); W. Strauch, 2006, personal commun.). Dashed line represents fracture zone inferred from the distributio n of thermal areas measured by Lewicki et al. (2003). Active degassing is occurring in Santiago crater, where vent widened noticeably in June 2006 and small lava extrusion was observed from new vent onto crater floor on 23 October 2006. Location of Masaya is shown in inset map.


39 Throughout the 9 month s tudy period, seismic tremors were recorded and incandescence was visible at night in Santiago crater. T wo periods of increased volcanic activity were noted by Instituto Nicaragense de Estudios Territoriales (INETER), in June and October 2006 (W. Strauch, 2006, personal commun ) In June vent widening with increased incandescence follow ed a rockfall on 25 May 2006 During 21 23 October 2006 a small lava lake formed in the crater following a peak in seismicity and an increase in SO 2 flux (W. Str auch, 2006, personal commun ). 2 .3 Method Monitoring equipment was installed in May 2006 near Comalito cinder cone (Fig ure 2 1). Four tubes were driven into the ground within a 1 m 2 area in a fumarole zone, to depths of 33 cm, 63 cm, 95 cm and 150 cm and a fifth was driven to 150 cm in a different fumarole area, located 7 m away. An ungrounded chromel alumel ( t ype K) thermocouple with a resolution of 0.1 C was inserted in each tube, which was then sealed with putty. Rainfall, barometric pressure, an d atmospheric temperature were also recorded adjacent to the thermocouples. Each sensor was sampled every ten seconds. Temperatures and atmospheric pressure were averaged every five minutes and rain measurements were cumulative over each five minute period Real time seismic amplitude measurements (RSAM) were calculated by INETER from seismic data recorded at a broadband seismometer near Santiago crater (Fig ure 2 1). 2 .4 Results Fumarole soil temperature time series for the sampl ing period are shown in


40 Figure 2 .2 a The coolest temperatures were consistently recorded by the deepest thermocouples, and the highest temperatures were recorded by the shallowest thermocouples (Fig ure 2 .2 a ). This is due to the presence of a thin (~5 20 cm) low permeability cl ay rich layer at the surface, causing accumulation and convection of hot gases in the shallow subsurface. Distinctive signals occur red during the two periods of anomalous crater activity in June and October. During these periods temperature increase d to a peak value over ~20 min utes followed by an exponential decrease asymptotically approaching a threshold temperature for as long as one day, and then a sudden drop to background temperature occurred (Figures 2.2 b d ). In individual fumaroles, the threshold temperatures were remarkably consistent across the numerous heating cycles. For example, the shallowest thermocouple (33 cm) recorded a background temperature of 72 C, a peak temperature of ~ 80 C and a threshold temperature of ~74 C. Analysis of the frequency spectrum of the soil temperature time series reveals unique characteristics during volcanic activity. At background state the time series were dominated by diurnal and semidiurnal variations. These periodicities disappeared during the two periods of volcanic activity, when series included a broad range of frequencies (Figures 2 3a and 2 .3 b). Thus, while the magnitudes of these anomalies were less than 5 C and their signatures were complex, the anomalous behaviors were easily identified on spect rograms.


41 F igure 2 2 ( a) Temperatures recorded at different depths in a fumarole zone. Note the increases in temperature in early June and late October, highlighted in gray, corresponding to periods of increased volcanic activity. (b ) Heating cycles recorded in the temperature signal during the first period of volcanic activity in June 2006. In several instances, sharp increases in temperature appear to correspond to rainfall events shown along the bottom of the graph. (c) The distinc tive signal formed by a sharp increase in temperature followed by an exponential decrease and then a sudden drop recorded at 33 cm depth. (d) One of these cycles recorded at 95 cm depth consists of (i) sudden increase in temperature due to an increase in p ressure; (ii) e xponential decrease in temperature while the system is stabilizing during higher gas flow; and (iii) rapid decrease in temperature to background (see Section 2.5)


42 I ndividual rainfall events often coincided with rapid increases in the soil temperatures of as much as 3 C. Of these temperature spikes, ~ 30% occurred during rainfall, although it was raining during only 1.35% of the sampling period (Fig ure 2 4d); thus te mperature spikes have a positive correlation with rainfall with >99% confidence (binomial interval estimate). During volcanic activity the sudden increase in soil temperature at the onset of cycles was often accompanied by rainfall (Fig ure 2 2b). In non eruptive periods, soil temperatures were characterized by a slight increase in average temperature (0.0010 C/day 0.0094 C/day), diurnal and semi diurnal variations (inverse to barometric pressure), and seasonal variations associated with rainfall (Fig u re s 2 2a, 2 4a, and 2 4 b). There was a peak in RSAM during the temperature excursion prior to the formation of the new lava lake in October 2006, but many larger RSAM spikes were observed when there were no changes in temperature or reported volcanic activity (Fig ure 2. 4c)


43 Figure 2 3 (a) Spectrograms showing dominant frequencies in the fumarole soil temperature time series. Each spectrogr am corresponds to a different thermocouple buried to the depth indicated. During the anomalies in June and October 2006, marked by red arrows, the fumarole temperature signals contained a much greater range of frequencies than during background activity, r evealed by greater intensity of high frequencies during these periods. Changes in frequency content in late August and mid November occurred when the datalogger was opened for maintenance. (b) Fourier transforms of the temperature record from 95 cm depth d uring each of the anomalies, in June and October 2006, and during background state, in December 2006. The dotted red lines highlight diurnal and semi diurnal frequency variations, which are only observed at the thermocouples during background conditions.


44 Figure 2 4 Time series of monitored data, with anomalous episodes highlighted in gray. (a) Average fumarole temperature and daily temperature fluctuations increased from the rainy season to the dry season. (b) Ground fumarole temper ature was recorded at 95cm depth but is representative of changes at all depths. (c) Peaks in RSAM sometimes corresponded to increases in temperature, particularly at the onset of the June and October anomalies. (d) Differential of temperature at 95cm dept h highlights rapid increases in fumarole temperature that show some correlation with rainfall. No decrease in temperature was noticeable after individual rainfall events.


45 2 .5 Discussion and Conclusion Changes in magmatic activity at Masaya during Octobe r and June 2006 appear to be reflected in both temperature variations in the Comalito fumaroles and mild volcanic activity within Santiago Crater. Widening of the active vent in June 2006 and a new lava extrusion on the floor of Santiago crater in October 2006 both corresponded to rapid fluctuations in soil temperatures near Comalito cinder cone. Anomalous temperatures during volcanic activity are particularly clear on spectrograms (Fig ure 2. 3a) and suggest that fumarole temperature time series can be a sen sitive monitoring tool in some circumstances. The nature of temperature anomalies during the episodes of crater activity provide clues about the mechanisms causing these anomalies. The temperature variations were rapid, with initial increases of up to 5 C in less than 20 minutes. As bulk hydraulic conductivity is meters per day at Masaya volcano (MacNeil et al., 2007), transfer of heat by advection of hot fluids from Santiago crater, or even from a closer source, to fumaroles near Comalito is not a plausi ble source of these temperature anomalies. The flow of groundwater toward the active vent, rather than away (MacNeil et al., 2007), makes direct advection of hot fluids from Santiago crater even less tenable. Instead, we suggest that rapid increases in te mperature may be a response to pore pressure increase in the fracture zone. A relevant observation is the durations of temperature anomalies during volcanic activity, which are consistently ~1500 min. The hydraulic conductivity ( K ) of fractured basaltic la vas in the edifice is ~ 1 10 3 m s 1 (MacNeil et al., 2007). Assuming specific storage ( S s ) of this aquifer is typical of fractured rock, ~ 1 10 5 m 1 the hydraulic diffusivity ( D = K / S s ) would be 100 m 2 s 1


46 The dimensionless response time ( ) of a hy draulic pressure wave through an aquifer is equal to the ratio Dt / L 2 where t is time and L is the distance from the pressure pulse source (perhaps beneath the volcanic vent) to the Comalito fumaroles. A full response occurs at a value of = ~1. If D = 1 00 m 2 s 1 and L = 3500 m, a full response time is ~2000 min. Given that the duration times of the events were always very close to this, ~1500 min (Fig ures 2 2b d), we think this duration time may reflect the hydraulics of the pressure pulse that is trav eling along the fracture zone from the source. We emphasize that all five fumarole soil temperature time series responded in the same way and the air temperature did not, although the magnitudes of soil temperature responses varied among the fumaroles. T hese temperature variations contained dist inctive, repetitive cycles (Figures 2 2b d) that are similar to those observed in pressure time series related to hydraulic fracturing (Enever et al., 1992). In light of this similarity we propose a similar mechan ism. Volcanic activity causes an increase in pore pressure throughout the nearby groundwater system. Along the fracture zone (Fig ure 2 1), this pore pressure increase dilates fractures or perhaps opens new fractures, enhancing permeability and resulting i n increased fluid transport to the surface. Rainfall also increases pore pressure (Wang and Sassa, 2003) and appears to sometimes trigger the temperature response, perhaps by creating fractures (Husen et al., 2007) and/or increasing pore pressure in the fr acture zone past some threshold value. Initially this change in pore pressure causes spikes in fumarole soil temperatures as water vapor is transported rapidly to the surface, but afterwards the temperature in the fumarole decreases exponentially as ascend ing fluids reach a new equilibrium and lose heat to the surrounding rock. At the end of the cycle, pore pressure drops due to a change in the volcanic system, fractures


47 close, and fluid flow rapidly returns to background conditions. In addition, pore press ure decrease at the end of the cycle could facilitate mixing between ascending fluids and cooler, shallow meteoric waters, as proposed at Merapi Volcano (Zimmer and Erzinger, 2003). In this model, each cycle of rapid temperature rise, flow at a new higher rate, and rapid decrease corresponds to pressure changes in the magmatic system. The distal fumarole soil temperature data suggest that the June activity was characterized by 13 cycles of pressure changes and the October activity by 11 cycles. Such cyclic pressurization of the magmatic system during episodes of volcanic activity may be characteristic of Masaya volcano.


48 CHAPTER 3 INTEGRATED GEOPHYSICAL AND HYDROTHERMAL MODELS OF FLANK DEGASSING AND FLUID FLOW AT MASAYA VOLCANO, NICARAGUA 3 .1 Introduction Surf ace diffuse degassing is a direct result of magma decompression and therefore changes in diffuse degassing can be related to changes in volcanic activity (Chiodini et al., 1998). However, the connection between surface degassing and the magmatic source region often involves groundwater which can strongly affect the location and magnitude of surface emanations that we can detect (Todesco, 2008). Volcanic terrains in particular are characterized by faults, fractures and per meability variations that strongly affect fluid flow paths and rates (Todesco, 1997; Caine et al., 1996; Evans et al., 2001; Manzocchi et al., 2008). These features can also change over time Geophysical investigations can help constrain these subsurface v ariations. Numerical simulations allow us to explore how these structures can affect fluid flow, and aid in relating surface diffuse degassing variations directly to their volcanic source. Volcanic activity perturbs groundwater flow by adding heat and gases to groundwater systems, manifest at the surface by flank degassing, development of springs, and flank fumaroles. Observation of diffuse degassing has been applied as a volcano monitoring tool with varying degrees of success at Mt. Merapi, Indones ia (Toutain et al.,


49 2009); Mt. Fuji, Japan (Notsu et al., 2006); Mt. Etna, Italy (Badalamenti et al., 2004); Phlegrean Fields, Italy (Granieri et al., 2003); Stromboli volcano, Italy (Finizola et al., 2002); Taal volcano, Philippines (Zlotnicki et al., 200 8); Ruapehu volcano, New Zealand (Werner et al., 2006); Nisyros, Greece (Brombach et al., 2001); Mammoth Mountain, USA (Rogie et al., 2001); and San Vicente volcano, El Salvador (Salazar et al., 2002). These studies all indicate that a good geological mod el of the subsurface is a vital prerequisite to interpreting changes in diffuse degassing in light of volcanic activity. In this study, therefore, we performed detailed geophysical investigations of flank degassing at Masaya volcano, Nicaragua, to constr ain numerical models of fluid transport. We collected magnetic data in a 133 m x 125 m low temperature fumarole zone on the flank of Masaya volcano, and used these data to infer local geological structures and permeability variations. These structures, and information about the hydrologic system gained from transient electromagnetic soundings (MacNeil et al., 2007), are used to constrain a TOUGH2 numerical simulation. Fluid flux output from the TOUGH2 model is compared with CO 2 flux and self potential (SP) measurements to create and refine a structural and hydrothermal model of the flank fumarole zone. The consistent spatial distribution of diffuse degassing is then used to constrain a TOUGH2 hydrothermal model of the entire fracture zone along the flank of the volcano. 3 .2 Methods for Delineating Flow Paths The geophysical techniques that we employ in this paper are each sensitive to a different rock or fluid property, and each technique provides unique information about the system. Magnetic data are par ticularly useful for detecting structures in volcanic


50 terrains, due to the high contrast in magnetic properties between basaltic lava, scoria, and alluvium (Stamatakos et al., 1997). Magnetic profiles and maps have been used successfully to infer geologica l features such as faults at a number of sites (e.g., Jones Cecil, 1995; Connor et al., 1997; La Femina et al., 2002), even in granitic areas where the magnetic contrast is relatively small (McPhee et al., 2004). On basaltic volcanoes, magnetic anomalies a ssociated with faults are often on the order of 100 1000 nT, two to three orders of magnitude larger than the magnetic anomalies induced through electrokinetic effects associated with fluid flow (Zlotnicki and Mouel, 1990; Adler et al., 1999). SP variation s result primarily from fluid flow, and therefore provide information about fluid flow paths and rates. Three different relevant mechanisms generate SP anomalies (Bedrosian et al., 2007): thermoelectric, electrokinetic and fluid disruption effects. Althoug h each can play a part, electrokinetic effects are theoretically larger (Corwin and Hoover, 1979) and most likely to dominate in a low temperature system like the flank fumaroles at Masaya. Electrokinetic interaction between moving pore fluid and the elect ric double layer at the fluid solid interface generates an electric potential gradient that is detectable with SP electrodes (Overbeek, 1952). In extensional tectonic environments volcanic gas and vapor often rise buoyantly along faults (Goff and Janik, 20 00), resulting in positive SP anomalies (Zablocki, 1976; Nishida et al., 1996; Michel and Zlotnicki, 1998). Negative SP anomalies around fissures and craters are interpreted as meteoric water recharge (e.g. Sasai et al., 1997). After water, CO 2 is the most abundant magmatic gas (Finizola et al., 2002), and variations in CO flux can therefore be used to deduce variations in total gas flux


51 (Chiodini et al., 2005). In low temperature hydrothermal systems ( 100 C) CO emissions likely reflect exsolution from groundwater or a hydrothermal aquifer (Evans et al., 2001). The permeability and porosity of the medium through which the gases travel sig nificantly affect the location and magnitude of surface outflux (Evans et al., 2001; Lewicki et al., 2004). Thus spatial variations in surface CO 2 flux can be used to infer fractures and other variations in the subsurface (e.g., Azzaro et al., 1998; Etiope et al., 1999; Finizola et al., 2002). 3 .3 Geologic Setting Masaya is a basaltic shield volcano located 20 km south of Managua, Nicaragua (Figure 3 1), and is one of the most persistently active volcanoes in Central America (McBirney, 1956; Stoiber et al., 1986; Rymer et al., 1998). It is bounded by a 12 x 5 km caldera that has been the site of Plinian basaltic eruptions during the last ~6ka (Williams, 1983 a ; van Wyk de Vries, 1993; Walker et al., 1993; Wehrmann et al., 2006; Kutterolf et al., 2007).


52 Figure 3 1 (a) Topographic image of Masaya volcano, showing the active crater, Santiago. W hite dots represent fumaroles zones (Lewicki et al., 2003). The volcano is bound ed by Masaya caldera The inset map shows the location of Masaya volcano within Nicaragua. (b) Digital elevation map of the NE flank of Masaya volcano highlighted b y the black box in (a). The geophysical surveys were carried out at the central fumaroles zone. The fracture zone that we modeled (dashed line) links all three fumaroles zones to the summit craters. Current activity is limited to an area of small pit craters within the caldera, including Masaya, Santiago, and Nindiri craters (Figure 3 1). Santiago crater is the most active of these and has been the site of very large noneruptive gas flux in recent decades (Stoiber et al., 1986; Horrocks et al., 1999; Burton et al., 2000; Delmelle et al., 2002; Duffell et al., 2003). Crater gas flux and composition are consistent over time (Horrocks et al., 1999) and imply a magma body of approximately 10 km 3 (Walker et al., 1993), fed by gas rich magma that ascends from depth, degasses, and sinks, producing an open, relatively stable system (Stix, 2007). Post caldera eruptions have also been widespread


53 elsewhere within the caldera (Walker et al., 1993; Rymer et al., 1998; Williams Jones et al., 2003). On the northeast flank of the volcano, outside of the pit crater network, there is a fault and fracture zone extending some kilometers where water vapor and CO 2 are emitted at moderate temperatures (40 80C; Lewicki et al., 2003; Pearson et al., 2008). This radial fracture zone on the northeast flank of Masaya volcano extends past Comalito cinder cone (Figure 3. 1) which apparently formed during the 1772 eruption. The fracture zone is approximately 1 00 m wide, as suggested by ground penetrating radar (GPR) and magnetic profiles, positive SP anomalies, and elevated CO 2 flux and temperature. In some areas along its length it also forms a topographic offset at the surface. There are three fumarole zones along the fracture in the 3 4 km between the summit area and Comalito (Figure 1; Lewicki et al., 2003; Pearson et al., 2008). The third fumarole zone, closest to Comalito cinder cone, hosts some of the highest known CO 2 fluxes from low temperature fumarole s (Lewicki et al., 2003). Carbon isotopes indicate that gases emitted by these fumaroles retain a magmatic component (St Amand et al. 1998; Lewicki et al., 2003). Changes in temperature at these flank fumarol es correspond to increased volcanic activity at the summit vent (Pearson et al., 2008). Within the caldera, the groundwater budget is controlled by the balance between meteoric recharge, evapotranspiration, and degassing, which occurs primarily in Santiag o crater. This budget, supported by numerical models of groundwater transport, indicates that groundwater flow in the area of Comalito and the radial fracture zone is toward Santiago crater, which acts as a groundwater sink due to vaporization (MacNeil et al., 2007). Thus, it appears unlikely that the source of heat or magmatic gases emitted from


54 flank fumarole zones is directly beneath Santiago crater. Rather, flank degassing 3 4 km from the active vent suggests the presence of a deeper, more widespread so urce of heat and gas within the caldera. Our study is primarily concerned with modeling the thermohydrologic conditions that give rise to flank fumarole zones. 3 .4 Data Collection and Results Geophysical studies of the flank fumarole area were completed over several years of field work. During this time there were no discernible changes in the location or nature of flank degassing, and volcanic unrest was limited to lava lake activity and crater collapse in Santiago crater. A total of 4508 ground magnetic measurements were gathered at the middle fumarole zone on the flank of the volcano in August 2007 (Figures 3 1 and 3 2a). They were recorded along approximately parallel transects covering an area of 133 m x 125 m. The readings were taken w ith a Geometrics, Inc. Portable Cesium Magnetometer G858, which has a sensitivity of 0.01 nT. Measurement locations were recorded using a GPS ( Leica Geosystems Inc. G S20 Professional Data Mapper ) with an accuracy of approximately 40 cm, using a GPS base st ation to improve accuracy (baseline < 1 km). Three SP and CO 2 flux profiles were completed in May 2006. The profiles were approximately parallel and between 95 and 105 m long (Figure 3 2a), depending on topography and dimensions of the fumarole zone SP measurements were referenced to a n electrode at the end of the line outside of the fumarole area and were recorded every 1 m using non polarizing Pb PbCl 2 electrodes and a high impedance voltmeter. CO 2 flux measurements were made every 2 m using a LI COR, Inc. Li 800 portable gas fluxmeter,


55 which has an error of + 6% and a range of 0 20000 ppm. An SP map was also made from data collected in 2003 and used to interpolate five transects covering th e same area as the magnetic measurements Magnetic measurements ranged between 32000 and 41000 nT, with an average of 35900 nT (Figure 3 2b; Table 3 1), consistent with the International Geomagnetic Reference Field (35400 nT, 37 inclination, 0.5 declin ation). There is a NE SW trending magnetic anomaly of 2300 nT in the eastern portion of the map area. This anomaly corresponds remarkably well to the location and trend of the topographic offset, and is interpreted as a fault. The topographic offset and ma gnitude of the magnetic anomaly become significantly larger downslope within the fumarole zone. The magnetic anomaly amplitude ranges from 1000 nT along profile d d' to 2300 nT at profile a a' (Figure 3 2b ; Table 3 1). SP measurements ranged be tween 151 and 160 mV, with peak to peak amplitudes of 135 to 210 mV along a single profile. There is a NE SW trending, 30 m wide positive anomaly of up to 140 mV along the center of the map (Figure 3 2 c). The anomaly decreases gradually to the west, mor e abruptly to the east. The eastern boundary reaches a minimum of 150 mV, and coincides with the topographic offset. The relative size of the anomalies remained approximately 200 mV between the map measurements in 2003 and the profile measurements in 2006 However, the absolute values were consistently smaller in 2006, with an average measurement of 31 mV in 2003 and 21 mV in 2006, possibly due to seasonal fluctuations in pore water content. Several SP profiles (Figures 3.3 and 3.4) show a gradual increas e and then a more rapid decrease


56 Figure 3 2 ( a) Topographic map of the central fumarole zone, showing offset extending to the NE, and locations of GPS and magnetic measurements (gray dots) and SP and CO 2 flux measurements (dashed lines 0 2). ( b) Map of magnetic results overlain on topo graphic contours. The NE trending positive magnetic anomaly corresponds to the topographic offset. Dashed lines (a e) show locations of magnetic profiles in Figure 3. (c) SP map overlain on topographic contours. A positive SP anomaly occurs NW of the topographic offset. Dashed lines (A E) show the locations of SP profiles in Figure 3. (d) Map of magnetic data plotted with SP contour lines.


57 from NW to SE (Lines 0, 1, 2 and B). The changes are concordant with fault locations estimated from topography and magnetic data. CO 2 flux measurements varied between 0 and 2200 gm 2 d 1 (Figure 3.4). Although Line 2 sh owed some agreement with SP data, CO 2 flux variations were generally more abrupt and less consistent than changes in SP. Elevated CO 2 flux was detected only NW of the topographic offset; there was a complete absence of CO 2 to the SE. Table 3 1 Statistical results from magnetic, CO 2 a nd SP surveys Observed Mag (nT) Max Min Ave SD Anomaly A 37402.5 35142.9 35843.3 568.9 2259.6 B 37239.3 35342.7 35832.8 458.4 1756.5 C 36885.0 35080.4 35798.2 453.9 1804.6 D 36911.8 35454.5 35490.3 304.7 1034.9 E 38384.4 35283.9 36383.4 744.8 2079.4 Map 40927.4 32640.7 35952.5 803.0 6000 Calculated Mag (nT) Max Min Ave SD Anomaly A 37358.0 35147.4 35844.2 548.8 2210.6 B 37252.8 35347.5 35832.7 430.4 1728.3 C 36731.1 35173.7 35796.5 397.4 1551.7 D 36905.7 35474.5 36011.0 305.4 1238.1 E 38262.1 35525.6 36384.4 647.0 1942.0 CO (g/m 3 /d) Max Min Ave SD Anomaly Line 0 1823.27 0 242.19 418.95 1823.3 Line 1 1935.17 0 375.06 523.81 1935.2 Line 2 2204.4 1.93 309.02 449.52 2204.4 SP (mV) Max Min Ave SD Anomaly A 139.79 39.37 43.69 50.22 168.3 B 129.00 26.41 44.17 45.99 143.0 C 131.29 89.58 30.53 53.86 163.9 D 122.32 40.48 20.41 43.27 138.3 E 133.06 55.00 14.51 44.33 138.7 Line 0 88.1 99.5 10.2 42.6 163.7 Line 1 59.4 150. 8 32.3 51.8 207.4 Line 2 74.5 134.0 17.3 54.1 153.2 Map 160.0 130.0 46.5 64.3 202.8


58 Figure 3 3 Magnetic and SP profiles. Red dots correspond to observed magn etism and black lines to calculated magnetism. The solid green lines are SP profiles captured from the map in Figure 3.2c. Shaded areas represent the structures estimated from magnetic data. The darker shaded unit corresponds to lava with a magnetic suscep tibility of of 2 x 10 2 (cgs units). The lighter shaded unit corresponds to an overlying veneer of scoria with a magnetic susceptibility of 1 x 10 6 Arrows suggest estimated normal faults. Broad (30 m) SP anomalies are observed to the NW side of the esti mated faults, with much smaller amplitude anomalies to the SE.


59 Figure 3 4 Profiles of SP (solid green lines) and CO 2 flux (blue dashes). Measurements were recorded along lines 0 2 in Figure 2a, and show positive anomalies in the NW part of the study area. 3 .5 Modeling and Interpretation 3 .5.1 Magnetic data The distribution and dip of faults in the study area was estimated from ground magnetic data using Geosoft Inc. Oasis Montaj softw are. The total intensity of magnetization is the sum of the intensity of induced and remnant magnetization (Plouff 1975), and these two magnetization values can be estimated separately. As we are only interested in the total vector of magnetization, we cre ated a model using bulk magnetizations. We assumed constant remnant magnetization because all of the potential magnetic field (~0.78 Ma; Cande and Kent, 1995), and tectonic r otations are thought to be minimal during this time. A peak to peak amplitude of over 1000 nT is consistent with vertical separation of basalt that has a strong magnetization (e.g., La Femina et al. 2002). We consider the


60 magnetic source to be a lava flo w unit of high bulk magnetization beneath scoria of low bulk magnetization (Figure 3 3b). Modeling is used to estimate depth to the top of lava beneath a veneer of scoria. We used vectors of magnetization of magnitude 2 x 10 2 for the lava flow and 1 x 1 0 6 for the scoria. This produces a plausible lava flow structure with a depth < 20 m and a variable depth to the top of the lava of +/ 10 m (Figure 3 3). This structural model and the observed topographic offset indicate the presence of SE dipping norma l faults in the map area. 3 .5. 2 Hydrothermal models Numerical modeling of the hydrothermal system was done using the Petrasim interface to the TOUGH2 code. TOUGH2 is a sophisticated program to model multiphase and multic omponent fluid flow, evolved from the MULKOM code (Pruess, 2004). TOUGH2 numerically simulates coupled non isothermal heat and fluid transport through porous media. Thermodynamic conditions according to the steam table equation (International Formation Com mittee, 1967) are used to calculate phase transitions, so that gases and liquids can be included. The effects of capillary pressure and relative permeability of each phase are calculated as a function of phase volumetric fraction. Full details of the model can be found in Pruess (1991). 3 .5. 2 .1 Vadose z one Our various geophysical data and models provide constraints on the shallow hydrothermal system. Transient electromagnetic soundings throughout the caldera indicate that the groundwater table is approx imately 250 m beneath the surface in the study area (MacNeil et al., 2007). Magnetic data and models are consistent with a faulted


61 terrain, dominated in this area by steeply dipping (60) normal faults. SP and CO 2 flux data indicate that these faults serve as barriers to flow, rather than enhancing flow along the fault plane. Shallow temperature measurements recorded over a period of years at five different depths at Comalito cinder cone show average temperatures of between 46 and 74C, and maximum temperat ures of between 61 and 82C (Figure 3 5 ). The average near surface maximum temperature over all the depths is 73C. Thus the hydrothermal system should be locally dominated by upward flow of fluids through the footwalls of these faults, with near surface temperatures of ~70C. Figure 3 5 Example time series of fumarole temperature recorded at Comalito cinder cone. The average near surface temperature is 70C, used as a constraint in the numerical modeling. Temperature anomalies cor respond to episodes of minor volcanic activity in Santiago crater (see Pearson et al. 2008). We created a 30 by 30 rectangular grid 500 m wid e and 250 m de ep (Figure 3 6 c). These dimensions correspond to the width of the fracture zone (and sl ightly beyond), and to the approximately 250 m deep water table (MacNeil et al., 2007). Within the grid, impermeable faults dip at approximately 60 (Figure 3 6 c). The entire model was initially at atmospheric pressure and contained moist air (99% air, 1% water vapor). The upper


62 Figure 3 6 (a) SP measurements recorded along profile C vertical fluid flux along a profile just below the surface, using the geometry in (c). (c) Geometry used in a TOUGH2 m odel to represent vadose zone conditions. Other parameters used in the model can be seen in Table 2. (d) Output of the TOUGH2 model after injecting heat at 3 W/m 2 for 250 years. Black arrows represent gas flow, and colored contours show temperature. Heat a nd fluid flux increase toward the footwall of each fault. (e) Output of TOUGH2 model after injecting heat at 4.5 W/m 2 and moist air at 1 x 10 4 kg/s for 250 years. Fluid flow and temperature increase significantly toward the footwall of each fault.


63 boundary of the model domain was fixed at atmospheric pressure, 20C, 99% air, and 99% porosity (to represent the atmosphere). Both of the vertical boundaries were impermeable to represent the limit of t he fractured/faulted zone. The interior of the model was initially 20C. The bottom boundary, corresponding to just above the water table, was at 100C. Heat and/or fluid were injected uniformly along the base of the model, with the rate of injection deter mined by comparison of model outputs with observed surface temperatures Bulk permeability was similarly determined. Other parameters (Table 3.2) were measured by Chiodini et al. (2005), MacNeil et al. (2007), or are typical of fractured basalt. The mode l was run for 250 years, the length of time since this area last saw surface volcanic activity. Table 3 2 Parameters used in numerical models Parameter Units Value Rock density kg/m 3 2000 Porosity 0.3 Permeability of vadose zone rock m 2 1 x 10 8 Permeability of fault m 2 1 x 10 20 Permeability of saturated zone rock m 2 1 x 10 10 Wet heat conductivity W/mC 1.43 Specific heat J/kgC 850 Enthalpy of water J/kg 4.19 x 10 5 Enthalpy of steam J/kg 1.27 x 10 5 After 250 years of heat injection moist air is circulating steadily, primarily controlled by the geometry of faults (Figure 3 6d ). Convection cells develop on each side of both faults, with the largest gas fluxes occurring between the two faults (Figure 3 6b ). Upward flow of gas increases towards the footwall of each fault, and becomes negligible over the faults. This effect is particularly distinct near the surface. Close to both vertical


64 boundaries gas flow s downward. Temperatures correlate with gas flux, increasing near the footwall of each fault. The maximum near surface temperature is 70C and the maximum temperature over the entire model domain is 175C. The maximum temperature occurs where hot gas es are trapped between the two faults. When fluid is injected in addition to heat, fluid flow is enhanced across the entire model. There is a more even surface temperature distribution, with a maximum near surface temperature of 70C and a maximum temperat ure of 99C overall. Gas flux and heat flow are enhanced in the footwall of the left hand fault compared to the model without mass injection (compare Figures 3 6d and e). The hanging wall of the right hand fault is isolated from fluid injection, and fl uid circulation there is the same in both models. Permeability and heat injection rate are both unknown, but can be partially constrained by the known surface temperature. The measured surface temperature of 70C can be modeled with a permeabilit y b etween 3.5 x 10 9 and 3 x 10 8 m 2 and heat injection rates of 2.4 to 9 W/m 2 for the heat injection only case (Figure 3 7a ). Within this range there is a threshold permeability for convection. Below this permeability value the relationship between heat injection and temperature is approximately linear. For convection dominated systems (k > ~ 3.5 x 10 m ) higher permeabilities require higher heat injection rates to attain temperatures of >70 C. In Figure 3. 6 we used an intermediate permeability of 1 x 10 8 m 2 and a heat injection rate of 3 W/m 2 to show representative gas flux within the system. These values are within the fully convective regime.


65 When both heat and fluid are injected, fluid flow is enhanced. Much lower permeabilities, and a wider r ange of permeabilities, can support a surface temperature of 70 C (Figure 3 7 b). To obtain the results shown in Figures 3 6 e and 3 7 b, moist air (99% air, 1% water) was injected along the base of the model at a rate of 1 x 10 4 kg/s. Above a high p ermeability threshold of about 10 m temperature decreases with increasing permeability. For a model permeability of 1 x 10 8 m 2 a larger heat injection rate (4.5 W/m 2 ) is required for the near surface temperature to reach 70C, compared to the heat in jection rate required without fluid injection (3 W/m 2 ). Figure 3 7 Near surface temperature as a function of permeability. Numbers labeling the curves are heat injection rates in W/m 2 a) Heat injection only; b) Injection of both h eat and moist air. The observed temperature of 70C (grey line) can be matched with a range of conditions.


66 The magnetic models suggest a shallow scoria layer underlain by lava. However, a uniform permeability distribution appears to give the best fit to the data. Adding a more permeable shallow layer causes the heat and gas flux to become more focused. A dding a less permeable layer results in more even temperature and gas flow rate distributions (Figure 3.8) We infer that although fractured lava and scoria have very different magnetic properties, they have similar permeabilities. Permeability and heat in put affect fluid velocities as well as the temperature and gas distribution at the surface. Maximum gas velocities vary from 5 x 10 4 m/s in the low permeability models to 2 x 10 1 m/s in a system with a high permeability of 10 7 m 2 and a heat input of 35. 5 W/m 2 This implies travel times ranging from 20 minutes to 5 days between the water table (250 m depth) and the surface. The velocity for our average model (permeability 1 x 10 8 m 2 heat injection 3 W/m 2 ) was 5.5 x 10 3 m/s, resulting in a minimum trave l time of 13 hours. The maximum surface gas flux from our average model was 1.4 x 10 2 kg/s and the maximum surface heat output was 105 W/m 2 comparable to the 91 W/m 2 measured at Comalito (Chiodini et al., 2005) Injecting fluid as well as heat did not appear to significantly affect these values.


67 Figure 3 8 Fluid flow and heat transfer when there is a strong permeability contrast between the scoria an d lava layers. Top: Overlying scoria is two orders of magnitude more permeable than lav a Bottom: Scoria layer is two orders of magnitude less permeable. 3 .5.3.2 Saturated z one The de gassing system at Masaya volcano is unusually stable and open, and the consistent spatial distribution of diffuse degassing can be used to infer spatial variations in fluid flux at depth. The three distinct fumarole zones along the fracture zo ne on the flank of Masaya volcano (Figure 3 1) are all characterized by stable and consistently


68 elevated temperature and CO 2 flux. S oil temperatures at the Comalito cinder cone fumarole range between 55 and 80C (Pearson et al., 2008), and are simi lar to those measured at the fumarole closest to the crater (38 to 83C). Measured CO 2 fluxes are similar at all three fumarole zones, varying between 2000 and 2255 g/m 2 /d. At Comalito cinder cone, measured gas fluxes are stable over remarkably long period s of time varying by less than 10% from average over 5 months. The uniformity in temperature and CO flux across the three fumarole zones suggests that they share a common source. Gases at Comalito cinder cone, 3.5 km from the active crater, retain a ma gmatic component (Lewicki et al., 2003), but as groundwater flow is from Comalito toward the active crater (MacNeil et al., 2007) there must either be separate magmatic sources, or one extensive source at depth. Lava lake formation in the active crater in 2006 corresponded to changes in temperature at the Comalito cinder cone fumaroles (Pearson et al., 2008), suggesting that there may in fact be an a re ally extensive magmatic source at depth beneath the caldera. We created a TOUGH2 model of the saturated zone, to determine whether variations in fluid flux at depth could result in the three distinct fumarole zones observed at the surface. A square 30 by 30 grid extended over 3500 m was used to represent the entire length of the fracture zone (Figure 3 9 a ). A depth of 3500 m was assumed for simplicity, as there is no information to constrain depth to the heat source. Heat and/or fluid was injected at various rates along the bottom boundary, and the entire model was fully saturated and under hydrostatic pressure (Table 3 2; Figure 3. 9 a ). Vertical boundaries were impermeable to reflect the limit of the fractured area. This model was also run for 250 years


69 Figure 3 9 (a) Geometry used for numerical model, representing a cross section of the flank fracture at Masaya (Figure 4.1b). The system is fully saturated. (b) Results show that convection of groundwater can create three zones of elevated fluid flux and tempera ture at the surface in a homogeneous fracture above a uniform heat source.


70 TOUGH2 modeling shows that injecting liquid evenly into the base of the model at 20 kg/s (2/3 kg/s into each cell) results in convection along the length of the fracture system ( Figure 3 9 b). With a permeability of 10 10 m 2 three distinct zones of elevated temperature and fluid flux result. The maximum temperature is 70C near the surface and 91C overall, with a near surface fluid flux of 4 kg/s. Similar results are found for a heat injection rate of 600 W/m 2 Higher and lower injection rates result in hotter and cooler temperatures, respectively. The number of plumes can also vary. Permeability and model and dimensions similarly affect maximum temperature and the number of plum es created, for example a model that is twice as wide has twice as many plumes. Constant heating at depth could potentially cause convection in both the vadose and saturated zones. Simple calculations can show whether convection is likely along the fra cture, either within the saturated or unsaturated zone. The onset of convection depends on the porous media Rayleigh number (Zhao et al., 2003): Substituting values appropriate for the vadose zone at Masaya volcano (Table 4. 3; Zhao et al., 2003; Chiodini et al. 2005; MacNeil et al. 2007) yields: for gases within the vadose zone. The critical Rayleigh number is given by (Zhao et al., 2003):


71 where H 1 is the length of the fracture, H 2 is its thickness, and H 3 its height. For t his two dimensional approximation to be valid, H 2 <

72 Table 3 4 Parameters used to determine whether convection will occur along the fracture at Masaya volcano Model H 1 H 2 H 3 Ra Ra crit Convect? V adose zone 3500 100 250 58 83 No Saturated zone 3500 100 3500 1.6x10 5 1.1x10 4 Yes 3 .6 Discussion Interpretation of magnetic trave rses on the flank of Masaya volcano suggest that SE dipping normal faults redirect upward fluid movement (e.g., Figure 3 3b). This has a strong effect on mass flux, soil gas concentrations, and fumarole temperature. Numerical modeling shows that relatively impermeable faults dipping at 60 focus flow in the footwalls, and inhibit flow in the hanging walls (Figures 3 6b 3 6d 3 9 ). Positive SP and CO 2 anomalies to the NW, in the footwall of the main fault, are in excellent agreement with this model, as is a complete absence of CO 2 and a small SP anomaly to the SE (Figures 3 3, 3 4, 3 6 ). The small negative SP anomaly to the SE could result from circulation in the hanging wall, as seen in the model, or from near surface permeable sedim ents and/or fractured bedrock channeling fluid away from the fault. However, the second possibility cannot explain the large positive SP anomaly observed in the footwall (Goff and Janik, 2000; Hochstein and Browne, 2000). Normal faults dipping at 60 and extending to the water table can explain the distinct distribution of CO 2 and SP anomalies observed at the surface. More steeply dipping faults would result in sharper and more localized anomalies than we observe (Figures 3. 4 and 3. 5a; Goff and Janik, 2000 ; Hochstein and Browne, 2000). For a water table at 250 m depth, a 60 normal fault would channel fluids to the surface from an


73 approximately 150 m wide area. This is generally consistent with the dimensions of the top of a single convection cell in our sa turated zone modeling (Figure 3 8c ). Figure 3 10 Conceptual model of the hydrothermal system at Masaya volcano, as determined by geophysical observations and numerical models. On a km scale, groundwater convection is a dominant control on fluid flux. Uniform heat injection at depth creates three zones of elevated fluid flux and temperature. Within the vadose zone, relatively impermeable faults redirect shallow fluid flow. The line marked with an open triangle represents the wa ter table. Numerical models are influenced by system geometry, heat injection rate, and rock permeability. Although the fracture geometry at Masaya volcano is constrained by surface features and GPR and magnetic surveys, heat injection rate and permeabili ty are unknown. Permeability of basalt is typically in the range of 10 14 to 10 9 m 2 (Freeze and Cherry, 1979; Ingebritsen and Scholl, 1993; Saar and Manga, 1999), with young, unaltered basalt in the range of 10 11 to 10 9 m 2 (Ingebritsen, 2009). Our infe rence of 3 x


74 10 9 to 3 x 10 8 m 2 within the vadose zone is reasonable for the highly fractured basalt of our study site. Our inferred saturated zone permeability of around 10 10 m 2 is within the range of measured basalt permeabilities worldwide. It is orde rs of magnitude lower than our inferred vadose zone values, consistent with the fact that the permeability of basalt decreases with depth (Ingebritsen and Scholl, 1993). The fracture zone is however still relatively permeable, possibly explaining why it ho sts the only extensive hydrothermal feature observable outside of the crater at Masaya volcano (MacNeil et al., 2007). Such focused degassing has also been observed at other volcanoes (Chiodini et al., 2001). Although faults are often assumed to be permeab le conduits to flow (e.g., Evans et al., 2001), the faults in our vadose zone models are relatively impermeable. Previous studies have shown that permeability can vary by several orders of magnitude within a fault (Manzocchi et al., 2008), and/or can creat e a predominantly low permeability zone with channels of high permeability (Marler and Ge, 2003; Fairley et al., 2003). Permeability is a function of fracture dilatancy, cementation, compaction, clay formation, deformation stress and temperature histories, micro cracking of dense rock, chemical alteration, diagenesis, and mineral precipitation (Goff and Gardner, 1994; Goff and Janik, 2000; Fisher and Knipe, 2002; Wibberley and Shimamoto, 2003; Bernabe et al., 2003). Low permeability fault gouge may be prese nt, as at Elkhorn Fault Zone, South Park Colorado, USA (Marler and Ge, 2003). Fault geometry can also be important. Caine et al. (1996) showed that a fault zone can be divided into fault core, fractured damage zone, and undeformed protolith. The fault core may be clay rich, brecciated and/or geochemically altered, resulting in a low permeability. If the fault core is large relative to the fault as a whole, the overall permeability will be low. In addition, Caine et al. showed


75 that small damage width of a fa ult (relative to its total width) will result in localized strain which inhibits flow, acting as a low permeability zone. There are therefore a number of ways that the faults at Masaya volcano could behave as low permeability zones that redirect fluid flow To recreate the surface temperatures observed at Masaya volcano we had to estimate both permeability and heat injection rate. When heat alone is injected into the base of the vadose zone, vadose zone permeabilities of 3 x 10 9 to 3 x 10 8 m 2 and heat i njection rates of 2.4 to 9 W/m 2 are required to match observed surface temperatures of ~70 C (Figure 3. 7a ). When both moist air and heat are injected into the base of the vadose zone, vadose zone permeabilities of <1 x 10 8 m 2 and heat injection rates o f 4.5 to 6 W/m 2 are required to match observed surface temperatures (Figure 3. 7b ). Our preferred heat injection rates of 2.4 to 9 W/m 2 are less than the measured heat output of 91 W/m 2 at Comalito cinder cone (Chiodini et al., 2005), but our model exten ds beyond the surficial degassing area, and maximum values are comparable. Gas velocities within the vadose zone suggest travel times on the order of hours to days. This is compatible with observations; fumaroles are known to respond to changes at depth, but there is generally a time lag. The maximum gas flux from our average model was 1.4 x 10 2 kg/s, orders of magnitude less than the 400 kg/s measured in the active crater several kilometers away by Burton et al. (2000). The maximum heat output was 105 W/ m 2 within 20% of that measured by Chiodini et al. (2005) at Comalito Cinder cone and representative of volcanoes worldwide (Harris and Stevenson, 1997). Numerical modeling shows that even with totally uniform fluid injection, three distinct and stable fu marole zones can develop along a fracture zone (Figure 3.8 c).


76 However, extremely high heat injection rates would be required for natural convection to take place. Injecting fluid results in forced convection with a simulated surface gas flux of 4 kg/s, tw o orders of magnitude smaller than that measured by Burton et al. (2000) in the active crater. The total heat output of 250 MW along the 1 m wide fracture zone is much greater than the 0.9 MW estimated for one fumarolic zone at Comalito (Chiodini et al., 2 005). This value is typical of much larger degassing fields (e.g., Karapiti geothermal field, Hochstein and Bromley, 2001; Hakone volcano, Iriyama and Oki, 1978; Solfatara volcano, Chiodini et al., 2001). Smaller heat fluxes can be created by our models, b ut result in lower temperatures. Our simplified two dimensional saturated zone model does ignore potentially important complications such as topography, variations in hydraulic properties, and lateral heterogeneities along and within the fracture. For exam ple, magnetic and GPR measurements suggest that the width of the fracture zone varies along its length. Work by Mheust and Schmittbuhl (2001) and Neuville et al. (2006) has shown that fracture roughness can enhance or inhibit flow. We think that these com plex heterogeneities play a role in channeling groundwater and gas flow, but that convection is also a fundamental control. 3 .7 Conclusions Geophysical measurements and numerical modeling show that fluid flux at Masaya volcano is controlled by faults on a local scale at very shallow levels (10s of meters), and by groundwater convection in a flank fracture on the kilometer scale. Magnetic measurements in 2007 detected a NE SW trending anomaly of 2300 nT to the east of our study area that coincides with sur face topographic offset. Modeling of this


77 data suggests relatively impermeable faults dipping SE at 60. Numerical models including the faults show fluid flux increasing toward the footwalls and becoming negligible over the faults. Small amounts of circula tion also occur in the hanging wall of one modeled fault This is in excellent agreement with elevated SP and CO 2 flux to the NW of the fault, and small SP values and an absence of CO 2 to the SE. Numerical models show that the three fumarolic zones observed within the 3 4 km fracture on the flank of Masaya volcano can be explained by convection of groundwater within the saturated zone. The number of convection cells and their dimensions are dependent on fracture geometry, rock permeability, an d heat injection rate. The picture is further complicated by heterogeneities within the fracture zone, including faults that control fluid flux within the vadose zone and result in complicated degassing patterns at a local scale. From our geophysical surve ys and models we conclude that within a generally high permeability fracture zone on the flank of Masaya volcano there are faults and other variations that create localized low permeability zones. Faults at the central fumarolic zone redirect fluid and hea t from a ~150 m wide diffuse source that results from groundwater convection below the water table. Relatively uniform heat and fluid transport from depth is focused by permeability variations at several scales in the upper ~3 km of the caldera section. Fl uids take hours to days to travel from the water table to the surface. A good understanding of the subsurface geology is therefore vital when attempting to understand the sources of both spatial and temporal variations in diffuse degassing.


78 CHAPTER 4 NUMERICAL MODELING OF FLANK FUMAROLE TEMPERATURE VARIATIONS RELATED TO VOLCANIC ACTIVITY 4 .1 Introduction Fumarole temperature measurements are used as a volcano monitoring tool around the globe. They have three distinct advantages over other methods; 1) fum aroles result directly from degassing of magma 2) measurements are relatively inexpensive 3 ) flank fumaroles can provide easier and safer access than crater si tes However, evidence of a clear link between volcanic activity and fumarole temperatures has proved elusive (e g. Barquero, 1988 ; Zimmer and Erzinger, 2003; De Gregorio et al., 2007). Observed variations in fuma role temperature have generally been attributed to changes in gas flux (Richter et al., 2004; Chiodini et al., 2007; De Gregorio et al., 2007), or to variations in mixing between magmatic gases, meteoric waters and hydrothermal waters (Tedesco et al., 1991 ; Todesco et al., 1997; Zimmer and Erzinger, 2003). However, these conceptual models are difficult to test and may be indeterminate Numerical modeling provides a powerful tool to relate measured surface temperature variations to cha nges in the hydrothermal system at depth. However, this has been limited by the lack of eruptive volcanic activity during sampling periods (e g. Tedesco et al., 1991; Connor et al., 1993 b ; Zimmer and Erzinger, 2003), or by sampling


79 rates on the ord er of hours ( De Gregorio et al., 2007) to weeks (Granieri et al., 2003) to months (Barquero, 1988). Steady state models of gas flow at fumaroles ha ve shown that surface temperatures are strongly dependent on source depth, mass flow rate, gas velocity, con duit/fracture geometry and surrounding rock properties ( Connor et al., 1993 b ; Stevenson, 1993 ), but have been unable to describe time varying properties, for example the dependence of fumarole temperature (and mass flow) on atmospheric pressure at Colima volca no, Mexico (Connor et al., 1993 b ). More detailed models of fumaroles have been created that attempt to simulate time varying parameters. Todesco et al. (2003) created a model of the hydrothermal system at Phlegrean Fie lds Italy based on gas composition which showed that permeability is a major control on surface diffuse degassing. Changes in the source fluid would also be reflected at the surface, but only after several years. Increased injection was shown to ca use rock deformation that is detectable much more rapidly than changes in the gas composition (Todesco et al., 2004). Of particular interest, however was the ability to model changes in the hydrothermal system that corresponded to variations in gravity sig nals recorded at the Phlegrean Fields. By varying the injection rate of gases, their composition, and the duration of injection, both the variable surface gas composition and gravity trends could be replicated (Todesco and Berrino, 2005). These models show ed that the changes observed can be explained by a pulsing source of hot fluid, which releases large amounts of CO 2 rich fluid and affects hydrothermal circulation (Todesco and Berrino, 2005). However, a more recent change in gravity and gas enrichment cou ld not be explained numerically by this model reflecting the complexity of the volcanic system (Todesco et al, 2006).


80 Here, we model variations in fumarole temperature recorded at Masaya volcano, Nicaragua. We recorded temperatures every five minutes at a low temperature flank fumarole, including during episodes of minor volcanic activity at the crater (Pearson et al., 2008). This provides an exceptional opportunity to study the details of the relationship between fumarole temperature and volcanic acti vity. We used the multiphase, multicomponent fluid flow code TOUGH2 to create a numerical model of the system, using constraints provided by previous geophysical studies (MacNeil et al., 2007; Pearson et al., submitted) During almost three years of continuous temperature monitoring we detected four episodes of very distinct temperature signals which corresponded with variation s in volcanic activity observed in the crater Comparison of these signals with TOUGH2 simulations of the hydrothermal system allows us to make quantitative estimates of heat transfer and gas flux, thereby improving ou r understandin g of the response of the hydrothermal system to changes in activity. 4 .2 Geologic s etting Masaya volcano in Nic aragua (Fig ure 4 1) is a shield volcano remarkable for its stable, open and persistently degassing system. The active vents lie within a 12 x 5 km caldera, formed at least in part during Plinian eruptions during the last ~6 ka (Williams, 1983 a ; van Wyk de Vries, 1993, Wehrmann et al., 2006; Kutterolf et al., 2007). The most recent voluminous eruption occurred in 1772, when lava flows covered the northern part of the caldera floor. Comalito cinder co ne is thought to have formed at th at time (Rymer et al., 1998). Degassing is currently largely localized at Santiago crater (Figure 4 1) where a small lava lake is sometimes visible ( Stoiber et al. 1986 ; W alker et al., 1993 ;


81 Delmelle et al. 2002 ) During the last decades activity has been characterized by small ballistic eruptions from Santiago crater and the occasional presence of the lava lake (H orrocks et al. 1999 ; Duffell et al., 2003; Stix, 2007 ). Figure 4 1 a) Topographic image of Masaya volcano. Grey lines correspond to groundwater contours at depth (in m asl ) with the potentiometric low centered about Santiago crater (MacNeil et al., 2007). Arrows show the inferred directions of groundwat er flow. The s tar represents Comalito cinder cone. Inset map shows the location of Masaya volcano in Nicaragua. b) Aerial photograph showing location of studied fumaroles adjacent to Comalito cinder cone on the flank of Masaya (Instituto Nicaraguense de Es tudios Territoriale s; W. Strauch, 2006, personal c ommun). Dashed line respresents the fracture zon e. To the northeast of the active summit complex, 3.5 km from and 200 m below Santiago crater, is Comalito cinder cone (Figure 4 .1) The two areas are linked by a 3 4 km fracture that is ~100 m wide as suggested by GPR, magnetics, and a topographic offset observed at some points along its length (Pearson et al., in review ). The fracture zone hosts outflow of CO 2 and elevated temperatures (~40 80C) at a number of distinct


82 sites along its length (Lewicki et al., 2003; Pearson et al., 2008). One of these fumarole zones, 25 0 m to the northwest of Comalito cinder cone, has some of the highest known CO 2 fluxes from low temperature fumaroles (Lewicki et al., 2003). Carbon isotopes indicate that these gases retain a magmatic component (St Amand, 1999; Lewicki et al., 2003). Pear son et al. (2008) used spectral analysis of time series to show that variations in fumarole temperature at this site correspond to changes in volcanic activity at Santiago crater (Figure 4 1 b ) The c aldera boundary has comparatively low permeability, which essentially isolates the local hydrologic system ( Figure 4 .1a; MacNeil et al., 2007). Meteoric recharge in the caldera is balanced by ground water vaporization at Santiago vent and a long the fracture zone and by evapotranspiration, which includ es evaporation from Laguna Masaya in the east As a result, groundwater in much of the caldera, includin g near the Comalito fumarole zone, flows toward the Santiago vent ( Figure 4 .1a; MacN eil et al., 2007) Therefore flank degassing is unlikely to result from direct flow of ga ses at shallow depth from Santiago crater, but rather from a deeper, more widespread source. Numerical modeling shows that the distinct fumarole zones along the fracture could result from convection of groundwater at depth, with gases being redirected loc ally by faults and similar permeability variations within the caldera section as they travel through the vadose zone to the surface (Pearson et al., in review ). Flank f umarole temperature s therefore result from a combination of volcanism, gro undwater flow and subsurface geology and are likely to have complex variations that require both conceptual and numerical modeling to explain.


83 4 .3 Fumarole t emperature d ata Fumarole temperatures were recorded to five different depths at Comalito cind er cone from May 2006 to November 2009 (Fig ure 4 2). During this time there were five distinct periods of temperature anomalies, identified primarily from interpretation of the frequency spectrum ( Pearson et al., 2008 ; see Appendix C ). The first two, in June and October 2006 (Figure 4 3 a and b) corresponded extremely well with a period of vent widening and elevated degassing in the crater, and with development of a small lava lake, respectively. Almost exactly one year later, in June and October 2007 very similar temperature anomalies were recorded (Figure 4 3 c and d) Volcanic activity is less obviously correlated, but in July 2007 the degassing changed to a more explosive, ash rich phase and in October 2007 new fumaroles formed in the active crater The fifth temperature excursion, on 29 May 2008, was around the same time as a n ash explosion immediately a fter Tropical Storm Alma made landfall less than 100 km away (Fig ure 4 2 ; Global Volcanism Report, 2010; Tenorio et al., 2010) It was the first tropical storm in recorded history to make landfall on the Pacific coast of Nicaragua (Knabb and Blake, 2008) The intense rainfall can be seen in F igure 4 3 e as can a temporary cooling of the fumarole zone that becomes less distinct with depth. The extremely rapid cooling that was only observed at the shallowest thermocouple and the one 7 m away likely resul ted from pooling of rainwater as sealant decay allowed water into the tubes containing the sensors. The thermocouple 7 m away was on a hill in a different area than the other sensors, and its temperature signals changed noticeably with respect to the oth er thermocouples after it was vandalized in July 2007.


84 Figure 4 2 Flank f umarole temperature measurements. Numbers represent the depth of each thermocouple in cm. Gr e y lines highlight the episodes with anomalous temperatures.


85 The fumarole temperature signal s during the first four anomalies ha ve a repetitive and distinct signature that provides constraints on its source. At background state the fumarole temperatures have diurnal and semi diurnal fluctuations, controlled by solar and tidal variations (Fig ure 4 3 f ). During the temperature anomalies the ir diurnal signal is lost. Instead, there is a rapid increase in temperature over 20 minutes, followed be an exponential decrease, and then a rapid decrease to background level (Fig ures 4 3 a d ). Each of these cycles lasts approximately one day, with 10 16 cycles per episode The cycles were observed at each thermocouple, although the amplitude of the responses varied between thermocouples Rainfall appears to be a contributing factor as 3 4 % of the rapid increases in temperature coincided with rainfall events (Fig ure 4 .3a d). However, 405 of the 611 temperature spikes did not occur during rainfall and there were 7237 rainfall events (out of 7443) that did not have associated temperature spike s Simil arly, significant rainfall events like Tropical Storm Alma and the first rainfall of the rainy season d id not have appear to be related to significant temperature increases (Figure 4 .3e and f). The timescales and amplitu des of the fumarole temperature variations their repetitive but non diurnal nature, and their consistent return to background level provide important clues about the change in the hydrothermal system during these anomalous episodes. For example, cha nge s in temperature were far too rapid (both increasing and decreasing over about 20 minutes) to be associated with advection of fluids from near the active crater. TOUGH2 models were developed to explore the heat and ma ss transfer conditions that could gi ve rise to these cycles.


86 Figure 4 3 Details of fumarole temperature measurements. Rainfall is shown along the bottom of each graph. Depths of the rmocouples are as in F igure 5. 2. a) Duri ng vent widening in the crater; b ) During lava extrusion in the crater; c ) Prior to the degassing becoming more ash rich ; d ) During new fumarole activity in the crater ; e) During Tropic al Storm Alma; f ) During the first rainfall of the season after several months with no rain 4 .4 Numerical modeling The Petrasim interface to the TOUGH2 code was used to simulate the fumarole system at Comalito cinder cone. TOUGH2 is a sophisticated program to model multiphase and multicomponent fluid flow, based on the MULKOM code (Pruess, 2004). TOUGH2 numerically simulates coupled non isothermal heat and fluid transport through porous media. Thermodynamic conditions according to the steam table equation s (International Formation Committee, 1967) are used to calculate phase transitions,


87 allowing gases and liquids to be included. The effects of capillary pressure and relative permeability of each phase are calculated as a function of phase volumetric fraction. Full details can be found in Pruess (1991). The distinctive signals recorded during the anomalous temperature episodes (Fig ures 4 3 a d ) provide constraints on the mechanisms to include in numerical models In each cycle the initial temperature increases were on the order of minutes, and the entire temperature response generally took less than one day. As flow of groundwater is on the order of meters per day (MacNeil et al., 2007) it would take months for groundwater to flow from the crater to the fumarole, even were it not for the fact t hat groundwater flows predominantly toward the crater rather than away. Therefore it is unlikely that the anomalies a re due to direct flow of gases or heated groundwater from near the crater Similarly, flow of fluids from a n e xtensive magma source at several kilometers depths would take a minimum of hours but more likely days to reach the su rface. For fumarole temperatures to increase by up to 10C over 20 minutes the source of these anomalies must be shallow, either within th e vadose zone or a t the groundwater table at 250 m depth (MacNeil et al., 2007). Numerical models allow us to explore possible sources of th ese local phenomen a Previous st udies at Comalito cinder cone have shown that within the high permeability fracture zone, the fumarole zone is a site of constantly elevated gas flux and temperature (Lewicki et al., 2003, Pearson et al., 2008). Gas and heat flux are concentrated along nea r vertical profiles controlled by relatively impermeable faults and other small scale permeability variations (Pearson et al., in review ). Therefore we created a model perpendicular to the NE trending fracture zone, focusing on the centra l 1 m wide


88 ( assumed for simplicity) high permeability zone and the surrounding few meters (Figure 4 .4a) As this is all within the fracture zone the entire area still hosts outflow of CO 2 and elevated temperatures, but the higher permeability means that flow is focused along th e central zone and decreases rapidly with distance from it. To simulate this cross fracture model we created a 30 by 30 grid, spanning 30 m in width and 250 m in depth (Fig ure 4 4a). 250 m is the depth to the water table near Comalito cinder cone based on transient electromagnetic soundings (MacNeil et al., 2007) and forms the lower boundary of our vadose zone model The model has a permeab ility of 1 x 10 12 m 2 except for a 1 m wide zone with permeability of 1 x 10 8 m 2 that corresponds to the fracture (Pearson et al., in review) The entire model was initially at atmospheric temperature and pressure, co ntaining 99% air to represent the moist air in the vadose zone. The vertical edges of the model space had no flow conditions. The upper boundary cells represent the atmosphere, and were multiplied by a volume factor of 10 50 so that they h ad a fixed temperature and pressure but allowed flow to pass from the ground into them. The entire lower boundary was the site of injection of uniform gas at 0.5 kg/s The temperature of the injected gas is an important variable, and is determined by its enthalpy. For simplicity, we used air to simulate volcanic gases. As the system is below 100C, steam is only present as moisture content within the air. Enthalpy is a thermodynamic measur e of the heat content of a system approximated by:

PAGE 100

89 Figure 4 4 a ) TOUGH2 model geometry Flow is concentrated vertically along the NE t rending, high permeability zone at the center of the fracture zone The model is perpendicular to the fracture zone, with grey cells correspond ing to th e high permeability core Flow decreases rapidly moving away from this zone Modeled s urface fumarole te mperature is from the fracture core cell below the boundary ( output cell ). b ) Conceptual model of changes in the system that produce the temperature curves to the right.

PAGE 101

90 where c pa is the specific heat capacity of air at constant pressure (1.006 x 10 3 J/kg between 100C and 100C), T is the air temperature, x is the humidity ratio (0.01 for our model), c pw is the specific heat capacity of water vapor at constant pressure (1.84 x 10 3 J/kg) and h we is the evaporation heat (2.501 x 10 6 J/kg). For our system at 70C, the enthalpy is theoretically 9.7 x 10 4 J/kg. However, enthalpy is strongly dependent on pressure and this equation is over simplified for our system. In TOUGH2 enthalpy is calculated from the steam table equations (International Formation Committee 1967) We used model based temperatures, varying the enthalpy until the near surface model temperature was similar to fumarole temperatures observed at Comalito. The average measured fumarole temperature of ~70C (Fig ure 4 2) could be replicated in the n ear surface cells of the model with an inject ed gas enthalpy of 1.496 x 10 5 J/ kg Other parameters we re taken from Chiodini et al. (2005), MacNeil et al. (2007) or we re typical of fractured basalt (Table 4 1). The model was run until it reached steady state, and these pressure and temperature conditions were then used as initial conditions in models where gas injection rate, enthalpy, rock permeability and source depth were varied to try to replicate the cyclic temperature time signal observed during vo lcanic activity (Figures 4 .3a d)

PAGE 102

91 Table 4 1 Parameters used in TOUGH2 model. Variable Value Units Density 2000 kg/m 3 Porosity 0.3 Permeability (fracture) 1 x 10 8 m 2 Permeability (surroundings) 1 x 10 12 m 2 Thermal conductivity 1.43 W/mC Specific heat capacity 850 J/kgC Air mass fraction 0.99 Initial temperature 20 C Initial pressure 1.013 x 10 5 Pa Background enthalpy 1.496 x 10 5 J/kg Background injection rate 0.5 kg/s 4 .5 Results The goal of our modeling is to explain the very consistent but complex fumarole temperature an omalies observed at Comalito cinder cone. As the signals are rap id and repeatable changes in the system must be local and elastic. To increase surface fumarole temperature by up to 10 C over 20 minutes, a significant amount of energy must be entering the near surface system. We modeled variations in the gas injection rate, enthalpy, source depth and rock permeability to determine if they could result in this magnitude of surface temperature increase. The models must also be able to explain t he gradual decrease in surface temper ature and the rapid return to background t emperature observed after the initial increase Varying the rate and enthalpy of gas injected into the base of the model (250 m depth) can result in surface temperature variations similar to those observed at Comalito

PAGE 103

92 fumarole (Fig ure 4 5). A background injection rate of 0.5 kg/s results in a model that attains steady state in 9.4 years. The rapid incre ase in temperature from approximately 72C to 79C in 20 minutes can be sim ulated by a sudden injection of hot gas at 50 kg/s, a hundred fold increase. The exponential decrease in temperature corresponds to an almost negligible gas flow rate of 0.05 kg/s or less. The rapid drop to background level is recreated by a slightly higher injection rate of 5 kg/s, but with background temperature gases. Figure 4 5 Changes in surface fumarole temperat ure are modeled with abrupt changes in mass flux and fluid enthalpy. The b ox shows where these abrupt changes occur.

PAGE 104

93 Each thermocouple at Comalito fumarole records different temperatures (Figures 4.2 and 4.3), which can be modeled with different enthalpies of the inje cted gas. The two deepest thermocouples, at 150 cm depth, measure average temperatures of ~62C, while the shallower thermocouples recorded temperatures between ~67 and 75C. This corresponds to a modeled enthalpy of injected gas between approximately 1.4 and 1.6 x 10 5 J/kg. For example, the 72 to 79C temperatures measured at the thermocouple at 66 cm depth correspond to enthalpies of 1.52 to 1.6 x 10 5 J/kg, respectively (Figure 4.5). The 70 to 76C temperatures observed at the thermocouple at 95 cm depth correspond to enthalpies of between 1.496 and 1.55 x 10 5 J/kg (Figure 4.6). The different amplitude fumarole temperature cycles observed at each thermocouple can therefore all be recreated by using different enthalpy ranges. The shape of the temperature cy cles is sensitive to both the enthalpy and the rate of gas being injected into the base of the model (Fig ure 4 6). With a low injection rate the system takes longer to heat up and with an injection rate of 1 kg/s it does not even reach the maximum t emperature of 76C within the 2 days modeled (Figure 4 .6b left panel ) With a low enthalpy the system reaches a lower equilibrium temperature for example 75C with an enthalpy of 1 .55 x 10 5 J/kg compared to 85 C with an enthalpy of 1.65 x 10 5 J/kg (Figure 4 .6c) For the final stage of the curve, a high injection rate results in a rapid change, with a n increase of up to 0.5 C and then a larger decrease of up to 2C corresponding to an injection rate of 50 kg/s (Fig ure 5. 6b, far right). With an injection rate of 1 kg/s or less the temperature does not increase first, and the return to background state is much more gradual taking over a day The min or increase in temperature prior to the rapid decrease is sometimes observed in the fumarole temperate measurements also

PAGE 105

94 Figure 4 6 Simulation of near surface temperature modeled with abrupt changes in gas enthalpy and injection rate a ) Curve that we are trying to re create D ashed gr e y line and b ox above show w he n injection rate and enthalpy changed. b ) Relationship between gas temperature and injection ra te. Numbers on the curves show injection rate in kg/s. Enthalpy of system at different times is shown in box es in (a). c ) Relationship between gas temperature and enthalpy. Numbers on the curve s correspond to enthalpy in J/kg. Gas injection rates are shown in the box es in (a).

PAGE 106

95 (Figure 4.3a and d). The timescale of the te mperature response is therefore controlled by the injection rate, while the amplitude of the response is controlled by the enthalpy of the injected gas. An optimum combination of the gas flux and enthalpy can be found to model each cycle at each thermocoupl e. Gas fluxes and temperatures are also affected by changes in permeability (Figure 4.7), although changes in permeability alone cannot cause the amplitude of temperature variations observed at Comalito cinder cone. If the fracture permeability decreases r elative to the surrounding rock, the near surface temperatures vary more rapidly and with higher amplitude. For a low permeability of 1 x 10 9 m 2 the initial temperature response is large, up to 78C, but the more gradual decrease afterwards is at a lower temperature, around 71C for a permeability of 1 x 10 9 m 2 compared to 75C with a permeability of 1 x 10 4 m 2 (Figure 4.7a). The effect is similar but with a smaller amplitude if the permeability of the surrounding rock is lowered relative to a constant fracture permeability (Figure 4.7b). For example, a fracture permeability of 1 x 10 11 m 2 corresponds to a near surface temperature of 72C during the gradual decrease in temperature; while a fracture permeability of 1 x 10 9 m 2 results in a near surface t emperature of 73C.

PAGE 107

96 Figure 4 7 Correlation between fumarole temperature and permeability. (a ) Numbers on the curves represent fractu re permeability as a negative power in m 2 (e.g. 4 corresponds to 1 x 10 4 m 2 ). Wall rock permeability is kept constant at 10 12 m 2 (b ) Numbers represent negative power of wall rock permeability in m 2 Fracture permeability is kept constant at 1 x 10 8 m 2 Other parameters are as in Figure 4.5. In both cases higher permeability causes a lower initial tempera ture response but a higher near surface temperature during periods of low gas flux This is probably due to associated changes in mixing between the fracture fluids and those in the surrounding rock. However, higher permeability of the fractured rock cor responds to a greater contrast with the surrounding rock, while higher permeability of the surrounding rock translates to a lower permeability contrast with the fracture. This suggests that the bulk permeability is more important than the relative permeabi lities between the two rock types. W hen the permeability of the surrounding rock is close to or greater than the fracture permeability ( 1 x 10 8 to 1 x 10 7 m 2 ) however, the rapid temperature response is lost. In both cases, if the permeability is reduced below the values modeled in Figure 5. 7 the injection rate or the number of injection cells must be reduced for the model to remain stable.

PAGE 108

97 The temperature response is found to be relatively insensitive to the depth to the source (Figure 4 .8). With a dee per source the gases have more time to react with surrounding fluids and the response is slightly enhanced, but the effect is generally less than 1C. Unless the water table rose by 200 m, it is unlikely that any change would be detectable at the surface. As the water table is unlikely to rise by so much, it appears that fluctuations in the water table are not a n important factor in either the timescales or the amplitudes of the fumarole temperature cycles Figure 4 8 Relationship between surface temperature and depth to the gas sourc e. Numbers represent depth in m. Even with an increase in the water table elevation of 150 m, the change is temperature is less than 1C. 4 .6 Discussion Numerical models of the system at Comalito cinder cone suggest that the distinctive temperature si gnals observed during volcanic a ctivity result from short term

PAGE 109

98 increases in gas flow rate and enthalpy. There is a trade off between the two; increased enthalpy requires a lower gas flow rate to reach the same temperature, but that slows the re sponse time. The most realistic model corresponds to an influx of hot gas that drops to a very low flux rate after about 20 minutes. We suggest that d uring this time the more stable fluid in the less permeable surrounding rock is kept out of the fracture and pressurized by the flow of hot gas (Fig ure 4 4b). When the flux of hot gas drops to essentially zero, the gases at background enthalpy then flood back in, causing an increase in flow rate and a second smaller temperature spike at the surface. The system then quickly reequilibriates to background temperature and gas flow rates with the rock both within the fracture and outside. The large number of repetitive cycles suggests that there is a non destructive so urce at depth causing the injection of hot fluid. Faulting or fracturing is therefore not a feasible mechanism A change in the heat source of the fluid is insufficient either to cause the magnitude and timescale of fumarole temperature variations directly or to increase buoyancy driven fluid flow. We suggest that either a sudden release of gases from a storage area within the volcanic edifice or a change in the pressure distribution is a more likely source for the changes in fumarole temperature that are related to volcanic activity. For example, d egassing of the magma source at depth could result in an abrupt change in the pressure distribution that cause s a pressure wave to travel through the water table. This in turn could enhance groundwater circulation and/or boiling increasing fluid flow and fumarole temperatures Slight differences in enthalpy and gas flow rates resulting from minor variations in either process could explain the variations between the different cycles within an episode. Constant magma degassing has been inferred from tremor at

PAGE 110

99 Masaya (Mtax ian et al., 1997), but no changes in seismic energy were identifiable related to the observed volcanic activity. The correlation of surface temperature with rainfall (Figures 4 .3a d) suggests that this may be an important factor in creating the temperature cycles As rainfall events fall within the same 5 minute sampling window as increasing surface temperature, the e ffect must be extremely shallow. For rainfall to percolate through the system, either affecting the water table o r fluid flow in the vadose zone, timescales on the order of hours to days are required. R ainfall is known to increase pore pressure (Wang and Sassa, 2003), and we instead hypothesize that rainfall affects the pressure within the system perhaps increasing it past some threshold value that allows increased tra nsport of fluids to the surface T he number of rainf all events without an associated increase in temperature and the lack of correlation between the amount of rainfall and magnitude of near surface temperature changes suggest that this is not the dominant mechanism however Chemical reactions may contri bute toward changes in soil temperature, particularly as there is a positive correlation of temperature with rainfall. The heat of wetting can be significant although the soil temperature time series associated with wetting are differ ent from those observed at Masaya ( Figure 4.9; Prunty and Bell, 2005) In addition, the lack of anomalies at the beginning of each rainy season after several months without rain (Fig ure 4 3 f ) indicates that this is not a major factor in causing observed temperature changes

PAGE 111

100 Figure 4 9 Heat of wetting curve according to Prunty and Bell (2005). The temperature does not drop as observed in the fumarole temperatures. Reactions of water, particularly with sulfur dioxide (known as scrubbing), are exothermic (Stiles and Felsing, 1926) and could give off heat and affect gas fluxes at Comalito cinder cone. Chiodini et al. (2005) did n ot detect any sulfur species from the fumaroles at Comalito although the magma at Masaya vol cano is known to be a significant sulfur source (Stoiber et al., 1986; Williams Jones et al., 2003) T he decrease in sulfur dioxide flux es in the crater prior to a phreatic explosion (Duffell et al., 2003) may suggest that scrubbing is occuring at Masaya volcano. Although this could release significant amounts of heat, we discard it as a dominant mechanis m b ecause some, but not all, cycles correspond to rainfall events. It is also unlikely that the SO 2 would have degassed from the magma and traveled through the saturated and vadose zones, only to react with water near the surface. Nevertheless we do not rule out the possibility that the

PAGE 112

101 changes in soil temperature may in part result from increased flux of SO 2 and exothermic reactions during scrubbing. Another variable that may affect gas flow rates, and therefore temperature, is permeability. After earthquakes increased fumarole temperature and gas flux ha ve been observed at a number of sites (Tedesco et al., 1991; Richter et al., 2004). There were no significant changes in seismicity corresponding to our temperature anomalies, but permeability can also change unrelated to seismicity. Wohletz and Heiken (19 92) found that hot gases could cause hydraulic fracturing increasing the effective permeability and therefore increasing gas fl ux Todesco et al. (2004) showed that hydrothermal flow cause s rock deformation and an associated change in permeability. Enever et al. (1992) showed that pressure curves resulting from hydraulic fracturing could have a very similar shape to our measured fumarole temperature anomalies In their study the system was closed, and fluid was injected. Fracturing of the rock resulted in an initial rapid increase in pressure. The following exponential decrease in pressure resulted from the fluid entering the crack and the crack tip propagating. When the fluid injection was turned off the pressure dropped rapidly. Although the system is not closed at Comalito cinder cone, it is possible that cracking of the rock surrounding the fracture could be enhancing the changes in temperature As the cycles are repetitive, with consistent thr eshold temperatures the process must be elastic and so it is unlikely that entirely new fractures are f orming in each case D ilation of existing fractures is however a plausible mechanism, although f luid injection and/or a change in pressure i s still required for this to occur

PAGE 113

102 A similar pattern in soil temperature variations has been observed in high temperature fumaroles at Colima volcano, Mexico at tributed to mixing of fumarole gas es and air (Connor et al., 1993 a ) This is analogous to our model, where fumarole gases are mixing with vadose zone fluids The background gas injection rate of 0.5 kg/s translates to a flux of 4 3 t/d at the Comalito fumarole This is in good agreement with the combined CO 2 and steam flux of 49 t/d measured by Chiodini et al. ( 2005) It is rather less than the 4 kg/s of SO 2 and 400 kg/s of water vapor emitted from the active crater at Masaya ( Galle et al., 199 3; Burton et al., 2001) During volcanic eruptions, SO 2 flux rates have been observed to increase by a factor of 10 at Mt Pinatubo (Daag et al., 1996) and Mt St Helens ( Casadevall et al., 1983) During the temperature anomalies at Masaya our models suggest that the total gas input increase d 100 fold, to a flux rate comparable with large degassing fields like Campi Flegrei (Chiodini et al., 2005). This is possible, but is probably overestimate d due to contributing factors like subsurface heterogeneity, groun dwater table depth and exothermic reactions. Geochemical flux monitoring is generally c arried out at a sampling rate of a maximum of one hour, and a s our anomalies only lasted 20 minutes it is possible that gas flux monitoring would not have sampled during this time. T he flux rate may a lso be overestimated due to heat contributions from chemical reactions and heat of wetting In our models, t he total gas input during a one day anomaly is 195 tons At background levels the da il y gas input is 43 tons sugge st ing an increase in mass of between 1520 and 2430 tons during a 10 to 16 day anomaly e pisode. SO 2 flux from Santiago crater has been observed to triple during increased magmatic degassing (Sto iber et al., 1986) A 4 .5 fold increase in total gas flux during temperature anomalies related to

PAGE 114

103 increased magma degassing is therefore entirely plausible. The change in source flux may be even smaller, as t he focusing effect of the fracture zone is only modeled to shallow dept hs, and the fracture may be concentrating fluids from greater depth Our study shows that changes in the volcanic system at Masaya volcano can be detected with fumarol e temperature monitoring, and should also be detectable with continuous gas flux monitoring. However, the changes are extremely abrupt and increases in temperature and gas flux only last 20 minutes, which is less than the sampling rate of most gas flux monitoring systems and some temperature monitoring systems. As the gradual de crease in temperature and the cor responding drop in gas flux last ed longer, approximately one day, it is possible that an increase in volcanic activity at Masaya volcano would actually be detected as a decrease in gas flux at the fumarole zone at Comalito cinder cone. 4 .7 Conclusion During volcanic activity at Masaya volcano, the Comalito fumarole shows a very distinctive cyclical temperature response that our modeling suggests is due to a rapid but t ransient influx of hot gas. At background s tate, fluid is injected along the base of the model at 0.5 kg/s Hot g as is then injected at 50 kg/s to create a rapid rise in fumarole temperature over 20 minutes. This prevents the gas in the less permeable surrounding rock from entering the fracture. The flux of hot gas drops to a rate of 0.05 kg/s or less, allowing the system to return to ward equilibrium. When it is almost at equilibrium the previously trapp ed gas in the surrounding rock re enters the fracture at a slightly higher

PAGE 115

104 injection rate of 5 kg/ s The rate then gradually drops as the system returns to background state. Our models suggest that a 100 fold increase in gas flux occurs during these temperature cycles although this only lasts for approximately 20 minut es We suggest that the increased gas flux results from a pressure pulse due to pulsing of the degassing magma which increases fluid circulation and possibly increase s permeability b y dilati ng fractures This is repeatable, resulting in a number of fumarole temperature cycles. The effect may also be enhanced by hydrofracturing or increased pore pressure due to rainfall. Wetting of the soil and chemical reactions also release heat which may co ntribute to the increases in temperature observed. However, increased gas flux appears to be the simplest physical explanation for the rapid fluctuations in fumarole temperature that we observe at Comalito cinder cone. We conclude that fumarole gas and temperature monitoring can be a useful way to detect changes in a volcano hydrologic system provid ed that a sufficient ly small ( << 1 hour) sampling rate is used Fum arole temperature monitoring can also be used to infer changes in gas flux at depth. A good conceptual and/or numerical model is however vital to be able to relate an y detected changes in fumarole activity to changes in volcanic activity.

PAGE 116

105 LIST OF REFERENCES Adler, P. M., and Le Mouel, J L. (1999). Electrokinetic and magnetic fields generated by flow through a fracture zone: a sensitivity study for La Fournaise volcano. Geophy sical Res earch Let ters 26 : 795 798. Giammanco, S., Parello, F., a nd Valenza, M. (2004). Magmatic gas leakage at Mount Etna (Sicily, Italy): relationships with the volcano tectonic structures, the hydrological pattern and the eruptive activity. In Bonaccorso, A., Calvari, S., Coltelli, M., Del Negro, C., Falsaperla, S.,editors, Mt. Etna: Volcano Laboratory pages 129 146. Geophysical Monograph 143 Azzaro, R., Branca, S., Giammanco, S., Gurrieri, S., Rasa, R., and Valenza M. (1998) New evidence for the form and extent of the Pernicana fault system (Mt. Etna) from structural and soil gas surveying J ournal of Volcanol ogy and Geotherm al Res earch 84 : 143 152. Badalamenti, B., Bruno, N., Caltabiano, T., Gangi, F. D., Giammanco, S., and Salerno G. (2 004) Continuous soil CO 2 and discrete plume SO 2 measurements at Mt. Etna (Italy) during 1997 2000; a contribution to volcano monitoring Bul letin of Volcanol ogy 66 : 80 89. Barquero, J. ( 1988 ). Changes in fumarole temperatures at Volcan Poas In Newhall, C .G., and Dzurisin, D., ed itors ., Historical unrest at large calderas of the world Volume 2 pages 87 3 87 8 U.S. Geological Survey Bulletin 1855 Baubron, J. C., Allard, P., Sabroux, J. C. Tedesco, D., and Toutain, J. P. (1991). Soil gas emanations as pre cursory indicators of volcanic eruptions. Jounral of the Geological Society, London 148:571 576. Bedrosian, P. A., Unsworth, M. J., and Johnston M. J. S. (2007) Hydrothermal circulation at Mount St. Helens determined by self potential measurements Jour nal of Volcanol ogy and Geotherm al Res earch 160 : 137 146. Bernabe, Y. Mok, U. and Evans B. (2003) Permeability porosity relationships in rocks subjected to various evolution processes Pure and Appl ied Geophys ics, 160 : 927 960.

PAGE 117

106 Brombach, T., Hunziker, J C., Chiodini, G., Cardellini, C., and Marini L. (2001) Soil diffuse degassing and thermal energy fluxes from the southern Lakki Plain, Nisyros (Greece) Geophysical Researc h Lett ers 28 :6 9 72 Burton, M. R., Oppenheimer, C., Horrocks, L. A., and Franci s, P. W. (2000). Remote sensing of CO 2 and H 2 O emission rates from Masaya volcano, N icaragua Geology, 28 : 915 918 Caine, J. S., Evans, J. P. and Forster C. B. (1996) Fault zone architecture and permeability structure Geology 24 : 1025 1028. Cande, S. C ., and Kent, D. V. (1995). Revised calibration of the geomagnetic polarity timescale for the Late Cretaceous and Cenozoic Journal of Geophys ical Res earch, 100 : 6093 6095. Casadevall, T., Rose, W., Gerlach, T., Greenland, L. P., Ewert, J., Wunderman, R. and Symonds, R. (1983). Gas Emissions and the Eruptions of Mount St. Helens Through 1982. Science 221:1383:1385. Chiodini, G., Cioni, R., Guidi, M., Raco, B., and Marini L. (1998) Soil CO 2 flux measurements in volcanic and geothermal areas Appl ied Geochem istry, 13 : 543 552. Chiodini, G., Frondini, F., Cardellini, C., Granieri, D., Marini, L., and Ventura G. (2001) CO 2 degassing and energy release at Solfatara Volcano, Campi Flegrei, Italy J ournal of Geophys ical Res earch 106 : 16213 16221. Chiodini, G., Gr anieri, D., Avino, R., Caliro, S., Costa, A., and Werner, C. (2005). Carbon dioxide diffuse degassing and estimation of heat release from volcanic and hydrothermal systems Journal of Geophysical Research, 110 : B08204 Connor, C. B., Lane, S. B., and Clemen t, B. M. (1993a). Structure and thermal characteristics of the summit dome, March 1990 March 1991; Volcan Colima, M exico. Geofisica Internacional, 32:643 657. Connor, C. B., Clement, B. M., Song, X., Lane, S. B., and West Thomas, J. ( 1993 b ). Continuous m onitoring of high temperature fumaroles on an active lava dome, Volcan Colima, Mexico; evidence of mass flow variation in response to atmospheric forcing Journal of Geophysical Research 98 : 19713 19722 Connor, C. B., Lane Magsino, S., Stamatakos, J. A., Martin, R. H., LaFemina, P. C., Hill, B. E., and Lieber S. (1997) Magnetic surveys help reassess volcanic hazards at Yucca Mountain, Nevada, Eos (Transactions, American Geophysical Union), 78 : 73 78.

PAGE 118

107 Corwin, R. F., and Hoover D. B. (1979) The self poten tial method in geothermal exploration Geophysics 44 : 226 245 Daag, A. S., Tubianosa, B. S., Newhall, C. G., Tugol, N. M., Javier, D., Dolan, M. T., Delos Reyes, P. J., Arboleda, R. A., Martinez, M. L., and Regalado, T., M. (1996). Monitoring sulfur dio xide emission at Mount Pinatubo In Newhall, C. G., and Punongbayan, R. S., editors, Fire and Mud pages 409 414 University of Washington Press and PHIVOLCS. De Gregorio, S., Madonia, P., Gurrieri, S., Giudice, G., and Inguaggiato, S. ( 2007 ). Contemporary total dissolved gas pressure and soil temperature anomalies recorded at Stromboli volcano (Italy) Geophysical Research Letters 34 : L08301 Delmelle, P., Stix, J., Baxter, P., Garcia Alvarez, J., and Barquero, J. (2002). Atmospheric dispersion, environm ental effects and potential health hazard associated with the low altitude gas plume of Masaya volcano, Nicaragua. Bulletin of Volcanology 64:423 434. Dobrin, M. B., and Savit, C. H. (1988). Introduction to geophysical prospecting McGraw Hill, New York 4th edition Duffell, H. J., Oppenheimer, C., Pyle, D. M., Galle, B., McGonigle, A.J.S., and Burton, M.R. ( 2003 ). Changes in gas composition prior to a minor explosive eruption at Masaya Volcano, Nicaragua Journal of Volcanology and Geothermal Research 1 26 : 327 339. Edge, A. B., and Laby, T. H. (1931). The Principles and Practices of Geophysical Prospecting Cambridge Univ. Press, New York. Elder, J. W. (1967). Transient conv ection in porous medium. Journal of Fluid Mechanics 27:609 623. Enever, J. R., Co rnet, F., and Roegiers, J. C. ( 1992 ). ISRM commission on interpretation of hydraulic fracture records International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts 29 : 69 72. Etiope, G., Beneduce, P., Calcara, M., Favali, P., Frugoni F., Schiattarella, M., and Smriglio G. (1999) Structural pattern and CO 2 CH 4 degassing of Ustica Island, southern Tyrrhenian Basin J ournal of Volcanol ogy and Geotherm al Res earch, 88 : 291 304. Evans, W. C., Sorey, M. L., Kennedy, B. M., Stonestrom, D. A., Rogie, J. D., and Shuster D. L. (2001) High CO 2 emissions through porous media: transport mechanisms and implications for flux measurements and fractionation Chem ical Geol ogy 177 : 15 29.

PAGE 119

108 Fairley, J., Heffner, J., and Hinds J. (2003). Geostatistical evaluation of permeability in an active fault zone. Geophysical Res earch Lett ers 30 : 1962 1966. Finizola, A., Sortino, F. Lenat, J., and Valenza M. (2002) Fluid circulation at Stromboli Volcano (Aeolian Islands, Italy) from self potential and CO 2 survey s J ournal of Volcanol ogy and Geotherm al Res earch 116 : 1 18. Fisher, Q. J., and Knipe R. J. (2002) The permeability of faults within siliciclastic petroleu m reservoirs of the North Sea and Norwegian Continental Shelf Marine and Petroleum Geology 18 : 106 3 1081. Freeze R. A., and Cherry, J. A. (1979) Groundwater Prentice Hall, Englewood Cliffs. Friedel, S., Byrdina, S., Jacobs, F., and Zimmer, M. ( 2004 ). Self potential and ground temperature at M erapi volcano prior to its crisis in the rainy season of 20 00 2001 Journal of Volcanology and Geothermal Research 134 : 149 168. Galle, B., Oppenheimer, C., Geyer, A., McGonigle, A. J. S., Edmonds, M., and Horrocks, L., (2003). A miniaturised ultraviolet spectrometer for remote sensing of SO 2 fluxes: a new tool fo r volcano surveillance. Journal of Volcanology and Geothermal Research, 119:241 254. Global Volcanism Program (2010). Masaya volcano monthly reports. http://www.volcano.si.edu/world/volcano.cfm?vnum=1404 10=&volpage=var accessed June 28 2010. Goff F., and Gardner J. N. (1994) Evolution of a Mineralized Geothermal System, Valles Caldera, New Mexico Econ omic Geology, 89 : 1803 1832. Goff F., and Janik C. J. (2000) Geothermal systems I n Encyclopedia of Volcanoes Sigurdsson, H. editor, pages 817 834 Els evier, Academic Press New York. Granieri, D., Chiodini, G., Marzocchi, W., and Avino, R. (2003). Continuous monitoring of CO 2 soil diffuse degassing at Phlegraean F ields ( I taly): Influence of environmental and volcanic parameters Earth and Planetary Scien ce Letters 212 : 167 179. Harris, A. J. L., and Stevenson D. S. (1997) Thermal observations of degassing open conduits and fumaroles at Stromboli and Vulcano using remotely sensed data J ournal of Volcanol ogy and Geothem al Res earch 76 : 175 198. Hashimoto, T. and Tanaka, Y. (1995). A large self potential anomaly on Unzen volcano, Shimabara Peninsula, Kyushu Island, J apan Geophysical Research Letters 22 : 191 194.

PAGE 120

109 Henry, H. R. (1964). Effects of dispersion on salt encroachment in coastal aquifers. U.S. Geological Survey Water Supply Paper, 1613 C :C71 C84. Hochstein, M. P., and Browne, P. R. (2000). Surface manifestations of geothermal systems with volcanic heat sources. In Encyclopedia of Volcanoes Sigurdsson, H. pages 835 855 Elsevier, Academic Press New York. Hochstein, M. P., and Bromley, C. J. (2001). Steam cloud characteristics and heat output of fumaroles. Geothermics, 30: 547 559. Horrocks, L., Burton, M. Francis, P., and Oppenheimer C. (1999) Stable gas plume composition measured by OP FTIR spectroscopy at Masaya Volcano, Nicaragua 19 98 1999 Geophys ical Res earch Let ters 26 : 3497 3500. Hurwitz, S., Kipp, K. L., Ingebritsen, S. E., and Reid, M. E. ( 2003 ). Groundwater flow, heat transport, and water table posi tion within volcanic edifices; I mplications for volcanic processes in the C asca de R ange Journal of Geophysical Research 108 : 2557 Husen, S., Bachmann, C., and Giardini, D. (2007). Locally triggered seismicity in the central Swiss Alps following the large rainfall event of August 2005. Geophysical Journal International 171:1126 113 4. Ingebritsen, S. E., and S c holl M. A. (1993) The hydrogeology of Kilauea volcano Geothermics, 22 : 255 270. Ingebritsen, S. E., Galloway, D. L., Colvard, E. M., Sorey, M. L., and Mariner, R. H. (2001). Time variation of hydrothermal discharge at selecte d sites in the western United States: Implications for monitoring. Journal of Volcanology and Geothermal Research 111:1 23. Ingebritsen, S. E. (2009). The hydrogeology of basaltic terrane. Geological Society of America Abstracts with Programs, 41: 36. Inte rnational Formation Committee (1967) A Formulation of the Thermodynamic Properties of Ordinary Water Substance ICF Secretariat, Dsseldorf. Jones Cecil, M. (1995). Structural controls of Holocene reactivation of the Meers fault, southwestern Oklahoma, from magnetic studies. Geological Society of America Buletin., 107 : 98 112.

PAGE 121

110 Knabb and Blake, ( 2010 ) Tropical Storm Alma Disc ussion Six: National Hurricane Center Retrieved 2010 03 16, http://www.nhc.noaa.gov/archive/2008/ep01/ep012008.discus.006.shtml? Kutterolf, S., Freundt, A., Prez, W., Wehrman, H., and Schminke, H U. (2007). Late Pleistocene to Holocene temporal successio n and magnitudes of highly explosive volcanic eruptions in west central Nicaragua. J ournal of Volcanol ogy and Geotherm al Res earch 163 : 55 82. La Femina, P., Connor, C. B., Stamatakos, J. A., and Farrell D. A. (2002) Imaging an Active Normal Fault in Allu vium by High Resolution Magnetic and Electromagnetic Surveys Environ mental and Eng ineering Geosci ences 8 : 193 207. Lehto, H. L. (2007). Self Potential Anomalies and CO2 Flux on Active Volcanoes: Insights from Time and Spatial Series at Masaya, Telica, and Cerro Negro, Nicaragua. Ms Thesis, University of South Florida, Tampa, Fl, USA. Lewicki, J. L., Connor, C., St Amand, K., Stix, J., and Spinner, W. ( 2003 ). Self potential, soil CO 2 flux, and temperature on Masaya volcano, Nicaragua Geophysical Research L etters 30 : 1817. Lewicki, J. L., Hilley, G. E., and Connor, C. (2004). The scaling relationship between self potential and fluid flow on Masaya Volcano, Nicaragua. Proceedings of the Eleventh international symposium on Water rock interaction, 11 : 153 156. MacNeil, R. E., Sanford, W. E., Connor, C. B., Sandberg, S. K., and Diez, M. ( 2007 ). Investigation of the groundwater system at Masaya Caldera, Nicaragua, using transient electromagnetics and numerical simulation Journal of Volcanology and Geothermal Rese arch 166 : 217 232. Manzocchi, T., Heath, A. E., Palananthakumar, B., Childs, C., and Walsh, J. J. (2008). Faults in conventional flow simulation models; a consideration of representational assumptions and geological uncertainties Petroleum Geosci ences 14 : 91 110. Marler, J., and Ge, S. (2003). The permeability of the Elkhorn fault zone, South Park, Colorado. Ground Water, 41 : 321 332. Mastin L. G. (1991) The roles of magma and groundwater in the phreatic eruptions at Inyo Craters, Long Valley Caldera, Cal ifornia Bulletin of Volcanology 53: 579 596. Mather, T. A., Pyle, D. M., Tsanev, V. I., McGonigle A. J. S., Oppenheimer, C., and Allen, A. G. ( 2006 ). A reassessment of current volcanic emissions from the Central American arc with specific examples fro m Nicaragua Journal of Volcanology and Geothermal Research 149 : 297 311.

PAGE 122

111 Mauri, G., Williams Jones, G., and Saracco, G. (2010). Depth determinations of shallow hydrothermal systems by self potential and multi scale wavelet tomography. Journal of Volcanolo gy and Geothermal Research 191 : 233 244. McBirney, A. R. ( 1956 ). The Nicaraguan volcano Masaya and its caldera Eos (Transactions, American Geophysical Union) 37 : 83 96. McPhee, D. K., Jachens, R. C., and Wentworth, C. M. (2004). Crustal structure across the San Andreas Fault at the SAFOD site from potential field and geologic studies. Geophysical Research Lett ers 31 : L12S03. Mheust, Y., and Schmittbuhl, J. (2001). Geometrical heterogeneities and permeability anisotropy of rough fractures. Journal of Geop hysical Res earch 106 : 2089 2102. Mtaxian, J., Lesage, P., and Dorel, J. (1997). Permanent tremor of Masaya volcano, Nicaragua: Wave field analysis and source location. Journal of Geophysical Research, 102:22529 22545. Michel, S., and Zlotnicki, J. (1998). Self Potential and Magnetic Surveying of La Fournaise Volcano (Runion Island): Correlations with Faulting, Fluid Circulation and Eruption. Journal of Geophysical Research, 103: 17845 17857. Nakada, S., Nagai, M., Kaneko, T., Nozawa, A., Suzuki Kamata, K., Nakada, S., and Druitt, T. ( 2005 ). Chronology and products of the 2000 eruption of Miyakejima Volcano, Japan Bulletin of Volcanology 67 : 205 218. Newhall C. G., Albano, S. E., Matsumoto, N., and Sandoval, T. (2001) Roles of groundwater in volcanic unr est. Journal of the G eological Society of the Philippines 56: 69 84. Neuville, A., R. Toussaint, and Schmittbuhl, J. (2006). Hydro thermal flow in a rough fracture. In EHDRA Scientific Conference Soultz sous Forts, France. Nishida, Y., Matsushima, N., Go to, A., Nakayama, Y., Oyamada, A., Utsugi, M., and Oshima H. (1996) Self potential studies in volcanic areas (3) Miyake jima, Esan and Usu J ournal of the Fac ulty of Sci e nce, Hokkaido Univ ersity, Ser ies 7 (Geophysics) 10 : 63 77. Notsu, K., Mori, T., Do V ale, S. C., Kagi, H, and Ito, T. (2006) Monitoring quiescent volcanoes by diffuse CO 2 degassing; case study of Mt. Fuji, Japan Pure and App lied Geophys ics 163 : 825 835. Okada, H., Nishimura, Y., Miyamachi, H., Mori, H., and Ishihara, K. (1990). Geophysi cal significance of the 1988 1989 explosive eruptions of Mt. Tokachi, Hokkaido, Japan. Bulletin of the Volcanological Society of Japan 35 : 175 203.

PAGE 123

112 Oldenburg, C. M., and Pruess, K. (1994) Numerical simulation of coupled flow and transport with TOUGH2: A v erification study. LBL 35273 Lawrence Berkeley National Labo ratory, Berkeley, CA. Overbeek, J. T. G. (1952). Electrochemistry of the double layer I n Colloid Science volume 1: irreversible systems, Kruyt, H. R., editor, p ages 115 193 Elsevier, Amsterdam. Pearson, S. C. P., Connor, C. B. and Sanford W. (2008) Rapid response of a hydrologic system to volcanic activity: Masaya volcano, Nicaragua Geology, 36 : 951 954. Pruess, K. (1991). TOUGH2: A general purpose numerical simulator for multiphase fluid and heat flow. LBL 29400 Lawrence Berkeley National Laboratory, Berkeley, CA. Pruess, K., Oldenburg, C., and Moridis, G. ( 1999 ). TOUGH2 user's guide, version 2.0. LBNL 43134 Lawrence Berkeley National Laboratory, Berkeley, CA. Pruess, K. ( 2004 ). The TOUGH c odes A family of simulation tools for multiphase flow and transport processes in permeable media Vadose Zone Journal 3 : 738 746. Prunty, L., and Bell, J. (2005). Soil temperature change over time during infiltration Soil Science Society of America Journa l 69 : 766 775. Radhakrishna Murthy, I. V., Sudhakar, K. S., and Rama Rao P. (2005) A new method of interpreting self potential anomalies of two dimensional inclined sheets Comput ational Geosci ences 31 : 661 665. Reid, M. E. (2004). Massive collapse of vo lcano edifices triggered by hydrothermal pressurization Geology 32: 373 376 Richter, G., Wassermann, J., Zimmer, M., and Ohrnberger, M. ( 2004 ). Correlation of seismic activity and fumarole temperature at the Mt. Merapi volcano (Indonesia) in 2000 Journ al of Volcanology and Geothermal Research 135 : 331 342. Rogie, J. D., Kerrick, D. M., Sorey, M. L., Chiodini, G., and Galloway D. L. (2001) Dynamics of carbon dioxide emission at Mammoth Mountain, California Earth and Planet ary Sci ences Let ters 188 : 535 541. Rymer, H., van Wyk de Vries, B., Stix, J., and Williams Jones, G. ( 1998 ). Pit crater structure and processes governing persistent activity at Masaya Volcano, Nicaragua Bulletin of Volcanology 59 : 345 355 Salazar, J. M. L., et al. (2002). Precursory diffuse carbon dioxide degassing signature related to a 5.1 magnitude earthquake in El Salvador, Central America. Earth and Planetary Sciences Letters, 205 : 81 89.

PAGE 124

113 Saar, M. O., and Manga M. (1999) Permeability porosity relationship in vesicular basalts Geophys ical Res earch Let ters 26 : 111 114. Sasai, Y., Zlotnicki, J., Nishida, Y., Yvetot, P., Morat, P., Murakami, H., Tanaka, Y., Ishikawa, Y., Koyama, S., and Sekiguchi W. (1997) Electromagnetic Monitoring of Miyake jima Volcano, Izu Bonin Arc, Japan: A Preliminary Report J ournal of Geomag netism and Geoelectr icity 49 : 1293 1316. Scott, K. M, Vallance, J. W., Kerle, N., Macas, J. L., Strauch, W., and Devoli G. (2005). Catastrophic precipitation triggered lahar at Casita volcano, Nicaragua: occurrence, bulking and transformation Earth Surface Processes and Landforms 30: 59 75. Shaw, A. M., Hilton, D. R., Fischer, T. P., Walker, J. A., and Alvarado, G. E. ( 2003 ). Contrasting He C relationships in Nicaragua and Costa Rica: Insights into C cycling through subduction zones Earth and Planetary Science Letters 214 : 499 513 Shibata, T., and Akita, F. (2001) Precursory changes in well water level prior to the March, 2000 eruption of Usu Volcano, Japan Geophys ical Res earch Let ters 28 : 1799 1802. Sparks, R. S J. (2003). Forecasting volcanic eruptions. Earth and Planetary Science Letters 210:1 15. St Amand K., Beaulieu A Stix J Gaonac'h H Lovejoy S Hebert R (1998) SO 2 flux, CO 2 and radon degassing at Masaya Caldera, Nicaragua. Program with Abs tracts Geological Association of Canada Joint Annual Meeting 23:178 Stamatakos, J. A., Connor, C. B., and Martin, R. H. (1997). Quaternary Basin Evolution and Basaltic Volcanism of Crater Flat, Nevada, from Detailed Ground Magnetic Surveys of the Lit tle Cones. Jou rnal of Geol ogy 105 : 319 330. Stevenson, D. S. (1993). Physical models of fumarolic flow. Journal of Volcanology and Geothermal Research 57:139 156. Stiles, A. G., and Felsing, W. A. (1926). The heat of solution of sulfur dioxide. Journal of the American Chemical Society 48 : 1543 1547. Stix, J. (2007). Stability and instability of quiescently active volcanoes; the case of Masaya, Nicaragua. Geology, 35 : 535 538. Stoiber, R. E., Rose, W. I., Lange, I. M., and Birnie, R. W., ( 1975 ). The cooling of Izalco Volcano (El Salvador). Geologisches Jahrbuch 13 : 193 205

PAGE 125

114 Stoiber, R. E., Williams, S. N., Huebert, B. J., Sato, M., Matsuo, S., and King, C. ( 1986 ). Sulfur and halogen gases at Masaya Caldera complex, Nicaragua; total fl ux and variations with t ime Journal of Geophysical Research 91 : 12 215 12231 Sudo, Y., and Hurst, A. W. ( 1998 ). Temperature changes at depths to 150 metres near the active crater of Aso Volcano; Preliminary analysis of seasonal and volcanic effect s. Journal of Volcanology and Ge othermal Research 81 : 159 172. Tanguy, J. C. (1994). The 1902 1905 eruptions of Montagne Pele, Martinique: anatomy and retrospection Journal of Volcanology and Geothermal Research 60:87 107. Tedesco, D., Toutain, J P., Allard, P., and Losno, R. ( 1991 ). Chemical variations in fumarolic gases at Vulcano Island (southern Italy); Seasonal and volcanic effects Journal of Volcanology and Geothermal Research 45 : 325 334. Tenorio, V., Saballos, A., Morales, A., Herrera, M., Alvarez, J., Mayorga, E., Muoz, A., Prez, P., Talavera, E. Acosta, A., Bodn M., Guzmn, C., Segura, F., and Zamora, J. (2010). INETER Boletines Mensuales y Anuales Volcn Masaya, 2006 2010 http://www.ineter.gob.ni/geofisica/sis/bolsis/bolsis.html accessed 28 June 2010. Todesco, M. (1997). Origin of fumarolic fluids at Vulcano (Italy); insights from isotope data and numerical modeling of hydrothermal circulation Journal of Volcanol ogy and Geotherm al Res earch 79 : 63 85. Todesco, M., Chiodini, G., and Macedonio, G. ( 2003 ). Monitorin g and modelling hydrothermal fluid emission at L a S olfatara ( P hlegrean F ields, I taly); an interdisciplinary approach to the study of diffuse degassing Journal of Volcanology and Geothermal Research 125 : 57 79. Todesco, M., Rutqvist, J., Chiodini, G., Prue ss, K., and Oldenburg, C.M ( 2004 ). Modeling of recent volcanic episodes at P hlegrean F ields ( I taly) : Geochemical variations and ground deformation Geothermics 33 : 531 547. Todesco, M. and Berrino, G. (2005). Modeling hydrothermal fluid circulation and gr avity signals at the Phlegraean Fields caldera. Earth and Planetary Science Letters 240:328 338. Todesco, M., Chiodini, G., and Berrino, G. (2006). Modeling of gas composition and gravity signals at the phlegrean fields caldera. Proceedings, TOUGH Symposi um 2006 Lawrence Berkeley National Laboratory, Berkeley, California, May 15 17. Todesco, M. (2008). Hydrothermal fluid circulation and its effect on caldera unrest. In Gottsmann, J., and Marti, J., editors, Caldera Volcanism 10, Analysis, Modelling and Re sponse p ages 393 416 Elsevier Science

PAGE 126

115 Toutain, J. P., Sortino, F., Baubron J. Richon, P., Surono, S. Sumarti, and Nonell A. (2009) Structure and CO 2 budget of Merapi volcano during inter eruptive periods Bull etin of Volcanol ogy 71 : 815 826 Van Wyk de Vries B (1993). Tectonics and magma evolution of Nicaraguan volcanic systems PhD dissertation Open University, Milton Keynes, UK. van Padang, M. ( 1963 ). The temperatures in the crater region of some I ndonesian volcanoes before the eruption Bul letin of Volcanology 26 : 319 336. Varley, N. R. and Armienta, M. A. ( 2001 ). The absence of diffuse degassing at P opocatpetl volcano, M exico Chemical Geology 177 : 157 173. Voss, C. I., and Souza, W. R. (1987). Variable density flow and solute transport si m ulations of regional aquifers containing a narrow freshwater saltwater transition zone. Water Resources Resea r ch, 23:1851 1866. Walker, J. A., Williams, S. N., Kalamarides, R. I., and Feigenson M. D. (1993) Shallow open system evolution of basaltic magm a beneath asubduction zone volcano: the Masaya Caldera Complex, Nicaragua J ournal of Volcanol ogy and Geotherm al Res earch 56 : 379 400. Wang, G., and Sassa, K. ( 2003 ). Pore pressure generation and movement of rainfall induced landslides; Effects of grain si ze and fine particle content Engineering Geology 69 : 109 125. Wehrmann, H., Bonadonna, C., Freundt, A., Houghton, B. F., and Kutterolf, S. (2006). Fontana Tephra: a basaltic Plinian eruption in Nicaragua I n Volcanic Hazards in Central America, Rose, W. I ., Bluth, G. J. S., Carr, M. J., Ewert, J., Patino, L. C., and Vallance., J. W., editors, pages 209 223 Geol ogical Soc iety of Amer ica Special Paper 412. Werner, C., Christenson, B. W., Hagerty, M. and Britten K. (2006) Variability of volcanic gas emiss ions during a crater lake heating cycle at Ruapehu Volcano, New Zealand J ournal of Volcanol ogy and Geotherm al Res earch 154 : 291 302 White, F. M. (198 4 ). Heat transfer. Addison Wesley, Reading, MA Williams, S. N. ( 1983 a ). Plinian airfall deposits of basa ltic composition Geology 11 : 211 214. Williams S. N. (1983b). Geology and Eruptive Mechanisms of the Masaya Caldera Complex, Nicaragua Ph .D. dissertation Dartmouth College, Hanover, N H, USA.

PAGE 127

116 Williams Jones, G., Rymer, H., and Rothery, D.A. ( 2003 ). Grav ity changes and passive SO 2 degassing at the M a saya caldera complex, Nicaragua Journal of Volcanology and Geothermal Research 123 : 137 160. Wohletz, K., and Heiken, G. (1992). Volcanology and geothermal energy Berkeley: University of California Press, pa ges 95 118. Yamashina, K., and Matsushima, T. ( 1999 ). Ground temperature change observed at Unzen Volcano associated with the 1990 1995 eruption Journal of Volcanology and Geothermal Research 89 : 65 71 Zimmer, M., and Erzinger, M. ( 2003 ). Continuous H 2 O CO 2 Rn 222 and temperature measurements on Merapi Volcano, Indonesia Journal of Volcanology and Geothermal Research 125 : 25 38. Zablocki, C. J. (1976). Mapping thermal anomalies on an active volcano by the self potential method, Kilauea, Hawaii Procee dings, 2nd UN Symposium of the development and use of geothermal resources, San Francisco, May 1975 1299 1309 Zhao, C., Hobbs, B., Mhlhaus, H., Ord, A., and Lin, G. (2003). Convective instability of 3 d fluid saturated geological fault zones heated from below Geophys ical J ournal In ternational 155 : 213 220. Zlotnicki, J., and Le Mouel J. L. (1990) Possible electrokinetic origin of large magnetic variations at La Fournaise volcano Nature, 343 : 633 636. Zlotnicki, J., Sasai, Y., Toutain, J., Villacorte, E., Bernard, A., Sabit, J., Gordon, J., Corpuz, E., Harada, M., Punongbayan, J., Hase, H., and Nagao, T. ( 2009 ). Combined electromagnetic, geoc hemical and thermal surveys of T aal volcano ( P hilippines) during the period 2005 2006 Bulletin of Volcanology 7 1 : 29 47.

PAGE 128

117 APPENDIX A Finite difference model of 2 D heat transfer

PAGE 129

118 % Matlab program to calculate 2D heat transfer by conduction and convection. Top and bottom boundaries are fixed temperature, right boundary is no flow and left boundary is convective. close all % Start with a clear workspace clear all Tnew=0; % Set the starting temperature at 0 nptx=25; % Set the number of points in the horizontal direction nptz=25; % Set the n umber of points in the vertical direction len=300; % Set the maximum number of iterations dx=(nptx 1)*0.02; % Calculate the grid spacing % Create a square mesh for the model domain T=200*ones(nptz,nptx,len+1); % Set the boundary conditions T(1,:, :)=10; % Fixed temperature at the top to represent air T(nptz,:,:)=800; % Fixed temperature at the bottom to represent heat source for i=2:nptz % Fixed heat flux condition T(i,:,:)=(T(nptz,:,:) T(1,:,:))/(nptz 1)+T(i 1,:,:); end % Start the iterati ons to calculate solution for j=1:len for q=1:1000000 % Calculate convective boundary % Call subroutine convect to calculate boundary condition [T(:,nptx,j),alpha,k,h,Tave]=convect(T(:,nptx 1,j),dx,T(nptz,1,1)); Bi=h*dx/k; % Calc how much of temp is transferred to adjacent cell T(1,nptx,j)=(1/(1+Bi))*[(T(1,1,j)+T(2,nptx,j))/2+(Tave*Bi)]; % Calculate new temperature from adjacent cells % Explicit, finite difference approximation for wall rock conduction [T,dt]= transient (T,nptx,nptz,j,al pha,k,h,dx); %Test for convergence within a cell diff=abs(Tnew T); if diff<1 break % If the temp difference is less than 1C, calculate the next cell e nd Tnew=T ; % Otherwise update the temp conditions and try agai n end if j>=2

PAGE 130

119 Appendix A (continued) difft=abs(T(:,:,j) T(:,:,j 1)); % Calculate temp diff between iterations if difft<1 % If less than 1C, simulation converged break end end end % Plot the results x=1:nptx; x=dx*x; figure [C,h] = contour((max(x) x),(x (1*dx)),T(:,:,j)); set(gca,'YDir','rev') xlabel('Distance from conduit (m)') ylabel('Depth from surface (m)') clabel(C,h); -----------------------------------------------------------------------------------------------------------% Subroutine convect.m % Subroutine called from conduct.m. Required inputs are temperature (Tw), length of cell (dx), and temperature at base (Ten) function [T,alpha,k,h,Tavg]=convect(Tw,dx,Ten) load steam.mat ASCII % Some constants, held in an exter nal file a=0.1; % Radius of conduit through which hot gas flows len=length(Tw); % Calculate how many cells there are in the boundary % Renumber the cells because they are upside down for k=1:len Twr(k)=Tw(len k+1); end Tw=Twr; %Calculate new cell values for m=1:len Tex=(Tw(m)+Ten)/2; % Exit temperature is average of previous temperature and input Tavg=exp((log(Tex)+log(Ten))/2); % Calculate average temperature Texold=0; % Make sure exit temperature of each c ell is 0 for 1 st calc for n=1:100000 % Call subroutine vars to find the properties of the fluid [rho,cp,k,viscw,Pr,alpha]=vars(Tw(m),steam); [rho,cp,k,visc,Pr,alpha]=vars(Tavg,steam); v=10; % Velocity of the fluid

PAGE 131

120 Append ix A (continued) k=3; % Thermal conductivity, set at 3 to match COMSOL model % Call subroutine nusselt to calculate convection heat transfer coefficient [h]=nusselt(a,rho,v,visc,viscw,Pr,k); % Calculate new values for the te mperature in each cell dif=Ten Tw(m); expt= 2*h*dx/(v*rho*a*cp); Tex=dif*exp(expt)+Tw(m); Tavg=exp((log(Tex)+log(Ten))/2); % See if cells in convective boundary have converged to a solution conv=abs(Tex Te xold); if conv<=0.001 break % If converged, exit if loop end Texold=Tex; % If not converged, update temperatures and start cycle again end Ten=Tex; % Update temperatures for next iteration T(len m +1)=Ten; end ---------------------------------------------------------------------------------------------------------% Subroutine vars.m % Subroutine to find the relevant properties of water/steam based on the temperature, called from convect.m funct ion [rho,cp,k,visc,Pr,alpha]=vars(T,steam); if T<325 i=1; elseif (T>=325) & (T<375) i=2; elseif (T>=375) & (T<425) i=3; elseif (T>=425) & (T<475) i=4; elseif (T>=475) & (T<525) i=5; elseif (T>=525) & (T<575) i=6; elseif (T>=575) & (T<650) i=7;

PAGE 132

121 Appendix A (continued) elseif (T>=650) & (T<750) i=8; elseif (T>=750) & (T<850) i=9; elseif (T>=850) & (T<950) i=10; elseif (T>=950) i=11; end rho=steam(i,2); % Density cp=steam(i,3); % Speci fic heat k=steam(i,4); % Thermal conductivity visc=steam(i,5); % Dynamic viscosity v =steam(i,6); % Kinematic viscosity Pr=steam(i,7); % Prandtl number alpha=steam(i,8); % Thermal diffusivity ---------------------------------------------------------------------------------------------------------% Subroutine nusselt.m % Subroutine to calculate the convection heat transfer coefficient function [h]=nusselt(a,rho,v,visc,viscw,Pr,k) D=2*a; % Diameter of pipe Re=rho*v*D/visc; er if Re<2200 Nu=8.66; % Nusselt number for laminar flow else Nu=0.027*(Re^0.8)*(Pr^(1/3))*((visc/viscw)^0.14); % Nusselt number for turbulent flow end h=Nu*k/D; % Convection heat transfer coefficient ---------------------------------------------------------------------------------------------------------% Subroutine transient.m % Subroutine to calculate conduction within model function [T,dt]=transient(T,nptx,nptz,j,alpha,k,h,dx) sigma=0.2; % Multiplication factor dt=s igma*(dx)^2/alpha; % Time step for l=2:nptz 1 for m=2:nptx 1 T(l,m,j+1)=sigma*[T(l,m 1,j)+T(l,m+1,j)+T(l 1,m,j)+T(l+1,m,j)] ...

PAGE 133

122 +(1 4*sigma)*T(l,m,j); end end -----------------------------------------------------------------------------------------------------------% File steam.mat, containing properties of water/steam ( White, 198 4 ) Temp Density Specific heat Thermal conductivity Dynamic viscosity Kinematic viscosity Prandtl number Thermal diffusivity K Kg/m 3 J/kgK W/mK PaS m 2 /s m 2 /s 300 0.0253 2041 0.0181 0.91e 5 36.1e 5 1.03 35.1e 5 350 0.258 2037 0.0222 1.12e 5 4.33e 5 1.02 4.22e 5 400 0.555 2000 0.0264 1.32e 5 2.38e 5 1.00 2.38e 5 450 0.491 1968 0.0307 1.52e 5 3.10e 5 0.95 3.17e 5 500 0.441 1977 0.0357 1.73e 5 3.92e 5 0.96 4.09e 5 550 0.401 1994 0.0411 1 .93e 5 4.82e 5 0.94 5.15e 5 600 0.367 2022 0.0464 2.13e 5 5.82e 5 0.93 6.25e 5 700 0.314 2083 0.0572 2.54e 5 8.09e 5 0.93 8.74e 5 800 0.275 2148 0.0686 2.95e 5 10.7e 5 0.92 11.6e 5 900 0.244 2217 0.078 3.36e 5 13.7e 5 0.95 14e.4e 5 1000 0.220 2288 0.0 87 3.76e 5 17.1e 5 0.99 17.3e 5 Figure A 1 Early numerical simulations, of magma conduit and wall rock next to it. Convective left boundary simulates flow along the conduit to the surface. M ATLAB results using my finite difference code (a) show excellent agreement with results using COMSOL multiphysics ( b), although the horizontal extent of the models is slightly diff erent.

PAGE 134

123 APPENDIX B Fumarole temperature data plot ted be year

PAGE 135

124 Figure B 1. Fumarole temperature data in 2006.

PAGE 136

125 Figure B 2 Fumarole temperature data in 2007.

PAGE 137

126 Figure B 3. Fumarole tempe rature data in 2008.

PAGE 138

127 Figure B 4. Fumarole temperature data in 2009.

PAGE 139

128 APPENDIX C Spectrograms

PAGE 140

129 Figure C 1. Top: Spectrogram from thermocouple at 33 cm depth. Bottom: Fumarole temperature time series at 33 cm depth.

PAGE 141

130 Figure C 2. Top: Spectrogram fro m thermocouple at 6 3 cm depth. Bottom: Fumarole temperature time series at 6 3 cm depth

PAGE 142

131 Figure C 3. Top: Spectrogram from thermocouple at 95 cm depth. Bottom: Fumarole temperature time series at 95 cm depth

PAGE 143

132 Figure C 4. Top: Spectrogram from thermocouple at 150 cm depth. Bottom: Fumarole temperature time series at 150 cm depth.

PAGE 144

133 Figure C 5. Top: Spectrogram from thermocouple 7 m away Bottom: Fumarole temperature time series 7 m away.

PAGE 145

134 Figure C 6. Top: Spectrogram of air temperature time se ries Bottom: Air temperature time series.

PAGE 146

ABOUT THE AUTHOR Sophie Pearson was born in Emsworth, Hampshire, England. S he received her MGeophys in Geophysical Sciences from Leeds University, UK i n 2005 which included a year studying at Penn State Univer sity, USA. A week after obtaining her degree she joined the Department of Geology at the Unversity of South Florida (USF) as a PhD student under the supervision of Dr. Charles Connor. She received a Presidential Fellowship from USF to support her throughou t her PhD, which included full funding for five years. As well as working at USF, Sophie participated in conferences and fieldwork in USA, Ecuador, Nicaragua, Taiwan, Japan, and Iceland. Her work has been published in international journals and in 2009 she received the Outstanding Service Award from the Department of Geology She has now accepted a permanent position as a Geothermal Geophysicist with GNS Science in New Zealand.


Download Options

Choose Size
Choose file type
Cite this item close


Cras ut cursus ante, a fringilla nunc. Mauris lorem nunc, cursus sit amet enim ac, vehicula vestibulum mi. Mauris viverra nisl vel enim faucibus porta. Praesent sit amet ornare diam, non finibus nulla.


Cras efficitur magna et sapien varius, luctus ullamcorper dolor convallis. Orci varius natoque penatibus et magnis dis parturient montes, nascetur ridiculus mus. Fusce sit amet justo ut erat laoreet congue sed a ante.


Phasellus ornare in augue eu imperdiet. Donec malesuada sapien ante, at vehicula orci tempor molestie. Proin vitae urna elit. Pellentesque vitae nisi et diam euismod malesuada aliquet non erat.


Nunc fringilla dolor ut dictum placerat. Proin ac neque rutrum, consectetur ligula id, laoreet ligula. Nulla lorem massa, consectetur vitae consequat in, lobortis at dolor. Nunc sed leo odio.