Developing a new passive diffusion sampling array to detect helium anomalies associated with volcanic unrest

Helium (He) concentration and 3He/4He anomalies in soil gas and spring water are potentially powerful tools for investigating hydrothermal circulation associated with volcanism and could perhaps serve as part of a hazards warning system. However, in opera

DEVELOPING A NEW PASSIVE DIFFUSION SAMPLING ARRAY TO DETECT HELIUM ANOMALIES ASSOCIATED WITH VOLCANIC UNREST by Brittany Elise Dame A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science in Geology Department of Geology and Geophysics The University of Utah $XJXVW 2014 Copyright © Brittany Elise Dame 2014 All Rights Reserved The University of Utah Graduate School STATEMENT OF THESIS APPROVAL The following faculty members served as the supervisory committee chair and Brittany Elise Dame members for the thesis of_________________________________________________. Dates at right indicate the members' approval of the thesis. D. Kip Solomon _________________________________________, Chair 4/16/2014 ________________ Date Approved John R. Bowman ________________________________________, Member 4/16/2014 ________________ Date Approved William C. Evans ________________________________________, Member 4/16/2014 ________________ Date Approved John M. Bartley The thesis has also been approved by__________________________________ Geology and Geophysics Chair of the Department of____________________________________ and by David B. Kieda, Dean of The Graduate School. ABSTRACT Helium (He) concentration and 3He/4He anomalies in soil gas and spring water are potentially powerful tools for investigating hydrothermal circulation associated with volcanism and could perhaps serve as part of a hazards warning system. However, in operational practice, He and other gases are often sampled only after volcanic unrest is detected by other means. To fully investigate the potential utility of He data in the context of volcanic unrest, samples would be needed days to weeks before any detectable seismic energy or geodetic deformation. Accurately recording precursory He anomalies using traditional collection methods would require high-cost sampling surveys with daily or weekly resolution. A new passive diffusion sampler array, intended to be collected after the onset of unrest, has been developed and tested as a relatively low-cost alternative. The samplers, each with a distinct equilibration time, passively record He concentration and isotope ratio in springs and soil gas. Helium gas diffuses through a semipermeable membrane into variably-sized gas reservoirs. Once collected and analyzed, the He concentrations in the samplers are used to deconvolve the time history of the He concentration and the 3He/4He ratio at the collection site. An array consisting of three samplers is sufficient to deconvolve both the magnitude and the timing of a step change in in situ concentration if the array is collected within 100 hours of the pulse. Volumetrically larger arrays may allow for longer time elapses between the end of a pulse and the array collection. As equilibration of each sampler is a diffusion process, the He effective diffusion coefficient for the semipermeable membrane was determined in aqueous and gaseous phases under varying temperatures. The diffusion model fits laboratory tests accurately and reliably predicts the timing and variation of dissolved He gas. In initial field studies, the array has captured an in situ He change and was used to deconvolve the He history. Passive diffusion sampler arrays appear to be an accurate and affordable alternative for determining He anomalies associated with volcanic unrest. iv TABLE OF CONTENTS ABSTRACT………………………………………………………………………. iii LIST OF FIGURES………………………………………………………………. vii ACKNOWLEDGEMENTS………………………………………………………. viii 1. INTRODUCTION…………………………………………………………….. . 1 2. DESCRIPTION AND OPERATION……………...………………………….. ... 4 3. THEORY………………………………………………………….…………… 7 4. TESTING AND RESULTS……………………………………………………. 14 4.1 Equilibration Constant and Calibration…………………….…………... 14 4.2 Lab/Model Validation……………………………………….…………. 19 5. MODEL UNCERTAINTY…………………………………………………….. 23 5.1 Model Experiments……………………………………………………... 23 5.2 Temperature Dependence of Helium Diffusion…………………………. 25 5.3 Environmental Effects on Diffusion……………………………………. 34 6. FIELD RESULTS……………………………………………………..………… 39 7. SUMMARY AND CONCLUSIONS…………………………………………… 46 Appendices A. ONE-STEP INVERSE MODEL MATLAB CODE………………………… 48 B. ONE-STEP FORWARD MODEL MATLAB CODE……………………...… 51 C. ONE-STEP FRECHET DERIVATIVE MATLAB CODE……………….…… 52 D. TWO-STEP INVERSE MODEL MATLAB CODE……………………...…….53 E. TWO-STEP FORWARD MODEL MATLAB CODE……………………..… 56 F. TWO-STEP FRECHET DERIVATIVE MATLAB CODE…………………… 57 REFERENCES………………………………………………………….…………… 58 vi LIST OF FIGURES 2.1. Sampler schematic………..…………………………………………………… 5 3.1 One-step model………………………………………………………………... 9 3.2 Two-step model………………………………………………………………... 9 4.1 Gas chromatograph inlet manifold…………………………………………….. 17 4.2 RSS error for the equilibration constant as a function of equilibration………… 20 4.3 Calibrated equilibration constants and curves for a sampling array……………. 21 4.4 One-step inverse model applied to experimental data………………………….. 22 5.1 Model error for varying pulse durations and time elapsed since the pulse………24 5.2 Model error as a function of equilibration constant variation……………………26 5.3 Copper (Cu) tubes fitted with semipermeable silicone membranes…………… 28 5.4 Temperature dependence of He diffusion……………………………………… 30 5.5 Equilibration constant (Ec) as a function of temperature…………………………32 5.6 Effect of temperature on two-step inverse model results……………………… 33 5.7 Environmental effects on the equilibration constant……………………………..38 6.1 R/Ra ratios for 2011 deployment……………………………..…..…..…..…..…. 40 6.2 One-step inverse modeling of Horseshoe Lake data………..…..…..…..…..…. 43 6.3 Seismic activity (M>0.5) in the Mammoth Mountain area …..…..…..…..…..….44 ACKNOWLEDGEMENTS Special thanks to my advisor, Dr. Kip Solomon, for letting me tackle this project on my own terms. With his guidance, I was allowed to take near full control of the project and gain insight into my strengths and weaknesses as a researcher and scientist. Also, thank you to Wil Mace and Alan Rigby for being available at any moment at any time of the day. Their lab expertise and willingness to help were invaluable. To my future husband, Mike, and my mother, Debby, thank you! Thank you for boundless love, support, and humor. Your jokes fix almost every problem. Almost. 1. INTRODUCTION The noble gas content of the mantle, crust, and atmosphere has been used extensively to study processes such as plate tectonics [Clarke et al., 1969], ocean circulation [Lupton and Craig, 1981], groundwater flow [Andrews, 1985; Torgersen and Clarke, 1985; Solomon, 2000], and hydrothermal circulation in volcanic terrain [Lupton et al., 1977; Craig et al., 1978]. Among the volatile species, helium (He) is particularly useful because of its primordial 3He isotope derived from the Earth's mantle. Most terrestrial 3He/4He variations can be explained by mixing of the three major endmembers: atmospheric He with 3He/4He = 1.4 x 10-6, radiogenic He with 3He/4He ≈10-8 to 10-7, and mantle He with 3He/4He ≈ 10-5 [Lupton, 1983]. Atmospheric He is nearly constant because of the rapid mixing time of the atmosphere and He's extended residence time. This uniform isotopic composition allows laboratories to use the atmospheric isotopic composition as a standard. Thus, He isotope variations are expressed relative to the 3He/4He ratio in air. In accordance with this nomenclature, units R/RA are used throughout this manuscript, where R = 3He/4He of a sample and RA = (3He/4He)AIR. Radiogenic He is produced by radioactive decay of uranium and thorium series elements. The 4He isotope is produced by alpha particle decay, while the 3He isotope is produced by lithium spallation. Radiogenic He is most abundant in the continental crust where U and Th are most concentrated and the isotopic ratio is dependent on the relative 2 concentrations of Li, U, and Th. Typical crustal rocks produce He with R/RA ≈ 0.01 - 0.1 [Morrison and Pine, 1955]. Mantle He is enriched with 3He that is assumed to be remnant from the time of Earth's formation [Clarke et al., 1969; Mamyrin et al., 1969; Lupton and Craig, 1975]. Helium-3 enrichments (R/RA = 5 - 30) are found in many tectonic environments including midocean ridges, subduction zones, hot spots and oceanic islands, and continental geothermal areas. As such, He concentration and isotopic ratio are powerful tools for investigating magmatic processes. Passively degassing, high temperature fumaroles and gaseous springs provide the opportunity to sample volatiles released directly from volcanic magma bodies. The 3 He/4He ratio has been used to differentiate He contribution between mantle and continental derived sources [Craig et al., 1978; Marty et al., 1989]. Variations in the He isotopic ratio may be useful to evaluate the state of volcanic activity and estimate magmatic intrusion mode, rate, and duration [Sorey et al., 1993; Sano et al., 1991]. Temporal variations in He concentration and 3He/4He may be of practical use in forecasting volcanic eruptions as these variations in hydrothermal gases provide evidence that new magma is entering the volcanic system [Sano et al., 1988]. In a study by Carapezza [2004], He concentration increased significantly in thermal waters prior to eruptive events at Stromboli in the absence of geophysical precursors. Increased diffusive He emission and isotope ratio was also documented by Padron et al. [2013] a month prior to a submarine eruption on El Hierro. This major He release also preceded increases in seismic energy associated with the volcanic unrest. These anomalies were only detected because sampling increased from monthly to daily and weekly intervals. 3 To investigate fully the potential utility of He data in forecasting volcanic unrest, samples may be needed days to weeks before any detectable seismic energy or geodetic deformation. Accurately recording precursory He anomalies using traditional collection methods requires high-cost sampling surveys with daily or weekly resolution. In an effort to efficiently and affordably evaluate He anomalies, a new passive diffusion sampler array, intended to be collected after the onset of unrest, has been developed and tested. Passive diffusion samplers have been used for over three decades to monitor work place chemical exposure [Brown et al., 1981], environmental pollutants [Kot et al., 2000; Stuer-Lauridsen, 2005], and dissolved gases in water [Sanford et al., 1996; Gardner and Solomon, 2009]. These methods employ equilibrium sampling for passive, in situ gas sampling. Based on passive sampling techniques, our sampling array was designed to passively record He concentration and isotope ratio in spring water and soil gas. At the sampling site, He gas diffuses through a semipermeable membrane into three variablysized gas reservoirs. Each of these reservoirs, or samplers, has a distinct equilibration time. As the He concentration or isotope ratio changes through time, each of the samplers records that change at a different rate. After volcanic unrest is detected by other (seismic or geodetic deformation) means, the sampling array is collected. Once analyzed, the He concentrations in the samplers are used to deconvolve the time history of the He concentration and the 3He/4He ratio at the collection site. The array is capable of deconvolving both the magnitude and the timing of a step change in in situ concentration. 2. DESCRIPTION AND OPERATION Figure 2.1 shows a schematic of the three samplers that make up a sampling array. Helium and other gasses diffuse across a semipermeable silicone membrane into the sampling reservoir. The gases are passively recorded by each sampler. The reservoir size (volume (V)), area of the membrane (A), thickness of the membrane (L), and the effective diffusion coefficient for He across the membrane (D) determine the equilibration time of each sampler. The small sampler is approximately 16 centimeters (cm) long and less than six cm wide including the pinch clamp. The medium and large samplers are ~27 and ~19 cm long, respectively. Similar to the passive diffusion samplers described by Gardner and Solomon [2009], the major sampling features include the sample volume where the gas is captured and stored and the semipermeable silicone membrane. The gas exchange area is a ~2.5 cm strip of 0.08 cm thick silicone rubber tubing stretched over 8 mm diameter stainless steel tubing fitted with fine stainless steel mesh. The mesh provides channels for gas movement along the tubing to a gas inlet port. The internal volume of the small, medium, and large samplers is ~3, ~14.5, ~68 cm3, respectively. The volume is gas tight and can be deployed in water (at depths up to 100 m), soil, or fumarole vents. Upon deployment, the array is initially equilibrated with ambient air. The samplers are submerged in the fluid or buried in an acrylonitrile-butadiene-styrene (ABS) or polyvinyl chloride (PVC) casing open at the bottom. A gas or water sample is taken to 5 Figure 2.1 Sampler schematic 6 determine the background He concentration and isotope ratio. The samplers begin to equilibrate with the surrounding gases in the medium immediately. The array is designed to be deployed until volcanic unrest is detected by seismic energy or geodetic deformation and to perform optimally should be recovered within 100 hours of the He pulse; recovery after longer periods of time is possible at the expense of resolution. During collection, the samplers are brought to the land surface where the pinch clamp is immediately closed and the gas is trapped in the nickel sampling tip. An additional water or gas sample should be collected at this time. The concentrations in each sampler, as well as those taken at deployment and collection, are used to deconvolve the history of the site concentration using inverse modeling. Prolonged exposure to harsh environmental conditions may affect the diffusive properties of the membrane. The effective diffusion coefficient is one factor that determines the equilibration time of each sampler and knowing this equilibration time is a necessity for the inverse modeling. Each sampler should be recalibrated after the deployment to determine its unique equilibration constant for He. 3. THEORY Equilibrium headspace sampling allows for passive in situ sampling of dissolved gasses in water. The in situ partial pressure of He is given by Henry's Law: 𝑃𝐻𝑒 = 𝐾𝐻𝑒 𝐶𝑤 (2.1) where KHe is Henry's coefficient for He and is dependent on temperature and Cw is the concentration of He in the water. The partial pressure inside the sampler (Psampler) is derived from the diffusive mass flux of gas through a thin membrane into a reservoir, assuming a steady state concentration profile within the membrane [Sanford et al., 1996] , and is given by 𝑑𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 −𝐷𝐴 = 𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 − 𝐾𝐻𝑒 𝐶𝑤 . 𝑑𝑡 𝑉𝐿 (2.2) The solution to equation (2.2) subject to the initial condition where Psampler =Po at time (t) equal zero, also given by Stanford et al. [1996], is 𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 = 𝑃𝑜 exp − 𝐷𝐴 𝐷𝐴 𝑡 + 𝐾𝐻𝑒 𝐶𝑤 1 − exp − 𝑡. 𝑉𝐿 𝑉𝐿 (2.3) For the condition where Po is equal to the gas phase concentration in the atmosphere (Patm), equation (2.3) can be written as 𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 = 𝑃𝑎𝑡𝑚 exp − 𝐷𝐴 𝑡 𝑉𝐿 𝐷𝐴 + 𝐾𝐻𝑒 𝐶𝑤 1 − exp − 𝑡. 𝑉𝐿 The equilibration constant for a sampler is defined by (2.4) 8 𝐸𝑐 = 𝐷𝐴 . 𝐿 (2.5) As the samplers can be used to sample both He in water or soil, the partial pressure of He is substituted for Henry's coefficient and the concentration in water to simplify the expression. Utilizing the equilibration constant and the partial pressure of He in the environment (Penv) equation (2.3) becomes 𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 = 𝑃𝑎𝑡𝑚 exp − 𝐸𝑐 𝐸𝑐 𝑡 + 𝑃𝑒𝑛𝑣 1 − exp − 𝑡. 𝑉 𝑉 (2.6) Two mathematical models were developed to describe hypothetical step changes in the environment. The simplest one-step model describes an increase in the partial pressure of He (Figure 3.1). At some time t1, the dissolved He concentration increases from some background partial pressure P1 to a new partial pressure P2. It stays constant until some later time t2 when the samplers are collected. A second two-step model was also developed to describe the same increase in He partial pressure from P1 to P2 at some time t1; however, in this model the partial pressure of He also decreases to P3 at some time t2 until the samplers are collected at a later time t3 (Figure 3.2). The superposition principle was applied to equation (2.6) to describe the partial pressure in a sampler as a function of the hypothetical one-step model: 𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 = 𝑃𝑎𝑡𝑚 exp − 𝐸𝑐 𝐸𝑐 𝑡2 + 𝑃1 1 − exp − 𝑡2 𝑉 𝑉 𝐸𝑐 + (𝑃2 − 𝑃1 ) 1 − exp− (𝑡2 − 𝑡1 ) . 𝑉 (2.7) This method was also used to create a forward operator that describes the partial pressure in the sampler as a function of the two-step model: 9 Figure 3.1 One-step model showing an increase in He partial pressure (black line) and each sampler's equilibration curve (gray line). The boxes represent the He partial pressure in each sampler at the collection time. Figure 3.2 Two-step model showing an increase in He partial pressure at t1 and slight decrease at t2 (black line). Each sampler's equilibration curve is shown in gray. The boxes represent the He partial pressure in each sampler at the collection time. 10 𝑃𝑠𝑎𝑚𝑝𝑙𝑒𝑟 = 𝑃𝑎𝑡𝑚 exp − 𝐸𝑐 𝐸𝑐 𝑡3 + 𝑃1 1 − exp − 𝑡3 𝑉 𝑉 + (𝑃2 − 𝑃1 ) 1 − exp− + (𝑃3 − 𝑃2 ) 1 − exp − 𝐸𝑐 (𝑡 − 𝑡1 ) 𝑉 3 (2.8) 𝐸𝑐 (𝑡 − 𝑡2 ) 𝑉 3 The objective of the one-step inverse model is to minimize the difference between the measured and calculated partial pressure of He in the samplers by solving equation (2.8) for the timing (t1) and magnitude (P2) of the increase in He partial pressure. The objective of the two-step inverse model follows the same minimization process by utilizing equation (2.8) to estimate the duration (t1 and t2) and magnitude (P2) of the He pulse. The system of nonlinear equations for both models can be written in operator notations as: 𝐴(𝑚) = 𝑑 (2.9) where m is the vector-column containing the model parameters and d is the vectorcolumn of data. The data for both models is the same: 𝑋𝐻𝑒 𝑠 𝑑 = 𝑋𝐻𝑒 𝑚 . 𝑋𝐻𝑒 𝑙 (2.10) The one-step model parameters are given by: 𝑃 𝑚 = 2 . 𝑡1 (2.11) Whereas the two-step model parameters are 𝑃2 𝑚 = 𝑡1 . 𝑡2 (2.12) 11 For each sampler volume and equilibration constant, the nonlinear forward operator (A) in equation (2.9) acts on the model parameters (m) and predicts the pressure of He present in each sampler. To solve this equation for the unknown model parameters (m) in both the forward operators, two iterative, regularized conjugate gradient models were written to find the minimum of the Tikhonov parametric functional [Zhdanov, 2002]: 𝑃𝛼 (𝑚) = 𝜑(𝑚) + 𝛼𝑠(𝑚) (2.13) 𝜑(𝑚) = (𝐴(𝑚) − 𝑑)𝑇 (𝐴(𝑚) − 𝑑) (2.14) where φ(m) is the misfit functional between the predicted data (A(m)) and the observed data (d). The stabilizing functional (s(m)) is introduced as the least-squares difference between the regularized solution and some a priori model (mapr). 𝑇 𝑠(𝑚) = 𝑚 − 𝑚𝑎𝑝𝑟 𝑚 − 𝑚𝑎𝑝𝑟 (2.15) The parametric functional, equation (2.13) , can be rewritten as 𝑃𝛼 (𝑚, 𝑑) = (𝐴(𝑚) − 𝑑)𝑇 (𝐴(𝑚) − 𝑑) 𝑇 (2.16) + 𝛼𝑚 − 𝑚𝑎𝑝𝑟 𝑚 − 𝑚𝑎𝑝𝑟 . The regularized descent method solves for a quasi-solution of the inverse problem by minimizing the parametric functional. 𝑃𝛼 (𝑚, 𝑑) = 𝑚𝑖𝑛. (2.17) To find the gradient direction, we calculate the first variation of the parametric functional. The operator A is a differentiable operator, such that 𝛿𝐴(𝑚) = 𝐹𝑚 𝛿𝑚 (2.18) 12 where Fm is the Frechet derivative of the operator matrix A. The first variation of the parametric functional is given by: 𝛿𝑃𝛼 (𝑚, 𝑑) = 2(𝛿𝑚)𝑇 𝐹𝑚 𝑇 (𝐴(𝑚) − 𝑑) 𝑇 (2.19) + 2𝛼(𝛿𝑚) 𝑚 − 𝑚𝑎𝑝𝑟 . The iterative process is based on the calculation of the regularized steepest descent directions: 𝑚𝑛+1 = 𝑚𝑛 − 𝛿𝑚 = 𝑚𝑛 − 𝑘𝑛𝛼 𝑙̃𝛼 (𝑚𝑛 ). The directions of steepest ascent ( (2.20) ) are selected according to the following algorithm. On the first step, the direction of ascent is determined by the least squares method: 𝑙̃𝛼 (𝑚𝑛 ) = 𝐹𝑚 𝑇 (𝐴(𝑚) − 𝑑) + 𝛼𝑚 − 𝑚𝑎𝑝𝑟 . (2.21) On the next step and all subsequent steps, the direction of ascent is defined by The steps ( 𝛼 𝛼 𝛼 𝛼𝑛−1 𝑙̃𝑛 𝑛 = 𝑙𝑛 𝑛 + 𝛽𝑛 𝑛 𝑙𝑛−1 . (2.22) ) are defined by a line search method that minimizes the misfit functional: (2.23) The following numerical scheme gives the conjugate gradient method used to solve the minimization of equation (2.16). 𝑟𝑛 = 𝐴(𝑚𝑛 ) − 𝑑 𝛼 𝑙𝑛 𝑛 = 𝑙 𝛼𝑛 (𝑚𝑛 ) = 𝐹𝑚𝑛 𝑇 𝑟𝑛 + 𝛼𝑛 𝑚𝑛 𝛼 𝛼 𝑙̃𝑛 0 = 𝑙𝑛 0 𝛼 𝛼 𝛼 𝛼𝑛−1 𝑙̃𝑛 𝑛 = 𝑙𝑛 𝑛 + 𝛽𝑛 𝑛 𝑙𝑛−1 (2.24) (2.25) (2.26) (2.27) 13 𝛼 𝛽𝑛 𝑛 = 𝛼 𝑘𝑛 𝑛 = 𝛼 𝑙𝑛 𝑛 𝛼 2 𝑛−1 𝑙𝑛−1 𝛼 2. (2.28) 𝛼 𝛼 𝑙̃𝑛 𝑛 , 𝑙𝑛 𝑛 2 𝛼 𝐹𝑚𝑛 𝑙̃𝑛 𝑛 + 𝛼𝐹𝑚𝑛 𝑙̃𝑛 𝑛 2 𝛼 𝛼 mn+1 = mn − 𝑘𝑛 𝑛 ̃l𝑛𝑛 (2.29) (2.30) The regularization method is adaptive. The parameter α is updated in the process of the iterative inversion: 𝛼𝑛 = 𝛼1 𝑞 𝑛−1 ; 𝑛 = 1,2,3, … , ; 0 < 𝑞 < 1. (2.31) The iteration is terminated when the percent error of the norm of the misfit to the norm of the data is less than 1%. 𝑃𝐸 = ‖rn ‖ ∗ 100% < 1% ‖𝑑‖ (2.32) A Monte Carlo approach is used to estimate the uncertainty of the model parameters associated with the measurements. A normal cumulative distribution given a 95% confidence interval is created for each of the input variables. The inverse model then iteratively solves the problem using a randomly generated (within 2σ) input variable. 4. TESTING AND RESULTS 4.1 Equilibration Constant and Calibration The sampler specific equilibration constant was determined experimentally by submerging each sampler in a water bath with an elevated He concentration. The measured He partial pressure inside the sampler was then used in an inverse model to calculate the equilibration constant of each sampler. A Monte Carlo simulation was used to estimate the uncertainty of the calculation. In order to create the He enriched bath at atmospheric pressure, a short length of silicone rubber tubing with a surface area (A) was submerged in a well-mixed volume of water (V) open to the atmosphere. After each sampler was removed from the system, the bath was allowed to equilibrate before the next set was submerged. A continuity equation was created to calculate the necessary equilibration time of the bath. Input of He to the system is defined as G (M/t), or the effusion rate of He, whereas He output can be described by 𝐾𝐴(𝐶 − 𝐶𝑎𝑡𝑚 ) (2.33) Where K is the gas transfer velocity (L/T), C is the He concentration in the tank, and Catm is He concentration in the atmosphere. The time rate of change of mass in the volume can be defined as 𝑉 𝑑𝐶 . 𝑑𝑡 (2.34) Therefore, the time rate of change of mass in the volume is equal to the output of He to 15 the system minus the input of He: 𝑉 𝑑𝐶 = 𝐾𝐴(𝐶 − 𝐶𝑎𝑡𝑚 ) − 𝐺. 𝑑𝑡 (2.35) Therefore, the transient equation where the concentration of He varies with time is given by 𝑑𝑡 𝑑𝐶 = . 𝑉 𝐾𝐴(𝐶 − 𝐶𝑎𝑡𝑚 ) − 𝐺 (2.36) Solving with the initial condition that C =Co at time equal to zero: exp − 𝐾𝐴𝑡 𝐾𝐴𝐶 − 𝐾𝐴𝐶𝑎𝑡𝑚 − 𝐺 = . 𝑉 𝐾𝐴𝐶𝑜 − 𝐾𝐴𝐶𝑎𝑡𝑚 − 𝐺 (2.37) If the initial concentration in the tank (Co) is assumed to be in equilibrium to the atmosphere (Catm), equation (2.37) is then solved for the concentration in the tank as a function of time (C(t)) 𝐶(𝑡) = 𝐶𝑎𝑡𝑚 − 𝐺 𝐾𝐴𝑡 1 − exp − . 𝐾𝐴 𝑉 (2.38) The gas transfer velocity for He was assumed to be 20 to 40 cm/hr [Holmen and Liss, 1984]. The effusion rate of He (G = 0.57 ccstp/hr at 2.2 PSI) for the silicone tubing was estimated based on the effusion rate of oxygen determined by the method of Wilson and Mackay [1995]. These approximations were used to estimate a maximum equilibration time of approximately 20 hours for the enriched He bath. Once the He bath had equilibrated for 24 hours, the samplers were submerged for a time period sufficient to reach 20 to 50% of the equilibrium He concentration. Equilibrium headspace samples from the water bath were taken with the advanced passive diffusion sampler designed by Gardner and Solomon [2009] to determine the equilibrium He partial pressure of the bath. Samplers were then collected and the contents 16 analyzed on a Shimadzu 8AIT Gas Chromatograph (GC) fitted with an AMP-7 amplifier on the thermal conductivity detector (TCD) and a 5A molecular sieve column eight meters long. Data were collected and processed using PeakSimple 3.56. The injection port and detector temperatures as well as the column temperature were maintained at 100oC and 40oC, respectively. Nitrogen (N2) was used as the carrier gas and was set to 0.9 and 3.8 kg/cm2 for carrier lines one and two, respectively. Carrier two goes to an unused 1.8 meter long column. The inlet manifold for the GC is designed to measure the total gas pressure of the sample. The sample is attached to the inlet manifold via an air-tight AN 37⁰ flared fitting (Figure 4.1). In the load position, gas travels via the black path. Carrier gas travels into port 4 and out of port 5 directly to the column. After the sample is attached, the tee valve is opened to a diaphragm vacuum pump and the manifold is evacuated. After closing the pump out valve, the sample clamp is opened allowing the sample gas to fill the manifold volume. Once the pressure of the manifold has been recorded, the manifold is changed to the inject position where gas travels via the gray flow path. The carrier gas then enters port 3 and sweeps the sample gas present in the sample loop (2.0 cc) out of port 5 and into the analytical column. A calibration curve was developed for the GC by injecting three known He standards. The total number of moles (𝑛𝑡𝑜𝑡𝑎𝑙 ) analyzed, those present in the loop, is given by the ideal gas law 𝑛𝑡𝑜𝑡𝑎𝑙 = 𝑃𝑚𝑎𝑛𝑖𝑓𝑜𝑙𝑑 𝑉𝑙𝑜𝑜𝑝 𝑅𝑇 (2.39) where Pmanifold is the pressure measured in the manifold, Vloop is the volume of the loop given by the manufacturer, R is the gas constant, and T is the temperature. The number of 17 Figure 4.1 Gas chromatograph inlet manifold. 18 moles of He analyzed (nHe) can then be determined by injecting a known mole fraction (XHe) of He: 𝑛𝐻𝑒 = 𝑛𝑡𝑜𝑡𝑎𝑙 𝑋𝐻𝑒 = 𝑃𝑚𝑎𝑛𝑖𝑓𝑜𝑙𝑑 𝑉𝑙𝑜𝑜𝑝 𝑋𝐻𝑒 . 𝑅𝑇 (2.40) The precision of the measurements using N2 as a carrier gas is ±5% for He. The minimum detection limit for helium using the manifold GC system is 1.2E-9 moles. An inverse model was used to solve for each sampler's equilibration constant. The equilibration constant (Ec) for a sampler is defined by 𝐸𝑐 = 𝐷𝐴 𝐿 (2.41) where D is the effective diffusion coefficient of He through the membrane, A is the area of the membrane, and L is the membrane thickness. The equilibration constant for each sampler was determined by minimizing the sum of the residuals squared (RSS) between the observed and the calculated partial pressure using equation (2.6) where (PHe(100%)) is the partial pressure measured from the advanced passive diffusion sampler. This sample was assumed to have completely equilibrated with the dissolved He in the tank. The total reservoir volume, both the sampler and sampling tip, was determined gravimetrically. An inverse normal cumulative distribution was calculated for each variable from its specified mean and variance. A Monte Carlo approach was then used to iteratively solve for the equilibration constant and estimate the uncertainty of the calculation. Long equilibration times for the medium and large samplers make it difficult to efficiently collect more than one sample to determine the equilibration constant. Fortunately, the mathematical construction of the exponential curve allows us to invert equation (2.6) for the equilibration constant using only one data point. To prove this 19 assertion, the forward model was used to calculate the expected He partial pressure as a function of time. The predicted He partial pressures at their respective time steps were then used in the RSS model to estimate the equilibration constant. An iterative Monte Carlo process was used to estimate the 95% confidence interval. Figure 4.2 shows the percent error between the forward and RSS models. The circles show the error between the forward model and the average of the 1000 iterations (averaged Ec).The gray window represents the error between forward model and the bounds of the 95% confidence interval. The RSS model most accurately predicts the equilibration constant when the samplers are approximately 75% equilibrated. The error between the forward and RSS model is less than four percent if the inversion is done with samples that are less than 75% equilibrated. Figure 4.3 shows the calculated equilibration curves for a single sampling array. The boxes are the normalized He partial pressure (%) that was measured with the GC. The gray lines are the equilibration curves for the iterations and the black line is the averaged equilibration curve. 4.2 Lab/Model Validation A laboratory test was conducted to test the model's capability of deconvolving an increase in He partial pressure. The sampling array was allowed to fully equilibrate to the atmospheric He concentration (2,000 hrs). The array was then submerged into the elevated He bath for 180 hours. The inverse model was applied to the experimental data collected from each of the samplers in the array. The samples were run on a gas chromatograph as previously described. The partial pressure for each sampler was used in the inverse model to estimate the He pulse magnitude (ΔPHe) and timing (t). The 20 Figure 4.2 RSS error for the equilibration constant (Ec) as a function of equilibration. The circles represent the inverse model error for the average equilibration constant. The gray window represents the 95% confidence interval for the inversion. 21 Figure 4.3 Calibrated equilibration constants and curves for a sampling array. The boxes are the measured values, the black lines are the averaged equilibration curves, and the gray lines are all equilibration curves generated from the Monte Carlo simulation. 22 background concentration was assumed to be the partial pressure of He in the atmosphere. The concentration at collection was determined from two additional passive diffusion samplers [Gardner and Solomon, 2009]. Figure 4.4 shows experimental conditions in black and the model iterations for the magnitude and timing of the He pulse in gray. The averaged He pulse magnitude and timing from the inverse model are within 5% of the imposed experimental conditions. 23 Figure 4.4 One-step inverse model applied to experimental data. Samplers that were equilibrated with the atmosphere were placed into a bath with elevated He concentration for 180 hours. The samplers were removed at 2600 hours. An inverse model was used to estimate the sampler emplacement time (2420 hours) and He concentration in the bath. The gray lines are the iterative results from the inverse model. The black line represents the known experimental conditions. The red lines are the samplers' equilibration curves and the red squares are the measured partial pressure of He in each sampler. 5. MODEL UNCERTAINTY 5.1 Model Experiments A numerical experiment was conducted to evaluate the two-step model's ability to estimate the magnitude of the He change (P2) and the timing of the He pulse (t1) as a function of pulse duration (t2-t1) and the elapsed time between when the pulse occurred and the sampler was collected (t3-t2). The forward model (equation 2.8) was used to calculate the anticipated partial pressure in each sampler under varied pulse intervals and collection times. The forward model assumed the sampling array had completely equilibrated with the in situ background concentration. These data simulations were input to the inverse model to determine the optimal timing conditions where the model predictions would be most accurate. The inverse model predicts the magnitude of the He change and the timing within 30% for all variables (P2, t1, and t2) when the He pulse lasts 20 to 70 hours and the samplers are collected within 100 hours after the pulse (t3-t2<100 hrs) (Figure 5.1). From this numerical experiment, optimal model conditions are defined for further model experiments as follows. The forward model (equation 2.8) assumes that the He partial pressure increased by 100% for 50 hours when it then decreased by 25%. The arrays are then collected 50 hours after the pulse decreased. It is important to know the model's sensitivity to variation in the equilibration constant because measured values of Ec vary, from 0.10 to 0.24. The range of variation 25 Figure 5.1 Model error for varying pulse durations and the time elapsed since the pulse. Gray boxes in the x,y plane represent conditions wherein the inverse model is capable of predicting all variables within 15% error. The black boxes represent the conditions the model can predict within 30% error. 26 is similar for all three sampler sizes and is large enough to make the use of the mean Ec in the inverse model ineffective. The observed variation in Ec is likely due to variability in the area and thickness of the membrane due to stretching during sampler construction. Using optimal model conditions, the calibrated equilibration constant was used in the forward model to estimate the hypothetical partial pressure in each sampler. Calibration of Ec in the laboratory is within ±4% as shown by Figure 5.2. The inverse model was run with 1% variation of Ec to show the model error associated with the laboratory calibrated equilibration constant. On five subsequent model runs, variation of Ec was increased to 10% and then by increments of 10%. When calibrated Ec values with 5% error are used in the inverse model, the percent error for all three model parameters is less than 10% (Figure 5.2). As the variation increases so does the error associated with the predicted model parameters. If variation of Ec exceeds 30%, the model's ability to accurately predict the He partial pressure decreases dramatically. The timing of the pulse also becomes increasingly difficult to predict; the percent error of t1 and t2 increases to ~55 and 30%, respectively. If an average of all the calibrated Ec values and the standard deviation of these values is used in the inverse model for all of the samplers, the model's predictions of the timing and magnitude of the He pulse are highly uncertain. The error for the timing (t1 and t2) and magnitude (P2) increases to 1.5 and 2,000 orders of magnitude, respectively. From these simulations, it is clear that the equilibration constant of each sampler must be calibrated in order for the model to provide accurate predictions. 27 Figure 5.2 Model error as a function of equilibration constant variation. The circles show the average model error for a particular variable. The windows represent the 95% confidence interval for the model error associated with the inversion. The partial pressure percent error is in a log scale. 28 5.2 Temperature Dependence of Helium Diffusion The effective diffusion coefficient (D) of a species is a measure of mobility and is a function of temperature (T) and the activation energy for diffusion (Q): 𝐷 = 𝐷𝑜 exp − 𝑄 𝑅𝑇 (2.42) where Do is the temperature independent pre-exponential and R is the gas constant [Arrhenius, 1889]. The activation energy and pre-exponential can be determined by plotting the effective diffusion coefficient against the reciprocal of T (in Kelvin). Temperature in and among volcanic springs and vents varies greatly. In order to determine the potential magnitude of the temperature effect on the inverse model, the temperature dependence of He diffusion in the semipermeable silicone membrane was determined experimentally. Copper tubes fitted with semipermeable silicone membranes were filled with pure He at atmospheric pressure. The tubes were cold welded and allowed to equilibrate with the atmosphere via the silicone membrane at room temperature (~21⁰C), 50⁰C, and 100⁰C for a specified period of time (Figure 5.3). Two membrane-fitted, copper tubes were constructed to allow sample collection at two different time steps for each temperature. To collect the gas samples, the copper tubes were pinch clamped. These gas samples were analyzed by the GC method described in section 2.3. Equation (2.3) was altered so that the partial pressure in the tube was equal to pure He that then equilibrated to the atmospheric partial pressure of He. 29 Figure 5.3 Copper (Cu) tubes fitted with semipermeable silicone membranes. Cu tubes were used to determine the temperature dependence of He diffusion. The tubes were filled with pure He, cold welded (1), and then allowed to equilibrate with the atmosphere for a specified period of time. The tubes were cold welded (2) to collect gas samples for analysis. 30 𝑃𝐻𝑒(𝑡𝑢𝑏𝑒) = 𝑃𝐻𝑒(𝑃𝑢𝑟𝑒) exp − 𝐷𝐴 𝑡 𝑉𝐿 𝐷𝐴 + 𝑃𝐻𝑒(𝑎𝑡𝑚) 1 − exp − 𝑡. 𝑉𝐿 (2.43) The area and thickness of the membranes were measured. The internal volume of each equilibration tube was calculated by measuring the internal diameter and length of each component (copper tubes and membrane). The effective diffusion coefficient for the membrane at each temperature was determined by minimizing the sum of the residuals squared (RSS) between the observed and the calculated partial pressure of He using equation (2.43). The RSS model was run four times for each temperature. On each run, a different combination of samples was used. For example, for a given temperature two samples were collected at time step one, 1A and 1B, and two were collected at time step 2, 2A and 2B. On the first run of the RSS model, samples 1A and 2A would be used to invert for the diffusion constant. On the second run, samples 1B and 2B would be used. The same Monte Carlo approach was used to estimate uncertainty. An Arrhenius equation was created to describe the temperature effect on diffusion. A trust-region algorithm was used to fit a nonlinear regression to the data and generate the 95% confidence interval (Figure 5.4). The best-fit Arrhenius equation for these data is given by: 𝐷 = 44.39 𝑐𝑚2 1000 ∙ exp −2.001 . ℎ𝑟 𝑇 (2.44) Estimations from the Arrhenius equation will be used in the inverse model to determine its ability to accurately predict the model parameters at various temperatures. In order to do that, it is necessary to convert the He gas effective diffusion coefficient to 31 Figure 5.4 Temperature dependence of He diffusion. Circles represent the effective diffusion coefficient for He gas at various temperatures from the RSS model. The solid line is the best fit Arrhenius equation for the data. The dotted lines show the 95% confidence interval of the regression. 32 an aqueous equilibration constant that can be used in the model. Recall that the aqueous equilibration constant is given by 𝐸𝑐 = 𝐷𝐴 . 𝐿 (2.1) The area and thickness of the membrane and conversion factor from diffusion in air to diffusion in water can be represented as a constant (C). 𝐸𝑐 = 𝐶 ∙ 𝐷 (2.2) Equation (2.44) was used to calculate the room temperature (21⁰C) effective diffusion coefficient in air. The equilibration constant for the sampler is known from lab calibration. The conversion factor was estimated for each sampler via RSS model. The lower, upper, and best fit equilibration constant for each sampler was then calculated at various temperatures for use in the inverse model (Figure 5.5). To evaluate the temperature effect, the model was tested under various temperatures. In the forward model, the lower, upper, and best fit Ec values at each temperature (0, 20, 50, 75, 100⁰C) were used to calculate the expected He partial pressure in the samplers. In the first set of inverse simulations, Panel A in Figure 5.6, the calibrated, room temperature Ec values were used in the inverse model. This experiment shows the model results associated with ignoring the temperature of the sampling site. In the second set of simulations, Panel B in Figure 5.6, the Ec values were temperature corrected. Error bars in this figure represent the error associated with the Monte Carlo iterations of that simulation. If the temperature of the sampling site is ignored, the model's capability of reproducing the parameters (P1, t1, and t2) greatly diminishes. At all temperatures, use of the upper limit Ec does not produce realistic predictions. The best fit limit Ec values 33 Figure 5.5 Equilibration constant (Ec) as a function of temperature. The boxes represent the best-fit Ec calculated from equation (2.44) and the bars show the lower and upper 95% confidence interval for that temperature. The circles represent calibrated Ec values at 20⁰C. 34 Figure 5.6 Effect of temperature on two-step inverse model results. The temperature corrected Ec value was used in the forward model. Panel A represents the inverse model results when the room temperature Ec value was used. Panel B represents the inverse results when Ec has been corrected for temperature. 35 predict the model parameters within 45% for 0 and 20⁰C. However, at higher temperatures the error of the predicted partial pressure within a given simulation increases indicating that the model becomes unstable under these conditions. If the equilibration constant is temperature corrected, use of the best fit Ec predicts the timing of the pulse within 20%. The partial pressure is more difficult to predict at higher temperatures above 20⁰C. Error for these higher temperatures ranges from 51% at 50⁰C to 130% at 100⁰C. It is probable that this error increase is due to the accelerated equilibration time. As the samplers take less time to equilibrate, the ‘optimal conditions' used for the simulations are too long to accurately record the pulse magnitude. If the temperature of the sampling site is ignored, the model's capability of reproducing the parameters (P1, t1, and t2) greatly diminishes. At all temperatures, use of the upper limit Ec does not produce realistic predictions. The best fit limit Ec values predict the model parameters within 45% for 0 and 20⁰C. However, at higher temperatures the error of the predicted partial pressure within a given simulation increases indicating that the model becomes unstable under these conditions. If the equilibration constant is temperature corrected, use of the best fit Ec predicts the timing of the pulse within 20%. The partial pressure is more difficult to predict at higher temperatures above 20⁰C. Error for these higher temperatures ranges from 51% at 50⁰C to 130% at 100⁰C. It is probable that this error increase is due to the accelerated equilibration time. As the samplers take less time to equilibrate, the ‘optimal conditions' used for the simulations are too long to accurately record the pulse magnitude. 36 5.3 Environmental Effects on Diffusion A potential problem for long term deployment of the samplers is the degradation of the membrane over time. Membrane degradation could affect the effective diffusion coefficient of He and thus the model's ability to deconvolve a concentration change. Factors that might alter the permeability of silicone membranes include colloidal, biological, organic, and mineral precipitation. Sampling conditions such as pH and temperature may also affect He permeability. In order to determine some of the effects of prolonged deployment on the equilibration constant, calibration tests were conducted on sampling arrays deployed on the flanks of Mammoth Mountain. Once the samplers were returned, they were calibrated to determine the equilibration constant via the method described in section 2.3. The membrane of each sampler was then replaced with a new membrane of the same area and thickness. The samplers were then recalibrated to determine the equilibration constant of the new membrane. Three arrays were deployed at routine United States Geological Survey (USGS) water and gas sampling sites as a part of the Volcano Hazards Program: Mammoth Mountain Fumarole, Artesian Soda Springs, and Horseshoe Lake. The Mammoth Mountain Fumarole is a steam vent at 3030 meters (m) on the north side of Mammoth Mountain. The USGS has collected temperature, steam flow rate, and chemical and isotopic composition of the fumarole gas since 1989 when distinct steam vents appeared. A plastic pipe was installed in August 1991 in the highest temperature vent to facilitate sample collection [Sorey et al., 1998]. The data collected through 1992 were reported by Sorey et al. [1993]. 37 The sampling array (array #3) was buried two feet to the northeast of the plastic pipe, 30 centimeters (cm) deep, in an acrylonitrile-butadiene-styrene (ABS) casing open at the bottom. At deployment, the steam vent was 81.6⁰C and the soil at the bottom of the array case was 44.5⁰C. The array was buried in the soil next to the steam vent on August 19, 2013 and retrieved on October 22, 2013. Artesian Soda Springs is a cold water spring (6-8⁰C) part of a system of springs along Lower Boundary Creek on the southwestern flank of the mountain. The high-flow spring is characterized by low pH (~5.4) and a strong CO2 signal (~20 mmol/l) [Evans et al., 2002]. The array (#4) was submerged approximately two feet below the water surface. It was deployed on August 20, 2013 and collected at October 22, 2013. Horseshoe Lake is the 30-hectare region of killed trees from which approximately 200 tons of CO2 per day was emitted at the time of deployment (J. Lewicki, unpublished data, 2014). Soils in this area show normal temperatures (14.3⁰C at deployment) and no visible signs of acid alteration. The array (#1) was buried 36 cm deep in an ABS casing open at the bottom. It was deployed on August 19, 2013 and collected on October 21, 2013. Figure 5.7 shows the equilibration constants for the samplers deployed at the various sampling sites contrasted against those same samplers with new, replaced membranes. At the soil gas deployment sites, there is no consistent pattern of membrane alteration. At Artesian Soda Springs, the equilibration constant for the used samplers is ~12% larger than the newly replaced membranes indicating the used sampler membranes may equilibrate more quickly than the new membranes. 38 Figure 5.7 Environmental effects on the equilibration constant. Red circles represent the equilibration constant for membranes that were deployed for two months while blue circles represent the equilibration constant for the same sampler but with a newly replaced membrane. 39 Microbial growth has been shown to reduce the performance of silicone samplers in some environments. One study utilizing silicone membranes for equilibrium headspace sampling in anoxic paddy soil observed microbial growth that decreased the permeability of H2 in the membrane within a matter of days [Krämer and Conrad, 1993]. On the other hand, Jacinthe [2001] did not detect any microbial growth on their semipermeable silicone membrane samplers after one year of deployment in water-saturated soils. Little to no microbial growth was observed on samplers deployed at any of the sites for nearly two months. In addition, the equilibration time for samplers deployed at Artesian Soda Springs decreased with respect to that of the new membranes. For the arrays deployed in soil, the ABS casing used to enclose the samplers may have prevented membrane alteration. Membrane fouling may be more likely if the samplers are deployed directly into the soil. The equilibration constants for the deployed samplers are well within the range of those not exposed to harsh environmental conditions and there seems to be little to no effect of these conditions on the equilibration constant. The differences between the equilibration constants for each sampler may be attributed to the total area and thickness of the membranes rather than environmental fouling. Although this test does not evaluate all the effects of the environment, it does give insight as to how to calibrate the sampling arrays. As it may be desirable to deploy some arrays for a much longer period of time, it is possible that longer times of exposure may have an effect on the equilibration constant. It is recommended that the sampling arrays be calibrated after they have been collected. This calibrated Ec should then be representative of the previous hundreds of hours in question. 6. FIELD RESULTS The arrays were deployed at various USGS Volcano Hazards Program monitoring sites (Manzanita Creek at Mount Lassen, Carbonate Springs at Mount St. Helens, and Medicine Lake Hot Spot at Medicine Lake Volcano) in 2011 and again in 2013. These sites have had a history of R/Ra values that are different from air and have been known to fluctuate with magmatic activity [Hilton, 1996; Evans et al., 2004; Wicks et al., 2002]. The samplers deployed in 2011 were not calibrated to determine the equilibration constant as its importance for inverse modeling was not known at the time. The results from these deployments do, however, demonstrate the array's capability to record the He concentration and isotope ratio at the site (Figure 6.1). For the following three sites, a copper tube sample was collected when the arrays were deployed. Manzanita Creek tributary spring lies approximately two miles west of a weak fumarolic area on the northern flank of Mount Lassen, the southernmost active volcano in the Cascade Range. This cold spring (~4⁰C) is part of an extensive monitoring network and from 2009 to 2012 was sampled hourly for bicarbonate flux [Ingebritsen et al., 2014]. The array was deployed for 79 days from July 11 to September 28, 2011. The R/Ra ratio for the smallest sampler seems to track that of the deployment very well (within 2%). Given the possible range of equilibration constants, the differences among the medium, large, and Cu tube ratios may be due to incomplete equilibration from atmosphere to the site's R/Ra ratio rather than tracking a change in the R/Ra history at the 41 Figure 6.1 R/Ra ratios for 2011 deployment. R/Ra ratios for sampling arrays versus those from Cu tube samples taken at deployment. The light gray crosshairs represent the ratio from the small sampler, medium gray crosshairs are those from the medium samplers, and dark gray crosshairs are those from the large samplers. Arrays were collected between 49 and 79 days after deployment. 42 site. In the same study by Ingebritsen et al., Carbonate Springs on Mount St. Helens was sampled hourly from 2009 to 2012 for chloride flux. The array was deployed during this sampling period for 84 days from July 11 to October 3, 2011. The results from this spring indicate no significant increase in the R/Ra ratio at this site. The entire array closely tracks (within 3.3%) the ratio determined at deployment indicating there was no significant ratio fluctuations over the deployment period. Medicine Lake Hot Spot is located in the caldera of Medicine Lake Volcano, California. It is continuously monitored by the USGS for temperature and is sampled intermittently. This array was deployed for 49 days from August 11 to September 29, 2011. Similar to the results from the Manzanita Creek tributary spring, the R/Ra ratio for the smallest sampler is within 5% of the copper tube ratio. The systematic shift of the medium and large sampler away from the 1:1 line is likely to due to incomplete equilibration from atmosphere to the site's R/Ra ratio rather than tracking a change in the R/Ra history at the site. The array deployed at Horseshoe Lake (HSL) in 2013 had significant R/Ra variability among the three samplers and the data were applied in the one-step inverse model. Unfortunately, the copper tube sample taken at collection leaked to the atmosphere. Therefore, for modeling purposes, the ratio at collection was assumed to be equal to that of the smallest sampler. This sampler is 90% equilibrated within 48 hours. The aqueous equilibration constant for the array was determined, at room temperature, by the method described; however, the array was buried in 14.3⁰C soil. The calibrated equilibration constant was corrected for nonaqueous conditions at 14⁰C using 43 equations 2.46 and 2.44. The 4He concentration did not change through time (<1% difference among the samplers). Therefore, the inverse model was run for the 3He concentration using the best fit, upper, and lower Ec values. The results using the lower and best fit Ec value were unrealistic given the data and are not reported. In order to demonstrate the range of probable histories, multiple simulations were run using Ec values ranging from the best fit to the upper limit Ec. The R/Ra ratio was then calculated using the results from the 3He data inverse and the average 4He concentration of the array. Figure 6.2 shows the model's predicitions of the previous ratio and the time at which the ratio increased. The model simulations predict an R/Ra value of 4.76±0.12 before the increase occurred 1.9±0.6 days before the array was collected. Seismic events (M>0.5) for the area from mid-June to December 2013 were accessed through the Northern California Earthquake Data Center (NCEDC) and are plotted on Figure 6.3. The USGS deems the documented seismic activity during this time period as ‘typical behavior' for the Long Valley caldera. Generally, the magnitude of this activity ranges from 0.1 (not shown) to ~1.5 M and occurs from 2 to 5 km and 9 to 11 km deep. All of the seismic activity plotted is within 5 km of HSL. A conceptual model of a sealed, low-temperature (~150⁰C), high-pressure gas reservoir, as shallow as 2 km deep, has been proposed to describe the relationship between seismic, chemical, and isotopic observations at Mammoth Mountain [Sorey et al., 1998]. They propose that, in 1989, a break in this high permeability seal was enlarged via intrusion or associated seismicity and is responsible for the change in the 3He/4He ratio to above 4 Ra that occurred from 1989 to mid-1990. Since then, the He isotope ratio 44 Figure 6.2 One-step inverse modeling of Horseshoe Lake data. The black box represents the R/Ra ratio at collection and was assumed to be equal to that of the small sampler at collection. The red boxes indicate the R/Ra ratio in each sampler at collection. The gray boxes represent the previous ratio and the timing of the change using a range of Ec (from the best fit to upper limit) in the one-step inverse model. 45 Figure 6.3 Seismic activity (M>0.5) in the Mammoth Mountain area. The circle overlain by a pentagon is a seismic event less than 1 km deep. The double circles represent seismic events that occurred from 2 to 4 km depth. The squared circles are those from 9 to 11 km deep. The color of each event represents the distance from HSL. 46 has fluctuated in response to shallow earthquake swarms and long period earthquakes at Mammoth Mountain [Sorey et al., 1998]. Via this mechanism, unlike the study by Padron [2013] at El Hierro, an increase in the ratio seems to occur after seismic activity. Sorey et al. [1993] estimated minimum volatile transport rates for He at Mammoth Mountain Fumarole after intrusive activity in 1989. Given their transport rate of 30 to 60 m/d and assuming the flow distance is equal to the depth of the seismic event, it could take 33 to 66 days for the volatiles to travel 2 km or 83 to 166 days for 5 km. If we combine these assumptions with the model results, seismic events 2 to 5 km deep occurring between May 6 and September 16 could be responsible for an increase in the He ratio. It is also possible that the volatile transport rate could be slower at HSL. The lack of a fumarolic vent may suggest a decrease in vertical permeability and subsequent longer travel times. Although it is not possible to discern precisely which seismic events are responsible for the increase in the R/Ra ratio from these data, it is definitive that the sampling array recorded an in situ ratio change. Continued use of the sampling array at this site could provide further insight to the relationship between seismic events, magmatic activity, and He isotope ratio fluctuations. 7. SUMMARY AND CONCLUSIONS Passive diffusion headspace samplers have been used successfully to determine the He isotope ratio for crustal degassing [Sheldon, 2002], groundwater recharge [Gascoyne and Sheppard, 1993], and hydrothermal circulation [Gardner et al., 2010] studies. Passive diffusion samples do not require laboratory extraction and are therefore more efficient to analyze. Samples collected via passive diffusion have been shown to be within 1% error of those analyzed from copper tubes [Gardner and Solomon, 2009]. In addition to other passive sampling methods, our array is useful for determining changes in the ratio through time and, given the ratio is static during deployment, for demonstrating a lack of significant fluctuations between sampling periods. The array is capable of passively recording He concentration and isotope ratio in spring water and soil gas. The unique feature of this array is that it can be deployed for long periods of time and is intended to be collected after volcanic unrest is detected by other means. An inverse model was built to reconstruct the timing and magnitude of in situ He pulses that occur before the samplers are collected. The two-step inverse model requires that additional water or gas samples be collected during deployment and collection. Once collected, each of the samplers in the array is calibrated to determine its unique equilibration time. Each sampler's unique equilibration time enables the array to record the He pulse at three different rates. When combined in the inverse model, the background, collection, and partial pressures of He in 48 each of the samplers is used to reconstruct the timing and magnitude of the He pulse. The associated two-step inverse model is capable of accurately (<30% error) reconstructing pulses that last between 20 and 70 hours when the samplers are collected within 100 hours of that pulse. Volumetrically larger sampling arrays or the addition of a larger sampler to the array may be capable of reconstructing pulses after a much longer time period from the end of the pulse to collection. Initial field results indicate the array and inverse model are capable of capturing and reconstructing in situ R/Ra ratio changes that occur before the array is collected. The sampling array provides an accurate and affordable alternative for capturing He anomalies that may be associated with volcanic unrest. APPENDIX A ONE-STEP INVERSE MODEL MATLAB CODE %data input var = xlsread('dataimport'); Ec = var(1,1:3); %[s m l] (cm^3/hr) Vs = var(3,1:3); %[s m l] (cm^3) d = var(5,1:3)'; %[s m l] (pressure torr) P2 = var(9,1); %partial pressure of He at collection (torr) Patm = 5.2e-6*646; %partial pressure of He in atmosphere (torr) t2 = var(7,1); %time elapsed since deployment (hrs) %generate values from the uniform distribution on the interval [0.025,0.975] a = 0.025; b = 0.975; r = a + (b-a).*rand(1000,1); %create inverse of the normal cumulative distribution function for each %varible Ec_r = ones(1000,1); for i=1:length(r); for j=1:3; Ec_r(i,j) = norminv(r(i),var(1,j),var(2,j)); end end V_r = ones(1000,3); for i=1:length(r); for j=1:3; V_r(i,j) = norminv(r(i),var(3,j),var(4,j)); end end d_r = ones(1000,3); for i=1:length(r); for j=1:3; d_r(i,j) = norminv(r(i),var(5,j),var(6,j)); end end d_r=d_r'; t2_r = ones(1000,1); for i=1:length(r); t2_r(i,1) = norminv(r(i),var(7,1),var(8,1)); 50 end P2_r = ones(1000,1); for i=1:length(r); P2_r(i,1) = norminv(r(i),var(9,1),var(10,1)); end %%Then, the starting model (m_0) is defined and the initial descent %%direction is calculated m_0 = [P2 t2]'; %create index and matrices to store variables from i-loop idx=1; idy=1; idz=1; m_iter=ones(1000,2); d_inv_iter=ones(3,1000); A_0_iter=ones(3,1000); %%i-loop to iteratively solve for m using cdf distribution of %%variables for i=1:1000; %calculate the frechet derivative matrix for m_0 frech_0 = frechetpressure(m_0(1,1),m_0(2,1),t2_r(i,1),P2_r(i,1),... Ec_r(i,1:3),V_r(i,1:3)); A_0 = Ampressure(m_0(1,1),m_0(2,1),t2_r(i,1),P2_r(i,1),Patm,... Ec_r(i,1:3),V_r(i,1:3)); %calculate residual vector for m r_0 = A_0-d_r(1:3,i); %find optimal model perturbation l_0 = frech_0'*r_0; ltilda_0= l_0; ktilda_0= norm(ltilda_0)^2/(norm(frech_0*ltilda_0))^2; del_m_0= -ktilda_0*ltilda_0; %calculate m1 to start for loop m = m_0+del_m_0; %set initial l vectors l_n_1 = l_0; ltilda_n_1 = ltilda_0; %%Start conjugate gradient descent method until the percent error of d is %%less than 1 for j=1:10^6; %calculate frechet derivative for m frech = frechetpressure(m(1,1),m(2,1),t2_r(i,1),P2_r(i,1),... Ec_r(i,1:3),V_r(i,1:3)); %calculate residual vector for m r = Ampressure(m(1,1),m(2,1),t2_r(i,1),P2_r(i,1),Patm,... Ec_r(i,1:3),V_r(i,1:3))-d_r(1:3,i); %calculate alpha alpha = norm(r)^2/norm(m)^2; %calculate steepest ascent vector l_n = frech'*r+alpha*m; %calculate beta beta = (norm(l_n))^2/(norm(l_n_1))^2; %calculate ltilda_n ltilda_n = l_n+beta*ltilda_n_1; 51 %find optimal model perturbation k = ltilda_n'*l_n/((norm(frech*ltilda_n))^2+alpha*(norm(ltilda_n))^2); del_m = -k*ltilda_n; %update model to m1 m = m+del_m; %set tolerance PE = norm(r)/norm(d_r(1:3,i))*100; if PE<1 break; end %reset variables for the next loop l_n_1=l_n; ltilda_n_1=ltilda_n; alpha = alpha*0.8; end %store values from iteration m_iter(idx,1:2)=m; idx=idx+1; d_inv_iter(1:3,idy)=d_r(1:3,i); idy=idy+1; A_0_iter(1:3,idz)=A_0; idz=idz+1; end %save m_iter file m_exp=m_iter; save('m_exp','m_exp'); APPENDIX B ONE-STEP FORWARD MODEL MATLAB CODE function A = Ampressure(P1,t1,t2,P2,Patm,Ec,V) A = [Patm*exp(-Ec(1,1)*t2/V(1,1))+P1*(1-exp(-Ec(1,1)*t2/V(1,1)))+(P2-P1)*(1-exp(-Ec(1,1)*(t2t1)/V(1,1))) Patm*exp(-Ec(1,2)*t2/V(1,2))+P1*(1-exp(-Ec(1,2)*t2/V(1,2)))+(P2-P1)*(1-exp(-Ec(1,2)*(t2t1)/V(1,2))) Patm*exp(-Ec(1,3)*t2/V(1,3))+P1*(1-exp(-Ec(1,3)*t2/V(1,3)))+(P2-P1)*(1-exp(-Ec(1,3)*(t2t1)/V(1,3)))]; end APPENDIX C ONE-STEP FRECHET DERIVATIVE MATLAB CODE function F = frechetpressure(P1,t1,t2,P2,Ec,V) F = [(1-exp(-Ec(1,1)*(t2)/V(1,1)))-(1-exp(-Ec(1,1)*(t2-t1)/V(1,1))) (-Ec(1,1)*(P2-P1)*(exp(-Ec(1,1)*(t2t1)/V(1,1))))/V(1,1) (1-exp(-Ec(1,2)*(t2)/V(1,2)))-(1-exp(-Ec(1,2)*(t2-t1)/V(1,2))) (-Ec(1,2)*(P2-P1)*(exp(-Ec(1,2)*(t2t1)/V(1,2))))/V(1,2) (1-exp(-Ec(1,3)*(t2)/V(1,3)))-(1-exp(-Ec(1,3)*(t2-t1)/V(1,3))) (-Ec(1,3)*(P2-P1)*(exp(-Ec(1,3)*(t2t1)/V(1,3))))/V(1,3)]; end APPENDIX D TWO-STEP INVERSE MODEL MATLAB CODE %data input var = xlsread('dataimport'); Ec = var(1,1:3); %[s m l] (cm^3/hr) Vs = var(3,1:3); %[s m l] (cm^3) d = var(5,1:3)'; %[s m l] (pressure torr) P1 = var(9,1); %partial pressure of He at deployment (torr) P3 = var(11,1); %partial pressure of He at collection (torr) Patm = 5.2e-6*646; %partial pressure of He in atmosphere (torr) tcoll = var(7,1); %time elapsed since deployment (hr) %generate values from the uniform distribution on the interval [0.025,0.975] a = 0.025; b = 0.975; r = a + (b-a).*rand(1000,1); %create inverse of the normal cumulative distribution function for each %varible Ec_r = ones(1000,1); for i=1:length(r); for j=1:3; Ec_r(i,j) = norminv(r(i),var(1,j),var(2,j)); end end V_r = ones(1000,3); for i=1:length(r); for j=1:3; V_r(i,j) = norminv(r(i),var(3,j),var(4,j)); end end d_r = ones(1000,3); for i=1:length(r); for j=1:3; d_r(i,j) = norminv(r(i),var(5,j),var(6,j)); end end d_r=d_r'; tcoll_r = ones(1000,1); for i=1:length(r); 55 tcoll_r(i,1) = norminv(r(i),var(7,1),var(8,1)); end P1_r = ones(1000,1); for i=1:length(r); P1_r(i,1) = norminv(r(i),var(9,1),var(10,1)); end P3_r = ones(1000,1); for i=1:length(r); P3_r(i,1) = norminv(r(i),var(11,1),var(12,1)); end %%Then, the starting model (m_0) is defined and the initial descent %%direction is calculated m_0 = [P3 tcoll tcoll]'; %create index and matrices to store variables from i-loop idx=1; idy=1; idz=1; m_iter=ones(1000,3); d_inv_iter=ones(3,1000); A_0_iter=ones(3,1000); %%i-loop to iteratively solve for m using cdf distribution of %%variables for i=1:1000; %calculate the frechet derivative matrix for m_0 frech_0 = frechetpressure(m_0(1,1),m_0(2,1),m_0(3,1),tcoll_r(i,1),P1_r(i,1),P3_r(i,1),... Ec_r(i,1:3),V_r(i,1:3)); A_0 = Ampressure(m_0(1,1),m_0(2,1),m_0(3,1),tcoll_r(i,1),Patm,P1_r(i,1),P3_r(i,1),... Ec_r(i,1:3),V_r(i,1:3)); %calculate residual vector for m r_0 = A_0-d_r(1:3,i); %find optimal model perturbation l_0 = frech_0'*r_0; ltilda_0= l_0; ktilda_0= norm(ltilda_0)^2/(norm(frech_0*ltilda_0))^2; del_m_0= -ktilda_0*ltilda_0; %calculate m1 to start for loop m = m_0+del_m_0; %set initial l vectors l_n_1 = l_0; ltilda_n_1 = ltilda_0; %%Start conjugate gradient descent method until the percent error of d is %%less than 1 for j=1:10^6; %calculate frechet derivative for m frech = frechetpressure(m(1,1),m(2,1),m(3,1),tcoll_r(i,1),P1_r(i,1),P3_r(i,1),... Ec_r(i,1:3),V_r(i,1:3)); %calculate residual vector for m r = Ampressure(m(1,1),m(2,1),m(3,1),tcoll_r(i,1),Patm,P1_r(i,1),P3_r(i,1),... Ec_r(i,1:3),V_r(i,1:3))-d_r(1:3,i); 56 %calculate alpha alpha = norm(r)^2/norm(m)^2; %calculate steepest ascent vector l_n = frech'*r+alpha*m; %calculate beta beta = (norm(l_n))^2/(norm(l_n_1))^2; %calculate ltilda_n ltilda_n = l_n+beta*ltilda_n_1; %find optimal model perturbation k = ltilda_n'*l_n/((norm(frech*ltilda_n))^2+alpha*(norm(ltilda_n))^2); del_m = -k*ltilda_n; %update model to m1 m = m+del_m; %set tolerance PE = norm(r)/norm(d_r(1:3,i))*100; if PE<1 break; end %reset variables for the next loop l_n_1=l_n; ltilda_n_1=ltilda_n; alpha = alpha*0.8; end %store values from iteration m_iter(idx,1:3)=m; idx=idx+1; d_inv_iter(1:3,idy)=d_r(1:3,i); idy=idy+1; A_0_iter(1:3,idz)=A_0; idz=idz+1; end %save m_iter file m_exp=m_iter; save('m_exp','m_mexp'); APPENDIX E TWO-STEP FORWARD MODEL MATLAB CODE function A = Ampressure(P2,t1,t2,t3,Patm,P1,P3,Ec,V) A = [Patm*exp(-Ec(1,1)*t3/V(1,1))+P1*(1-exp(-Ec(1,1)*t3/V(1,1)))+(P2-P1)*(1-exp(-Ec(1,1)*(t3t1)/V(1,1)))+(P3-P2)*(1-exp(-Ec(1,1)*(t3-t2)/V(1,1))) Patm*exp(-Ec(1,2)*t3/V(1,2))+P1*(1-exp(-Ec(1,2)*t3/V(1,2)))+(P2-P1)*(1-exp(-Ec(1,2)*(t3t1)/V(1,2)))+(P3-P2)*(1-exp(-Ec(1,2)*(t3-t2)/V(1,2))) Patm*exp(-Ec(1,3)*t3/V(1,3))+P1*(1-exp(-Ec(1,3)*t3/V(1,3)))+(P2-P1)*(1-exp(-Ec(1,3)*(t3t1)/V(1,3)))+(P3-P2)*(1-exp(-Ec(1,3)*(t3-t2)/V(1,3)))]; end APPENDIX F TWO-STEP FRECHET DERIVATIVE MATLAB CODE function F = frechetpressure(P2,t1,t2,t3,P1,P3,Ec,V) F = [(1-exp(-Ec(1,1)*(t3-t1)/V(1,1)))-(1-exp(-Ec(1,1)*(t3-t2)/V(1,1))) (-Ec(1,1)*(P2-P1)*(exp(Ec(1,1)*(t3-t1)/V(1,1))))/V(1,1) (-Ec(1,1)*(P3-P2)*(exp(-Ec(1,1)*(t3-t2)/V(1,1))))/V(1,1) (1-exp(-Ec(1,2)*(t3-t1)/V(1,2)))-(1-exp(-Ec(1,2)*(t3-t2)/V(1,2))) (-Ec(1,2)*(P2-P1)*(exp(-Ec(1,2)*(t3t1)/V(1,2))))/V(1,2) (-Ec(1,2)*(P3-P2)*(exp(-Ec(1,2)*(t3-t2)/V(1,2))))/V(1,2) (1-exp(-Ec(1,3)*(t3-t1)/V(1,3)))-(1-exp(-Ec(1,3)*(t3-t2)/V(1,3))) (-Ec(1,3)*(P2-P1)*(exp(-Ec(1,3)*(t3t1)/V(1,3))))/V(1,3) (-Ec(1,3)*(P3-P2)*(exp(-Ec(1,3)*(t3-t2)/V(1,3))))/V(1,3)]; end REFERENCES Andrews, J. (1985), The isotopic composition of radiogenic helium and its use to study groundwater movement in confined aquifers, Chemical Geology, 49(1), 339-351. Arrhenius, S. (1889), Über die Reaktionsgeschwindigkeit bei der Inversion von Rohrzucker durch Säuren, Zeitschrift für physikalische Chemie, 4, 226-248. Brown, R., J. Charlton, and K. Saunders (1981), The development of an improved diffusive sampler, The American Industrial Hygiene Association Journal, 42(12), 865-869. Carapezza, M. L., S. Inguaggiato, L. Brusca, and M. Longo (2004), Geochemical precursors of the activity of an open-conduit volcano: The Stromboli 2002-2003 eruptive events, Geophysical Research Letters, 31(7), L07620. Clarke, W. B., M. Beg, and H. Craig (1969), Excess 3He in the sea: Evidence for terrestrial primodal helium, Earth and Planetary Science Letters, 6(3), 213-220. Craig, H., J. Lupton, and Y. Horibe (1978), A mantle helium component in circumPacific volcanic gases: Hakone, the Marianas, and Mt. Lassen, Terrestrial Rare Gases (1978), 3-16. Craig, H., J. Lupton, J. Welhan, and R. Poreda (1978), Helium isotope ratios in Yellowstone and Lassen Park volcanic gases, Geophysical Research Letters, 5(11), 897-900. Evans, W. C., M. C. van Soest, R. H. Mariner, S. Hurwitz, S. E. Ingebritsen, C. W. Wicks, and M. E. Schmidt (2004), Magmatic intrusion west of Three Sisters, central Oregon, USA: The perspective from spring geochemistry, Geology, 32(1), 69-72. Evans, W. C., M. L. Sorey, A. C. Cook, B. M. Kennedy, D. L. Shuster, E. M. Colvard, L. D. White, and M. A. Huebner (2002), Tracing and quantifying magmatic carbon discharge in cold groundwaters: Lessons learned from Mammoth Mountain, USA, Journal of Volcanology and Geothermal Research, 114(3-4), 291-312. Gardner, P., and D. K. Solomon (2009), An advanced passive diffusion sampler for the determination of dissolved gas concentrations, Water Resources Research, 45(6), W06423. 60 Gardner, W. P., D. D. Susong, D. K. Solomon, and H. P. Heasler (2010), Using noble gases measured in spring discharge to trace hydrothermal processes in the Norris Geyser Basin, Yellowstone National Park, U.S.A, Journal of Volcanology and Geothermal Research, 198(3-4), 394-404. Gascoyne, M., and M. I. Sheppard (1993), Evidence of terrestrial discharge of deep groundwater on the Canadian Shield from helium in soil gases, Environmental Science & Technology, 27(12), 2420-2426. Hilton, D. R. (1996), The helium and carbon isotope systematics of a continental geothermal system: Results from monitoring studies at Long Valley caldera (California, U.S.A.), Chemical Geology, 127(4), 269-295. Holmen, K. I. M., and P. Liss (1984), Models for air-water gas transfer: An experimental investigation, Tellus B, 36B(2), 92-100. Ingebritsen, S. E., et al. (2014), Hydrothermal monitoring in a quiescent volcanic arc: Cascade Range, northwestern United States, Geofluids, doi: 10.1111/gfl.12079, in press. Jacinthe, P. A., and P. M. Groffman (2001), Silicone rubber sampler to measure dissolved gases in saturated soils and waters, Soil Biology and Biochemistry, 33(7-8), 907912. Kot, A., B. Zabiegała, and J. Namieśnik (2000), Passive sampling for long-term monitoring of organic pollutants in water, TrAC Trends in Analytical Chemistry, 19(7), 446-459. Krämer, H., and R. Conrad (1993), Measurement of dissolved H< sub> 2</sub> concentrations in methanogenic environments with a gas diffusion probe, FEMS microbiology ecology, 12(3), 149-158. Lupton, J., and H. Craig (1975), Excess< sup> 3</sup> He in oceanic basalts: Evidence for terrestrial primordial helium, Earth and Planetary Science Letters, 26(2), 133139. Lupton, J., R. Weiss, and H. Craig (1977), Mantle helium in hydrothermal plumes in the Galapagos Rift, Nature, 267, 603-604. Lupton, J. E. (1983), Terrestrial inert gases-Isotope tracer studies and clues to primordial components in the mantle, Annual Review of Earth and Planetary Sciences, 11, 371-414. Lupton, J. E., and H. Craig (1981), A major helium-3 source at 15 S on the East Pacific Rise, Science, 214(4516), 13-18. Mamyrin, B., T. IN, G. Anufriev, and K. IL (1969), Anomalic isotopic composition of helium in volcanic gases, Doklady Akademii Nauk SSSR, 184(5), 1197. 61 Marty, B., A. Jambon, and Y. Sano (1989), Helium isotopes and CO2 in volcanic gases of Japan, Chemical Geology, 76(1-2), 25-40. Morrison, P., and J. Pine (1955), RADIOGENIC ORIGIN OF THE HELIUM ISOTOPES IN ROCK, Annals of the New York Academy of Sciences, 62(3), 7192. Padrón, E., et al. (2013), Diffusive helium emissions as a precursory sign of volcanic unrest, Geology, 41(5), 539-542. Ponko, S. C., and C. O. Sanders (1994), Inversion for P and S wave attenuation structure, Long Valley caldera, California, Journal of Geophysical Research: Solid Earth (1978-2012), 99(B2), 2619-2635. Sanford, W. E., R. G. Shropshire, and D. K. Solomon (1996), Dissolved gas tracers in groundwater: Simplified injection, sampling, and analysis, Water Resources Research, 32(6), 1635-1642. Sano, Y., Y. Nakamura, K. Notsu, and H. Wakita (1988), Influence of volcanic eruptions on helium isotope ratios in hydrothermal systems induced by volcanic eruptions, Geochimica et Cosmochimica Acta, 52(5), 1305-1308. Sano, Y., K. Sano, J.-i. Notsu, G. Ishibashi, H. Igarashi, and H. Wakita (1991), Secular variations in helium isotope ratios in an active volcano: Eruption and plug hypothesis, Earth and Planetary Science Letters, 107(1), 95-100. Sheldon, A. L. (2002), Diffusion of radiogenic helium in shallow groundwater: Implications for crustal degassing, Ph.D. thesis, 171-171 p. pp, The University of Utah, Ann Arbor. Solomon, D. K. (2000), 4He in groundwater, in Environmental Tracers in Subsurface Hydrology, edited, pp. 425-439, Springer. Sorey, M. L., B. M. Kennedy, W. C. Evans, C. D. Farrar, and G. A. Suemnicht (1993), Helium isotope and gas discharge variations associated with crustal unrest in Long Valley Caldera, California, 1989-1992, Journal of Geophysical Research: Solid Earth, 98(B9), 15871-15889. Sorey, M. L., W. C. Evans, B. M. Kennedy, C. D. Farrar, L. J. Hainsworth, and B. Hausback (1998), Carbon dioxide and helium emissions from a reservoir of magmatic gas beneath Mammoth Mountain, California, Journal of Geophysical Research: Solid Earth, 103(B7), 15303-15323. Stuer-Lauridsen, F. (2005), Review of passive accumulation devices for monitoring organic micropollutants in the aquatic environment, Environmental Pollution, 136(3), 503-524. 62 Torgersen, T., and W. Clarke (1985), Helium accumulation in groundwater, I: An evaluation of sources and the continental flux of crustal< sup> 4</sup> He in the Great Artesian Basin, Australia, Geochimica et Cosmochimica Acta, 49(5), 12111218. Wicks, C. W., D. Dzurisin, S. Ingebritsen, W. Thatcher, Z. Lu, and J. Iverson (2002), Magmatic activity beneath the quiescent Three Sisters volcanic center, central Oregon Cascade Range, USA, Geophysical Research Letters, 29(7), 26-21-26-24. Wilson, R. D., and D. M. Mackay (1995), A method for passive release of solutes from an unpumped well, Ground Water, 33(6), 936-945. Zhdanov, M. S. (2002), Preface, in Methods in Geochemistry and Geophysics, edited by S. Z. Michael, pp. XIX-XXIII, Elsevier.