## Publications

### REFEREED PUBLICATIONS IN JOURNALS

56. | A.M.M. Leal, D.A. Kulik, W.R. Smith, M.O. Saar An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling, Pure and Applied Chemistry, 89/5, pp. 597-643, 2017. Abstract We present an overview of novel numerical methods for chemical equilibrium and kinetic calculations for complex non-ideal multiphase systems. The methods we present for equilibrium calculations are based either on Gibbs energy minimization (GEM) calculations or on solving the system of extended law of mass-action (xLMA) equations. In both methods, no a posteriori phase stability tests, and thus no tentative addition or removal of phases during or at the end of the calculations, are necessary. All potentially stable phases are considered from the beginning of the calculation, and stability indices are immediately available at the end of the computation to determine which phases are actually stable at equilibrium. Both GEM and xLMA equilibrium methods are tailored for computationally demanding applications that require many rapid local equilibrium calculations, such as reactive transport modeling. The numerical method for chemical kinetic calculations we present supports both closed and open systems, and it considers a partial equilibrium simplification for fast reactions. The method employs an implicit integration scheme that improves stability and speed when solving the often stiff differential equations in kinetic calculations. As such, it requires compositional derivatives of the reaction rates to assemble the Jacobian matrix of the resultant implicit algebraic equations that are solved at every time step. We present a detailed procedure to calculate these derivatives, and we show how the partial equilibrium assumption affects their computation. These numerical methods have been implemented in Reaktoro (reaktoro.org), an open-source software for modeling chemically reactive systems. We finish with a discussion on the comparison of these methods with others in the literature. / Download |

55. | A.J. Luhmann, B.M. Tutolo, B.C. Bagley, D.F.R. Mildner, W.E. Seyfried Jr., M.O. Saar Permeability, porosity, and mineral surface area changes in basalt cores induced by reactive transport of CO2-rich brine, Water Resources Research, 53, pp. 1-20, 2017. Abstract Four reactive flow-through laboratory experiments (two each at 0.1 mL/min and 0.01 mL/min flow rates) at 150°C and 150 bar (15 MPa) are conducted on intact basalt cores to assess changes in porosity, permeability, and surface area caused by CO2-rich fluid-rock interaction. Permeability decreases slightly during the lower flow rate experiments and increases during the higher flow rate experiments. At the higher flow rate, core permeability increases by more than one order of magnitude in one experiment and less than a factor of two in the other due to differences in preexisting flow path structure. X-ray computed tomography (XRCT) scans of pre- and post-experiment cores identify both mineral dissolution and secondary mineralization, with a net decrease in XRCT porosity of ∼0.7%–0.8% for the larger pores in all four cores. (Ultra) small-angle neutron scattering ((U)SANS) data sets indicate an increase in both (U)SANS porosity and specific surface area (SSA) over the ∼1 nm to 10 µm scale range in post-experiment basalt samples, with differences due to flow rate and reaction time. Net porosity increases from summing porosity changes from XRCT and (U)SANS analyses are consistent with core mass decreases. (U)SANS data suggest an overall preservation of the pore structure with no change in mineral surface roughness from reaction, and the pore structure is unique in comparison to previously published basalt analyses. Together, these data sets illustrate changes in physical parameters that arise due to fluid-basalt interaction in relatively low pH environments with elevated CO2 concentration, with significant implications for flow, transport, and reaction through geologic formations. / Download |

54. | A.J. Luhmann, B.M. Tutolo, C. Tan, B.M. Moskowitz, M.O. Saar, W.E. Seyfried, Jr. Whole rock basalt alteration from CO2-rich brine during flow-through experiments at 150°C and 150 bar, Chemical Geology, 453, pp. 92-110, 2017. Abstract Four flow-through experiments at 150 °C were conducted on intact cores of basalt to assess alteration and mass transfer during reaction with CO2-rich fluid. Two experiments used a flow rate of 0.1 ml/min, and two used a flow rate of 0.01 ml/min. Permeability increased for both experiments at the higher flow rate, but decreased for the lower flow rate experiments. The experimental fluid (initial pH of 3.3) became enriched in Si, Mg, and Fe upon passing through the cores, primarily from olivine and titanomagnetite dissolution and possibly pyroxene dissolution. Secondary minerals enriched in Al and Si were present on post-experimental cores, and an Fe2O3-rich phase was identified on the downstream ends of the cores from the experiments at the lower flow rate. While we could not specifically identify if siderite (FeCO3) was present in the post-experimental basalt cores, siderite was generally saturated or supersaturated in outlet fluid samples, suggesting a thermodynamic drive for Fe carbonation from basalt-H2O-CO2 reaction. Reaction path models that employ dissolution kinetics of olivine, labradorite, and enstatite also suggest siderite formation at low pH. Furthermore, fluid-rock interaction caused a relatively high mobility of the alkali metals; up to 29% and 99% of the K and Cs present in the core, respectively, were preferentially dissolved from the cores, likely due to fractional crystallization effects that made alkali metals highly accessible. Together, these datasets illustrate changes in chemical parameters that arise due to fluid-basalt interaction in relatively low pH environments with elevated CO2. / Download |

53. | J.M. Myre, E. Frahm, D.J. Lilja, M.O. Saar TNT-NN: A Fast Active Set Method for Solving Large Non-Negative Least Squares Problems, (in press). no_download |

52. | A.M.M. Leal, D. Kulik, G. Kosakowski, M.O. Saar Computational methods for reactive transport modeling: An extended law of mass-action, xLMA, method for multiphase equilibrium calculations, Advances in Water Resources, 96, pp. 405-422, 2016. Abstract We present a numerical method for multiphase chemical equilibrium calculations based on a Gibbs energy minimization approach. The method can accurately and efficiently determine the stable phase assemblage at equilibrium independently of the type of phases and species that constitute the chemical system. We have successfully applied our chemical equilibrium algorithm in reactive transport simulations to demonstrate its effective use in computationally intensive applications. We used FEniCS to solve the governing partial differential equations of mass transport in porous media using finite element methods in unstructured meshes. Our equilibrium calculations were benchmarked with GEMS3K, the numerical kernel of the geochemical package GEMS. This allowed us to compare our results with a well-established Gibbs energy minimization algorithm, as well as their performance on every mesh node, at every time step of the transport simulation. The benchmark shows that our novel chemical equilibrium algorithm is accurate, robust, and efficient for reactive transport applications, and it is an improvement over the Gibbs energy minimization algorithm used in GEMS3K. The proposed chemical equilibrium method has been implemented in Reaktoro, a unified framework for modeling chemically reactive systems, which is now used as an alternative numerical kernel of GEMS. / Download |

51. | B.M. Tutolo, D.F. Mildner, C.V. Gagnon, M.O. Saar, W.E. Seyfried Nanoscale constraints on porosity generation and fluid flow during serpentinization, Geology, 44/2, pp. 103-106, 2016. Abstract Field samples of olivine-rich rocks are nearly always serpentinized—commonly to completion—but, paradoxically, their intrinsic porosity and permeability are diminishingly low. Serpentinization reactions occur through a coupled process of fluid infiltration, volumetric expansion, and reaction-driven fracturing. Pores and reactive surface area generated during this process are the primary pathways for fluid infiltration into and reaction with serpentinizing rocks, but the size and distribution of these pores and surface area have not yet been described. Here, we utilize neutron scattering techniques to present the first measurements of the evolution of pore size and specific surface area distribution in partially serpentinized rocks. Samples were obtained from the ca. 2 Ma Atlantis Massif oceanic core complex located off-axis of the Mid-Atlantic Ridge and an olivine-rich outcrop of the ca. 1.1 Ga Duluth Complex of the North American Mid-Continent Rift. Our measurements and analyses demonstrate that serpentine and accessory phases form with their own, inherent porosity, which accommodates the bulk of diffusive fluid flow during serpentinization and thereby permits continued serpentinization after voluminous serpentine minerals fill reaction-generated porosity. / Download |

50. | A.M.M. Leal, D.A. Kulik, M.O. Saar Enabling Gibbs energy minimization algorithms to use equilibrium constants of reactions in multiphase equilibrium calculations, Chemical Geology, 437, pp. 170-181, 2016. Abstract The geochemical literature provides numerous thermodynamic databases compiled from equilibrium constants of reactions. These databases are typically used in speciation calculations based on the law of mass action (LMA) approach. Unfortunately, such LMA databases cannot be directly used in equilibrium speciation methods based on the Gibbs energy minimization (GEM) approach because of their lack of standard chemical potentials of species. Therefore, we present in this work a simple conversion approach that calculates apparent standard chemical potentials of species from equilibrium constants of reactions. We assess the consistency and accuracy of the use of apparent standard chemical potentials in GEM algorithms by benchmarking equilibrium speciation calculations using GEM and LMA methods with the same LMA database. In all cases, we use PHREEQC to perform the LMA calculations, and we use its LMA databases to calculate the equilibrium constants of reactions. GEM calculations are performed using a Gibbs energy minimization method of Reaktoro — a unified open-source framework for numerical modeling of chemically reactive systems. By comparing the GEM and LMA results, we show that the use of apparent standard chemical potentials in GEM methods produces consistent and accurate equilibrium speciation results, thus validating our new, practical conversion technique that enables GEM algorithms to take advantage of many existing LMA databases, consequently extending and diversifying their range of applicability. / Download |

49. | T.A. Buscheck, J.M. Bielicki, T.A. Edmunds, Y. Hao, Y. Sun, J.B. Randolph, M.O. Saar Multifluid geo-energy systems: Using geologic CO2 storage for geothermal energy production and grid-scale energy storage in sedimentary basins, Geosphere, 12/3, pp. 1-19, 2016. Download |

48. | B.M. Tutolo, X.-Z. Kong, W.E. Seyfried Jr., M.O. Saar High performance reactive transport simulations examining the effects of thermal, hydraulic, and chemical (THC) gradients on fluid injectivity at carbonate CCUS reservoir scales, International Journal of Greenhouse Gas Control, 39, pp. 285-301, 2015. Abstract Carbonate minerals and CO2 are both considerably more soluble at low temperatures than they are at elevated temperatures. This inverse solubility has led a number of researchers to hypothesize that injecting low-temperature (i.e., less than the background reservoir temperature) CO2 into deep, saline reservoirs for CO2 Capture, Utilization, and Storage (CCUS) will dissolve CO2 and carbonate minerals near the injection well and subsequently exsolve and re-precipitate these phases as the fluids flow into the geothermally warm portion of the reservoir. In this study, we utilize high performance computing to examine the coupled effects of cool CO2 injection and background hydraulic head gradients on reservoir-scale mineral volume changes. We employ the fully coupled reactive transport simulator PFLOTRAN with calculations distributed over up to 800 processors to test 21 scenarios designed to represent a range of reservoir depths, hydraulic head gradients, and CO2 injection rates and temperatures. In the default simulations, 50 °C CO2 is injected at a rate of 50 kg/s into a 200 bar, 100 °C calcite or dolomite reservoir. By comparing these simulations with others run at varying conditions, we show that the effect of cool CO2 injection on reservoir-scale mineral volume changes tends to be relatively minor. We conclude that the low heat capacity of CO2 effectively prevents low-temperature CO2 injection from decreasing the temperature across large portions of the simulated carbonate reservoirs. This small thermal perturbation, combined with the low relative permeability of brine within the supercritical CO2 plume, yields limited dissolution and precipitation effects directly attributable to cool CO2 injection. Finally, we calculate that relatively high water-to-rock ratios, which may occur over much longer CCUS reservoir lifetimes or in materials with sufficiently high brine relative permeability within the supercritical CO2 plume, would be required to substantially affect injectivity through thermally-induced mineral dissolution and precipitation. Importantly, this study shows the utility of reservoir scale-reactive transport simulators for testing hypotheses and placing laboratory-scale observations into a CCUS reservoir-scale context. / Download |

47. | B.M. Tutolo, A.J. Luhmann, X.-Z. Kong, M.O. Saar, W.E. Seyfried Jr. CO2 sequestration in feldspar-rich sandstone: Coupled evolution of fluid chemistry, mineral reaction rates, and hydrogeochemical properties, Geochimica Et Cosmochimica Acta, 160, pp. 132-154, 2015. Abstract To investigate CO / Download_{2} Capture, Utilization, and Storage (CCUS) in sandstones, we performed three 150 °C flow-through experiments on K-feldspar-rich cores from the Eau Claire formation. By characterizing fluid and solid samples from these experiments using a suite of analytical techniques, we explored the coupled evolution of fluid chemistry, mineral reaction rates, and hydrogeochemical properties during CO_{2} sequestration in feldspar-rich sandstone. Overall, our results confirm predictions that the heightened acidity resulting from supercritical CO_{2} injection into feldspar-rich sandstone will dissolve primary feldspars and precipitate secondary aluminum minerals. A core through which CO_{2}-rich deionized water was recycled for 52 days decreased in bulk permeability, exhibited generally low porosity associated with high surface area in post-experiment core sub-samples, and produced an Al hydroxide secondary mineral, such as boehmite. However, two samples subjected to ?3 day single-pass experiments run with CO_{2}-rich, 0.94 mol/kg NaCl brines decreased in bulk permeability, showed generally elevated porosity associated with elevated surface area in post-experiment core sub-samples, and produced a phase with kaolinite-like stoichiometry. CO_{2}-induced metal mobilization during the experiments was relatively minor and likely related to Ca mineral dissolution. Based on the relatively rapid approach to equilibrium, the relatively slow near-equilibrium reaction rates, and the minor magnitudes of permeability changes in these experiments, we conclude that CCUS systems with projected lifetimes of several decades are geochemically feasible in the feldspar-rich sandstone end-member examined here. Additionally, the observation that K-feldspar dissolution rates calculated from our whole-rock experiments are in good agreement with literature parameterizations suggests that the latter can be utilized to model CCUS in K-feldspar-rich sandstone. Finally, by performing a number of reactive transport modeling experiments to explore processes occurring during the flow-through experiments, we have found that the overall progress of feldspar hydrolysis is negligibly affected by quartz dissolution, but significantly impacted by the rates of secondary mineral precipitation and their effect on feldspar saturation state. The observations produced here are critical to the development of models of CCUS operations, yet more work, particularly in the quantification of coupled dissolution and precipitation processes, will be required in order to produce models that can accurately predict the behavior of these systems. |

46. | N. Garapati, J.B. Randolph, M.O. Saar Brine displacement by CO2, energy extraction rates, and lifespan of a CO2-limited CO2-Plume Geothermal (CPG) system with a horizontal production well, Geothermics, 55, pp. 182-194, 2015. Abstract Several studies suggest that CO2-based geothermal energy systems may be operated economically when added to ongoing geologic CO2 sequestration. Alternatively, we demonstrate here that CO2-Plume Geothermal (CPG) systems may be operated long-term with a finite amount of CO2. We analyze the performance of such CO2-limited CPG systems as a function of various geologic and operational parameters. We find that the amount of CO2 required increases with reservoir depth, permeability, and well spacing and decreases with larger geothermal gradients. Furthermore, the onset of reservoir heat depletion decreases for increasing geothermal gradients and for both particularly shallow and deep reservoirs. / Download |

45. | B.M. Tutolo, A.T. Schaen, M.O. Saar, W.E. Seyfried Jr. Implications of the redissociation phenomenon for mineral-buffered fluids and aqueous species transport at elevated temperatures and pressures, Applied Geochemistry, 55, pp. 119-127, 2015. Abstract Aqueous species equilibrium constants and activity models form the foundation of the complex speciation codes used to model the geochemistry of geothermal energy production, extremophilic ecosystems, ore deposition, and a variety of other processes. Researchers have shown that a simple three species model (i.e., Na+, Cl?, and NaCl(aq)) can accurately describe conductivity measurements of concentrated NaCl and KCl solutions at elevated temperatures and pressures (Sharygin et al., 2002). In this model, activity coefficients of the charged species (e.g., Na+, K+, Cl?) become sufficiently low that the complexes must redisocciate with increasing salt concentration in order to meet equilibrium constant constraints. Redissociation decreases the proportion of the elements bound up as neutral complexes, and thereby increases the true ionic strength of the solution. In this contribution, we explore the consequences of the redissociation phenomenon in albite–paragonite–quartz (APQ) buffered systems. We focus on the implications of the redissociation phenomenon for mineral solubilities, particularly the observation that, at certain temperatures and pressures, calculated activities of charged ions in solution remain practically constant even as element concentrations increase from <1 molal to 4.5 molal. Finally, we note that redissociation has a similar effect on pH, and therefore aqueous speciation, in APQ-hosted systems. The calculations and discussion presented here are not limited to APQ-hosted systems, but additionally apply to many others in which the dominant cations and anions can form neutral complexes. / Download |

44. | B.M. Adams, T.H. Kuehn, J.M. Bielicki, J.B. Randolph, M.O. Saar A comparison of electric power output of CO2 Plume Geothermal (CPG) and brine geothermal systems for varying reservoir conditions, Applied Energy, 140, pp. 365-377, 2015. Abstract In contrast to conventional hydrothermal systems or enhanced geothermal systems, CO2 Plume Geothermal (CPG) systems generate electricity by using CO2 that has been geothermally heated due to sequestration in a sedimentary basin. Four CPG and two brine-based geothermal systems are modeled to estimate their power production for sedimentary basin reservoir depths between 1 and 5km, geothermal temperature gradients from 20 to 50°Ckm-1, reservoir permeabilities from 1×10-15 to 1×10-12m2 and well casing inner diameters from 0.14m to 0.41m. Results show that CPG direct-type systems produce more electricity than brine-based geothermal systems at depths between 2 and 3km, and at permeabilities between 10-14 and 10-13m2, often by a factor of two. This better performance of CPG is due to the low kinematic viscosity of CO2, relative to brine at those depths, and the strong thermosiphon effect generated by CO2. When CO2 is used instead of R245fa as the secondary working fluid in an organic Rankine cycle (ORC), the power production of both the CPG and the brine-reservoir system increases substantially; for example, by 22% and 20% for subsurface brine and CO2 systems, respectively, with a 35°Ckm-1 thermal gradient, 0.27m production and 0.41m injection well diameters, and 5×10-14m2 reservoir permeability. / Download |

43. | A.J. Luhmann, M. Covington, J. Myre, M. Perne, S.W. Jones, C.E. Alexander Jr., M.O. Saar Thermal damping and retardation in karst conduits, Hydrology and Earth System Sciences, 19/1, pp. 137-157, 2015. Abstract Water temperature is a non-conservative tracer in the environment. Variations in recharge temperature are damped and retarded as water moves through an aquifer due to heat exchange between water and rock. However, within karst aquifers, seasonal and short-term fluctuations in recharge temperature are often transmitted over long distances before they are fully damped. Using analytical solutions and numerical simulations, we develop relationships that describe the effect of flow path properties, flow-through time, recharge characteristics, and water and rock physical properties on the damping and retardation of thermal peaks/troughs in karst conduits. Using these relationships, one can estimate the thermal retardation and damping that would occur under given conditions with a given conduit geometry. Ultimately, these relationships can be used with thermal damping and retardation field data to estimate parameters such as conduit diameter. We also examine sets of numerical simulations where we relax some of the assumptions used to develop these relationships, testing the effects of variable diameter, variable velocity, open channels, and recharge shape on thermal damping and retardation to provide some constraints on uncertainty. Finally, we discuss a multitracer experiment that provides some field confirmation of our relationships. High temporal resolution water temperature data are required to obtain sufficient constraints on the magnitude and timing of thermal peaks and troughs in order to take full advantage of water temperature as a tracer. / Download |

42. | B.M. Tutolo, A.J. Luhmann, X.-Z. Kong, M.O. Saar, W.E. Seyfried Jr. Experimental observation of permeability changes in dolomite at CO2 sequestration conditions, Environmental Science and Technology, 48/4, pp. 2445-2452, 2014. Abstract Injection of cool CO2 into geothermally warm carbonate reservoirs for storage or geothermal energy production may lower near-well temperature and lead to mass transfer along flow paths leading away from the well. To investigate this process, a dolomite core was subjected to a 650 h, high pressure, CO2 saturated, flow-through experiment. Permeability increased from 10–15.9 to 10–15.2 m2 over the initial 216 h at 21 °C, decreased to 10–16.2 m2 over 289 h at 50 °C, largely due to thermally driven CO2 exsolution, and reached a final value of 10–16.4 m2 after 145 h at 100 °C due to continued exsolution and the onset of dolomite precipitation. Theoretical calculations show that CO2 exsolution results in a maximum pore space CO2 saturation of 0.5, and steady state relative permeabilities of CO2 and water on the order of 0.0065 and 0.1, respectively. Post-experiment imagery reveals matrix dissolution at low temperatures, and subsequent filling-in of flow passages at elevated temperature. Geochemical calculations indicate that reservoir fluids subjected to a thermal gradient may exsolve and precipitate up to 200 cm3 CO2 and 1.5 cm3 dolomite per kg of water, respectively, resulting in substantial porosity and permeability redistribution. / Download |

41. | A.J. Luhmann, X.-Z. Kong, B.M. Tutolo, N. Garapati, B.C. Bagley, M.O. Saar, W.E. Seyfried Jr. Experimental dissolution of dolomite by CO2-charged brine at 100oC and 150 bar: Evolution of porosity, permeability, and reactive surface area, Chemical Geology, 380, pp. 145-160, 2014. Abstract Hydrothermal flow experiments of single-pass injection of CO2-charged brine were conducted on nine dolomite cores to examine fluid–rock reactions in dolomite reservoirs under geologic carbon sequestration conditions. Post-experimental X-ray computed tomography (XRCT) analysis illustrates a range of dissolution patterns, and significant increases in core bulk permeability were measured as the dolomite dissolved. Outflow fluids were below dolomite saturation, and cation concentrations decreased with time due to reductions in reactive surface area with reaction progress. To determine changes in reactive surface area, we employ a power-law relationship between reactive surface area and porosity (Luquot and Gouze, 2009). The exponent in this relationship is interpreted to be a geometrical parameter that controls the degree of surface area change per change in core porosity. Combined with XRCT reconstructions of dissolution patterns, we demonstrate that this exponent is inversely related to both the flow path diameter and tortuosity of the dissolution channel. Even though XRCT reconstructions illustrate dissolution at selected regions within each core, relatively high Ba and Mn recoveries in fluid samples suggest that dissolution occurred along the core’s entire length and width. Analysis of porosity–permeability data indicates an increase in the rate of permeability enhancement per increase in porosity with reaction progress as dissolution channels lengthen along the core. Finally, we incorporate the surface area–porosity model of Luquot and Gouze (2009) with our experimentally fit parameters into TOUGHREACT to simulate experimental observations. / Download |

40. | B.M. Adams, T.H. Kuehn, J.M. Bielicki, J.B. Randolph, M.O. Saar On the importance of the thermosiphon effect in CPG (CO2-Plume geothermal) power systems, Energy, 69, pp. 409-418, 2014. Abstract CPG (CO2 Plume Geothermal) energy systems use CO2 to extract thermal energy from naturally permeable geologic formations at depth. CO2 has advantages over brine: high mobility, low solubility of amorphous silica, and higher density sensitivity to temperature. The density of CO2 changes substantially between geothermal reservoir and surface plant, resulting in a buoyancy-driven convective current – a thermosiphon – that reduces or eliminates pumping requirements. We estimated and compared the strength of this thermosiphon for CO2 and for 20 weight percent NaCl brine for reservoir depths up to 5 km and geothermal gradients of 20, 35, and 50 °C/km. We found that through the reservoir, CO2 has a pressure drop approximately 3–12 times less than brine at the same mass flowrate, making the CO2 thermosiphon sufficient to produce power using reservoirs as shallow as 0.5 km. At 2.5 km depth with a 35 °C/km gradient – the approximate western U.S. continental mean – the CO2 thermosiphon converted approximately 10% of the energy extracted from the reservoir to fluid circulation, compared to less than 1% with brine, where additional mechanical pumping is necessary. We found CO2 is a particularly advantageous working fluid at depths between 0.5 km and 3 km. / Download |

39. | B.M. Tutolo, X.-Z. Kong, W.E. Seyfried Jr., M.O. Saar Internal consistency in aqueous geochemical data revisited: Applications to the aluminum system, Geochimica et Cosmochimica Acta, 133, pp. 216-234, 2014. Abstract Internal consistency of thermodynamic data has long been considered vital for confident calculations of aqueous geochemical processes. However, an internally consistent mineral thermodynamic data set is not necessarily consistent with calculations of aqueous species thermodynamic properties due, potentially, to improper or inconsistent constraints used in the derivation process. In this study, we attempt to accommodate the need for a mineral thermodynamic data set that is internally consistent with respect to aqueous species thermodynamic properties by adapting the least squares optimization methods of Powell and Holland (1985). This adapted method allows for both the derivation of mineral thermodynamic properties from fluid chemistry measurements of solutions in equilibrium with mineral assemblages, as well as estimates of the uncertainty on the derived results. Using a large number of phase equilibria, solubility, and calorimetric measurements, we have developed a thermodynamic data set of 12 key aluminum-bearing mineral phases. These data are derived to be consistent with Na+ and K+ speciation data presented by Shock and Helgeson (1988), H4SiO4(aq) data presented by Stefánsson (2001), and the Al speciation data set presented by Tagirov and Schott (2001). Many of the constraining phase equilibrium measurements are exactly the same as those used to develop other thermodynamic data, yet our derived values tend to be quite different than some of the others’ due to our choices of reference data. The differing values of mineral thermodynamic properties have implications for calculations of Al mineral solubilities; specifically, kaolinite solubilities calculated with the developed data set are as much as 6.75 times lower and 73% greater than those calculated with Helgeson et al. (1978) and Holland and Powell (2011) data, respectively. Where possible, calculations and experimental data are compared at low T, and the disagreement between the two sources reiterates the common assertion that low-T measurements of phase equilibria and mineral solubilities in the aluminum system rarely represent equilibrium between water and well-crystallized, aluminum-bearing minerals. As an ancillary benefit of the derived data, we show that it may be combined with high precision measurements of aqueous complex association constants to derive neutral species activity coefficients in supercritical fluids. Although this contribution is specific to the aluminum system, the methods and concepts developed here can help to improve the calculation of water–rock interactions in a broad range of earth systems. / Download |

38. | N. Garapati, J.B. Randolph, J.L. Valencia Jr., M.O. Saar CO2-Plume Geothermal (CPG) Heat Extraction in Multi-layered Geologic Reservoirs, Energy Procedia, 63, pp. 7631-7643, 2014. Abstract CO2-Plume Geothermal (CPG) technology involves injecting CO2 into natural, highly permeable geologic units to extract energy. The subsurface CO2 absorbs heat from the reservoir, buoyantly rises to the surface, and drives a power generation system. The CO2 is then cooled and reinjected underground. Here, we analyze the effects of multi-layered geologic reservoirs on CPG system performance by examining the CO2 mass fraction in the produced fluid, pore-fluid pressure buildup during operation, and heat energy extraction rates. The produced CO2 mass fraction depends on the stratigraphic positions of highly permeable layers which also affect the pore-fluid pressure drop across the reservoir. / Download |

37. | T.A. Buscheck, J.M. Bielicki, M. Chen, Y. Sun, Y. Hao, T.A. Edmunds, M.O. Saar, J.B. Randolph Integrating CO2 Storage with Geothermal Resources for Dispatchable Renewable Electricity, Energy procedia, 63, pp. 7619-7630, 2014. Abstract We present an approach that uses the huge fluid and thermal storage capacity of the subsurface, together with geologic CO2 storage, to harvest, store, and dispatch energy from subsurface (geothermal) and surface (solar, nuclear, fossil) thermal resources, as well as energy from electrical grids. Captured CO2 is injected into saline aquifers to store pressure, generate artesian flow of brine, and provide an additional working fluid for efficient heat extraction and power conversion. Concentric rings of injection and production wells are used to create a hydraulic divide to store pressure, CO2, and thermal energy. Such storage can take excess power from the grid and excess/waste thermal energy, and dispatch that energy when it is demanded, enabling increased penetration of variable renewables. Stored CO2 functions as a cushion gas to provide enormous pressure-storage capacity and displaces large quantities of brine, which can be desalinated and/or treated for a variety of beneficial uses. Geothermal power and energy-storage applications may generate enough revenues to justify CO2 capture costs. / Download |

36. | B.M. Adams, T.H. Kuehn, J.B. Randolph, M.O. Saar The reduced pumping power requirements from increasing the injection well fluid density, Transactions – Geothermal Resources Council, 37, pp. 667-672, 2013. no_download |

35. | J.B. Randolph, M.O. Saar, J.M. Bielicki Geothermal energy production at geologic CO2 sequestration sites: Impact of thermal drawdown on reservoir pressure, Energy Procedia, 37, pp. 6625-6635, 2013. Abstract Recent geotechnical research shows that geothermal heat can be efficiently mined by circulating carbon dioxide through naturally permeable rock formations — a method called CO2 Plume Geothermal — the same geologic reservoirs that are suitable for deep saline aquifer CO2 sequestration or enhanced oil recovery. This paper describes the effect of thermal drawdown on reservoir pressure buildup during sequestration operations, revealing that geothermal heat mining can decrease overpressurization by 10% or more. / Download |

34. | R. Gottardi, P.-H. Kao, M.O. Saar, Ch. Teyssier Effects of permeability fields on fluid, heat, and oxygen isotope transport in extensional detachment systems, Geochemistry, Geophysics, Geosystems, 14/5, pp. 1493-1522, 2013. Abstract [1] Field studies of Cordilleran metamorphic core complexes indicate that meteoric fluids permeated the upper crust down to the detachment shear zone and interacted with highly deformed and recrystallized (mylonitic) rocks. The presence of fluids in the brittle/ductile transition zone is recorded in the oxygen and hydrogen stable isotope compositions of the mylonites and may play an important role in the thermomechanical evolution of the detachment shear zone. Geochemical data show that fluid flow in the brittle upper crust is primarily controlled by the large-scale fault-zone architecture. We conduct continuum-scale (i.e., large-scale, partial bounce-back) lattice-Boltzman fluid, heat, and oxygen isotope transport simulations of an idealized cross section of a metamorphic core complex. The simulations investigate the effects of crust and fault permeability fields as well as buoyancy-driven flow on two-way coupled fluid and heat transfer and resultant exchange of oxygen isotopes between meteoric fluid and rock. Results show that fluid migration to middle to lower crustal levels is fault controlled and depends primarily on the permeability contrast between the fault zone and the crustal rocks. High fault/crust permeability ratios lead to channelized flow in the fault and shear zones, while lower ratios allow leakage of the fluids from the fault into the crust. Buoyancy affects mainly flow patterns (more upward directed) and, to a lesser extent, temperature distributions (disturbance of the geothermal field by ~25°C). Channelized fluid flow in the shear zone leads to strong vertical and horizontal thermal gradients, comparable to field observations. The oxygen isotope results show δ18O depletion concentrated along the fault and shear zones, similar to field data. / Download |

33. | S.D.C. Walsh, M.O. Saar Developing extensible lattice-Boltzmann simulators for general-purpose graphics-processing units, Communications in Computational Physics, 13/3, pp. 867-879, 2013. Abstract Lattice-Boltzmann methods are versatile numerical modeling techniques capable of reproducing a wide variety of fluid-mechanical behavior. These methods are well suited to parallel implementation, particularly on the single-instruction multiple data (SIMD) parallel processing environments found in computer graphics processing units (GPUs). Although recent programming tools dramatically improve the ease with which GPUbased applications can be written, the programming environment still lacks the flexibility available to more traditional CPU programs. In particular, it may be difficult to develop modular and extensible programs that require variable on-device functionality with current GPU architectures. This paper describes a process of automatic code generation that overcomes these difficulties for lattice-Boltzmann simulations. It details the development of GPU-based modules for an extensible lattice-Boltzmann simulation package – LBHydra. The performance of the automatically generated code is compared to equivalent purposewritten codes for both single-phase,multiphase, andmulticomponent flows. The flexibility of the new method is demonstrated by simulating a rising, dissolving droplet moving through a porous medium with user generated lattice-Boltzmann models and subroutines. / Download |

32. | X.-Z. Kong, M.O. Saar DBCreate: A SUPCRT92-based program for producing EQ3/6, TOUGHREACT, and GWB thermodynamic databases at user-defined T and P, Computers and Geosciences, 51, pp. 415-417, 2013. Abstract SUPCRT92 is a widely used software package for calculating the standard thermodynamic properties of minerals, gases, aqueous species, and reactions. However, it is labor-intensive and error-prone to use it directly to produce databases for geochemical modeling programs such as EQ3/6, the Geochemist’s Workbench, and TOUGHREACT. DBCreate is a SUPCRT92-based software program written in FORTRAN90/95 and was developed in order to produce the required databases for these programs in a rapid and convenient way. This paper describes the overall structure of the program and provides detailed usage instructions. / Download |

31. | X.-Z. Kong, M.O. Saar Numerical study of the effects of permeability heterogeneity on density-driven convective mixing during CO2 dissolution storage, Int. J. Greenhouse Gas Control, 19, pp. 160-173, 2013. Abstract Permanence and security of carbon dioxide (CO2) in geologic formations requires dissolution of CO2 into brine, which slightly increases the brine density. Previous studies have shown that this small increase in brine density induces convective currents, which greatly enhances the mixing efficiency and thus CO2 storage capacity and rate in the brine. Density-driven convection, in turn, is known to be largely dominated by permeability heterogeneity. This study explores the relationship between the process of density-driven convection and the permeability heterogeneity of an aquifer during CO2 dissolution storage, using high-resolution numerical simulations. While the porosity is kept constant, the heterogeneity of the aquifer is introduced through a spatially varying permeability field, characterized by the Dykstra-Parsons coefficient and the correlation length. Depending on the concentration profile of dissolved CO2, we classify the convective finger patterns as dispersive, preferential, and unbiased fingering. Our results indicate that the transition between unbiased and both preferential and dispersive fingering is mainly governed by the Dykstra-Parsons coefficient, whereas the transition between preferential and dispersive fingering is controlled by the permeability correlation length. Furthermore, we find that the CO2 dissolution flux at the top boundary will reach a time-independent steady state. Although this flux strongly correlates with permeability distribution, it generally increases with the permeability heterogeneity when the correlation length is less than the system size. / Download |

30. | A.J. Luhmann, X.-Z. Kong, B.M. Tutolo, K. Ding, M.O. Saar, W.E. Seyfried Jr. Permeability reduction produced by grain reorganization and accumulation of exsolved CO2 during geologic carbon sequestration: A new CO2 trapping mechanism, Environmental Science and Technlogy, 47/1, pp. 242-251, 2013. Abstract Carbon sequestration experiments were conducted on uncemented sediment and lithified rock from the Eau Claire Formation, which consisted primarily of K-feldspar and quartz. Cores were heated to accentuate reactivity between fluid and mineral grains and to force CO2 exsolution. Measured permeability of one sediment core ultimately reduced by 4 orders of magnitude as it was incrementally heated from 21 to 150 °C. Water-rock interaction produced some alteration, yielding sub-?m clay precipitation on K-feldspar grains in the core’s upstream end. Experimental results also revealed abundant newly formed pore space in regions of the core, and in some cases pores that were several times larger than the average grain size of the sediment. These large pores likely formed from elevated localized pressure caused by rapid CO2 exsolution within the core and/or an accumulating CO2 phase capable of pushing out surrounding sediment. CO2 filled the pores and blocked flow pathways. Comparison with a similar experiment using a solid arkose core indicates that CO2 accumulation and grain reorganization mainly contributed to permeability reduction during the heated sediment core experiment. This suggests that CO2 injection into sediments may store more CO2 and cause additional permeability reduction than is possible in lithified rock due to grain reorganization. / Download |

29. | M. Covington, A.F. Banwell, J. Gulley, M.O. Saar Quantifying the effects of glacier conduit geometry and recharge on proglacial hydrograph form, Journal of Hydrology, 414-415, pp. 59-71, 2012. Abstract The configuration of glacier hydrological systems is often inferred from proxy data, such as hydrographs, that are collected in proglacial streams. Seasonal changes in the peakedness of hydrographs are thought to reflect changes in the configuration of the subglacial drainage system. However, the amount of information that proglacial hydrographs contain about drainage system configurations depends critically on the degree to which the drainage systems modify recharge hydrographs. If the drainage system does not modify recharge hydrographs, then proglacial hydrographs primarily reflect the recharge conditions produced by supraglacial inputs. Here, we develop a theoretical framework to determine the circumstances under which glacier drainage systems can modify recharge hydrographs and the circumstances under which recharge pulses pass through glaciers unchanged. We address the capability of single conduits, simple arborescent conduit networks, and linked cavity systems to modify diurnal recharge pulses. Simulations of discharge through large sets of such systems demonstrate that, unless large reservoirs or significant constrictions are present, the discharge hydrographs of simple glacial conduit systems are nearly identical to their recharge hydrographs. Conduit systems tend not to modify hydrographs because the changes in storage within englacial and subglacial conduit networks on short time scales are typically small compared to their ability to transmit water. This finding suggests that proglacial hydrographs reflect a variety of factors, including surface melt rate, surface water transfer, and subglacial water transfer. In many cases the influence of suglacial processes may be relatively minor. As a result, the evolution of proglacial hydrographs cannot be used unambiguously to infer changes in the structure or efficiency of englacial or subglacial hydrological systems, without accurate knowledge of the nature of the recharge hydrograph driving the flow. / Download |

28. | J.B. Randolph, B.M. Adams, T.H. Kuehn, M.O. Saar Wellbore heat transfer in CO2-based geothermal systems, Geothermal Resources Council (GRC) Transactions, 36, pp. 549-554, 2012. no_download |

27. | M.D. Covington, A.J. Luhmann, C. Wicks, M.O. Saar Process length scales and longitudinal damping in karst conduits, Journal of Geophysical Research – Earth Surface, 117, F01025, 2012. Abstract [1] Simple mathematical models often allow an intuitive grasp of the function of physical systems. We develop a mathematical framework to investigate reactive or dissipative transport processes within karst conduits. Specifically, we note that for processes that occur within a characteristic timescale, advection along the conduit produces a characteristic process length scale. We calculate characteristic length scales for the propagation of thermal and electrical conductivity signals along karst conduits. These process lengths provide a quantitative connection between karst conduit geometry and the signals observed at a karst spring. We show that water input from the porous/fractured matrix is also characterized by a length scale and derive an approximation that accounts for the influence of matrix flow on the transmission of signals through the aquifer. The single conduit model is then extended to account for conduits with changing geometries and conduit flow networks, demonstrating how these concepts can be applied in more realistic conduit geometries. We introduce a recharge density function, ϕR, which determines the capability of an aquifer to damp a given signal, and cast previous explanations of spring variability within this framework. Process lengths are a general feature of karst conduits and surface streams, and we conclude with a discussion of other potential applications of this conceptual and mathematical framework. / Download |

26. | S.C. Alexander, M.O. Saar Improved characterization of small u for Jacob pumping test analysis methods, Ground Water, 50/2, pp. 256-265, 2012. Abstract Numerous refinements have been proposed to traditional pumping test analyses, yet many hydrogeologists continue to use the Jacob method due to its simplicity. Recent research favors hydraulic tomography and inverse numerical modeling of pumping test data. However, at sites with few wells, or relatively short screens, the data requirements of these methods may be impractical within physical and fiscal constraints. Alternatively, an improved understanding of the assumptions and limitations of Theis and, due to their widespread usage, Jacob analyses, leads to improved interpretations in data-poor environments. A fundamental requirement of Jacob is a “small” value of u = f(r2/t), with radial distance, r, and pumping time, t. However, selection of a too stringent (i.e., too low) maximum permissible u-value, umax, results in rejection of usable data from wells beyond a maximum radius, rmax. Conversely, data from small radii, less than rmin, where turbulent- and vertical-flow components arise, can result in acceptance of inappropriate data. Usage of drawdown data from wells too close to the pumping well, and exclusion of data from wells deemed too far, can cause unrealistic aquifer transmissivity, permeability, and storativity determinations. Here, data from an extensive well field in a glacial-outwash aquifer in north-central Minnesota, USA, are used to develop a new estimate for umax. Traditionally quoted values for umax range from 0.01 to 0.05. Our proposed value for Jacob distance-drawdown analyses is significantly higher with umax up to 0.2, resulting in larger allowable rmax-values and a higher likelihood of inclusion of additional wells in such pumping test analyses. / Download |

25. | M.D. Covington, A.J. Luhmann, F. Gabrovsek, M.O. Saar, I. Willis, C.M. Wicks Mechanisms of heat exchange between water and rock in karst conduits, Water Resources Research, 47, W10514/10, 2011. Abstract [1] Previous studies, motivated by understanding water quality, have explored the mechanisms for heat transport and heat exchange in surface streams. In karst aquifers, temperature signals play an additional important role since they carry information about internal aquifer structures. Models for heat transport in karst conduits have previously been developed; however, these models make different, sometimes contradictory, assumptions. Additionally, previous models of heat transport in karst conduits have not been validated using field data from conduits with known geometries. Here we use analytical solutions of heat transfer to examine the relative importance of heat exchange mechanisms and the validity of the assumptions made by previous models. The relative importance of convection, conduction, and radiation is a function of time. Using a characteristic timescale, we show that models neglecting rock conduction produce spurious results in realistic cases. In contrast to the behavior of surface streams, where conduction is often negligible, conduction through the rock surrounding a conduit determines heat flux at timescales of weeks and longer. In open channel conduits, radiative heat flux can be significant. In contrast, convective heat exchange through the conduit air is often negligible. Using the rules derived from our analytical analysis, we develop a numerical model for heat transport in a karst conduit. Our model compares favorably to thermal responses observed in two different karst settings: a cave stream fed via autogenic recharge during a snowmelt event, and an allogenically recharged cave stream that experiences continuous temperature fluctuations on many timescales. / Download |

24. | J.B. Randolph, M.O. Saar Impact of reservoir permeability on the choice of subsurface geothermal heat exchange fluid: CO2 versus water and native brine, Geothermal Resources Council (GRC) Transactions, 35, pp. 521-526, 2011. no_download |

23. | J.B. Randolph, M.O. Saar Combining geothermal energy capture with geologic carbon dioxide sequestration, Geophysical Research Letters, 38, L10401, 2011. Abstract [1] Geothermal energy offers clean, renewable, reliable electric power with no need for grid-scale energy storage, yet its use has been constrained to the few locations worldwide with naturally high geothermal heat resources and groundwater availability. We present a novel approach with the potential to permit expansion of geothermal energy utilization: heat extraction from naturally porous, permeable formations with CO2 as the injected subsurface working fluid. Fluid-mechanical simulations reveal that the significantly higher mobility of CO2, compared to water, at the temperature/pressure conditions of interest makes CO2 an attractive heat exchange fluid. We show numerically that, compared to conventional water-based and engineered geothermal systems, the proposed approach provides up to factors of 2.9 and 5.0, respectively, higher geothermal heat energy extraction rates. Consequently, more regions worldwide could be economically used for geothermal electricity production. Furthermore, as the injected CO2 is eventually geologically sequestered, such power plants would have negative carbon footprints. / Download |

22. | M.A. Davis, S.D.C. Walsh, M.O. Saar Statistically reconstructing continuous isotropic and anisotropic two-phase media while preserving macroscopic material properties, Physical Review E, 83, 026706, 2011. Abstract We propose a method to generate statistically similar reconstructions of two-phase media. As with previous work, we initially characterize the microstructure of the material using two-point correlation functions (a subset of spatial correlation functions) and then generate numerical reconstructions using a simulated annealing method that preserves the geometric relationships of the material’s phase of interest. However, in contrast to earlier contributions that consider reconstructions composed of discrete arrays of pixels or voxels alone, we generate reconstructions based on assemblies of continuous, three-dimensional, interpenetrating objects. The result is a continuum description of the material microstructure (as opposed to a discretized or pixelated description), capable of efficiently representing large disparities in scale. Different reconstruction methods are considered based on distinct combinations of two-point correlation functions of varying degrees of complexity. The quality of the reconstruction methods are evaluated by comparing the total pore fraction, specific surface area of the percolating cluster, pore fraction of the percolating cluster, tortuosity, and permeability of the reconstructions to those of a set of reference assemblies. Elsewhere it has been proposed that two-phase media could be statistically reproduced with only two spatial correlation functions: the two-point probability function (the probability that two points lie within the same phase) and the lineal path function (the probability that a line between two points lies entirely within the same phase). We find that methods employing the two-point probability function and lineal path function are improved if the percolating cluster volume is also considered in the reconstruction. However, to reproduce more complicated geometric assemblies, we find it necessary to employ the two-point probability, two-point cluster, and lineal path function in addition to the percolating cluster volume to produce a generally accurate statistical reconstruction. / Download |

21. | J.B. Randolph, M.O. Saar Coupling carbon dioxide sequestration with geothermal energy capture in naturally permeable, porous geologic formations: Implications for CO2 sequestration, Energy Procedia, 4, pp. 2206-2213, 2011. Abstract Carbon dioxide (CO2) sequestration in deep saline aquifers and exhausted oil and natural gas fields has been widely considered as a means for reducing CO2 emissions to the atmosphere as a counter-measure to global warming. However, rather than treating CO2 merely as a waste fluid in need of permanent disposal, we propose that it could also be used as a working fluid in geothermal energy capture, as its thermodynamic and fluid mechanical properties suggest it transfers geothermal heat more efficiently than water. Energy production and sales in conjunction with sequestration would improve the economic viability of CO2 sequestration, a critical challenge for large-scale implementation of the technology. In addition, using CO2 as the working fluid in geothermal power systems may permit utilization of lower temperature geologic formations than those that are currently deemed economically viable, leading to more widespread utilization of geothermal energy. Here, we present the results of early-stage calculations demonstrating the geothermal energy capture potential of CO2-based geothermal systems and implications of such energy capture for the economic viability of geologic CO2 sequestration. / Download |

20. | M.O. Saar Review: Geothermal heat as a tracer of large-scale groundwater flow and as a means to determine permeability fields, special theme issue on Environmental Tracers and Groundwater Flow, editor-invited peer-reviewed contribution, Hydrogeology Journal, 19, pp. 31-52, 2011. Abstract A review of coupled groundwater and heat transfer theory is followed by an introduction to geothermal measurement techniques. Thereafter, temperature-depth profiles (geotherms) and heat discharge at springs to infer hydraulic parameters and processes are discussed. Several studies included in this review state that minimum permeabilities of approximately 5 × 10−17 < kmin <10−15 m2 are required to observe advective heat transfer and resultant geotherm perturbations. Permeabilities below kmin tend to cause heat-conduction-dominated systems, precluding inversion of temperature fields for groundwater flow patterns and constraint of permeabilities other than being / Download |

19. | J. Myre, S.D.C. Walsh, D.J. Lilja, M.O. Saar Performance analysis of single-phase multiphase, and multicomponent lattice-Boltzmann fluid flow simulations on GPU clusters, Concurrency and Computation: Practice and Experience, 23, pp. 332-350, 2010. Abstract Abstract The lattice-Boltzmann method is well suited for implementation in single-instruction multiple-data (SIMD) environments provided by general purpose graphics processing units (GPGPUs). This paper discusses the integration of these GPGPU programs with OpenMP to create lattice-Boltzmann applications for multi-GPU clusters. In addition to the standard single-phase single-component lattice-Boltzmann method, the performances of more complex multiphase, multicomponent models are also examined. The contributions of various GPU lattice-Boltzmann parameters to the performance are examined and quantified with a statistical model of the performance using Analysis of Variance (ANOVA). By examining single- and multi-GPU lattice-Boltzmann simulations with ANOVA, we show that all the lattice-Boltzmann simulations primarily depend on effects corresponding to simulation geometry and decomposition, and not on the architectural aspects of GPU. Additionally, using ANOVA we confirm that the metrics of Efficiency and Utilization are not suitable for memory-bandwidth-dependent codes. Copyright © 2010 John Wiley & Sons, Ltd. / Download |

18. | S. Dasgupta, M.O. Saar, R.L. Edwards, C.-C. Shen, H. Cheng, C.E. Alexander Jr. Three thousand years of extreme rainfall events recorded in stalagmites from Spring Valley Caverns, Minnesota, Earth and Planetary Science Letters, 300, pp. 46-54, 2010. Abstract Annual layer analysis in two stalagmites collected from Spring Valley Caverns, southeastern Minnesota, reveals hydrological response of the cave to extreme rainfall events in the Midwest, USA. Cave-flooding events are identified within the two samples by the presence of detrital layers composed of clay sized particles. Comparison with instrumental records of precipitation demonstrates a strong correlation between these cave-flood events and extreme rainfall observed in the Upper Mississippi Valley. A simple model is developed to assess the nature of rainfall capable of flooding the cave. The model is first calibrated to the last 50-yr (1950–1998 A.D.) instrumental record of daily precipitation data for the town of Spring Valley and verified with the first 50 yr of record from 1900 to 1949 A.D. Frequency analysis shows that these extreme flood events have increased from the last half of the nineteenth century. Comparison with other paleohydrological records shows increased occurrence of extreme rain events during periods of higher moisture availability. Our study implies that increased moisture availability in the Midwestern region, due to rise in temperature from global warming could lead to an increase in the occurrence of extreme rainfall events. / Download |

17. | S.D.C. Walsh, M.O. Saar Interpolated lattice-Boltzmann boundary conditions for surface reaction kinetics, Physical Review E, 82, 066703, 2010. Abstract This paper describes a method for implementing surface reaction kinetics in lattice Boltzmann simulations. The interpolated boundary conditions are capable of simulating surface reactions and dissolution at both stationary and moving solid-fluid and fluid-fluid interfaces. Results obtained with the boundary conditions are compared to analytical solutions for first-order and constant-flux kinetic surface reactions in a one-dimensional half space, as well as to the analytical solution for evaporation from the surface of a cylinder. Excellent agreement between analytical and simulated results is obtained for a wide range of diffusivities, lattice velocities, and surface reaction rates. The boundary model’s ability to represent dissolution in binary fluid mixtures is demonstrated by modeling diffusion from a rising bubble and dissolution of a droplet near a flat plate. / Download |

16. | J.B. Randolph, M.O. Saar Coupling geothermal energy capture with carbon dioxide sequestration in naturally permeable, porous geologic formations: A comparison with enhanced geothermal systems, Geothermal Resources Council (GRC) Transactions, 34, pp. 433-438, 2010. no_download |

15. | S.D.C. Walsh, M.O. Saar Macroscale lattice-Boltzmann methods for low-Peclet-number solute and heat transport in heterogeneous porous media, Water Resources Research, 46, W07517, 2010. Abstract [1] This paper introduces new methods for simulating subsurface solute and heat transport in heterogeneous media using large-scale lattice-Boltzmann models capable of representing both macroscopically averaged porous media and open channel flows. Previous examples of macroscopically averaged lattice-Boltzmann models for solute and heat transport are only applicable to homogeneous media. Here, we extend these models to properly account for heterogeneous pore-space distributions. For simplicity, in the majority of this paper we assume low Peclet number flows with an isotropic dispersion tensor. Nevertheless, this approach may also be extended to include anisotropic-dispersion by using multiple relaxation time lattice-Boltzmann methods. We describe two methods for introducing heterogeneity into macroscopically averaged lattice-Boltzmann models. The first model delivers the desired behavior by introducing an additional time-derivative term to the collision rule; the second model by separately weighting symmetric and anti-symmetric components of the fluid packet densities. Chapman-Enskog expansions are conducted on the governing equations of the two models, demonstrating that the correct constitutive behavior is obtained in both cases. In addition, methods for improving model stability at low porosities are also discussed: (1) an implicit formulation of the model; and (2) a local transformation that normalizes the lattice-Boltzmann model by the local porosity. The model performances are evaluated through comparisons of simulated results with analytical solutions for one- and two-dimensional flows, and by comparing model predictions to finite element simulations of advection isotropic-dispersion in heterogeneous porous media. We conclude by presenting an example application, demonstrating the ability of the new models to couple with simulations of reactive flow and changing flow geometry: a simulation of groundwater flow through a carbonate system. / Download |

14. | M. Covington, C.M. Wicks, M.O. Saar A dimensionless number describing the effects of recharge and geometry on discharge from simple karst aquifers, Water Resources Research, 45, W11410, 2009. Abstract [1] The responses of karstic aquifers to storms are often used to obtain information about aquifer geometry. In general, spring hydrographs are a function of both system geometry and recharge. However, the majority of prior work on storm pulses through karst has not studied the effect of recharge on spring hydrographs. To examine the relative importance of geometry and recharge, we break karstic aquifers into elements according to the manner of their response to transient flow and demonstrate that each element has a characteristic response timescale. These fundamental elements are full pipes, open channels, reservoir/constrictions, and the porous matrix. Taking the ratio of the element timescale with the recharge timescale produces a dimensionless number, γ, that is used to characterize aquifer response to a storm event. Using sets of simulations run with randomly selected element parameters, we demonstrate that each element type has a critical value of γ below which the shape of the spring hydrograph is dominated by the shape of the recharge hydrograph and above which the spring hydrograph is significantly modified by the system geometry. This allows separation of particular element/storm pairs into recharge-dominated and geometry-dominated regimes. While most real karstic aquifers are complex combinations of these elements, we draw examples from several karst systems that can be represented by single elements. These examples demonstrate that for real karstic aquifers full pipe and open channel elements are generally in the recharge-dominated regime, whereas reservoir/constriction elements can fall in either the recharge- or geometry-dominated regimes. / Download |

13. | S.D.C. Walsh, M.O. Saar, P. Bailey, D.J. Lilja Accelerating geoscience and engineering system simulations on graphics hardware, Computers and Geosciences, 35/12, pp. 2353-2364, 2009. Abstract Many complex natural systems studied in the geosciences are characterized by simple local-scale interactions that result in complex emergent behavior. Simulations of these systems, often implemented in parallel using standard central processing unit (CPU) clusters, may be better suited to parallel processing environments with large numbers of simple processors. Such an environment is found in graphics processing units (GPUs) on graphics cards. This paper discusses GPU implementations of three example applications from computational fluid dynamics, seismic wave propagation, and rock magnetism. These candidate applications involve important numerical modeling techniques, widely employed in physical system simulations, that are themselves examples of distinct computing classes identified as fundamental to scientific and engineering computing. The presented numerical methods (and respective computing classes they belong to) are: (1) a lattice-Boltzmann code for geofluid dynamics (structured grid class); (2) a spectral-finite-element code for seismic wave propagation simulations (sparse linear algebra class); and (3) a least-squares minimization code for interpreting magnetic force microscopy data (dense linear algebra class). Significant performance increases (between 10×× and 30×× in most cases) are seen in all three applications, demonstrating the power of GPU implementations for these types of simulations and, more generally, their associated computing classes. / Download |

12. | S.D.C. Walsh, H. Burwinkle, M.O. Saar A new partial-bounceback lattice-Boltzmann method for fluid flow through heterogeneous media, Computers and Geosciences, 35/6, pp. 1186-1193, 2009. Abstract Partial-bounceback lattice-Boltzmann methods employ a probabilistic meso-scale model that varies individual lattice node properties to reflect a material’s local permeability. These types of models have great potential in a range of geofluid, and other science and engineering, simulations of complex fluid flow. However, there are several different possible approaches for formulating partial-bounceback algorithms. This paper introduces a new partial-bounceback algorithm and compares it to two pre-existing partial-bounceback models. Unlike the two other partial-bounceback methods, the new approach conserves mass in heterogeneous media and shows improvements in simulating buoyancy-driven flow as well as diffusive processes. Further, the new model is better-suited for parallel processing implementations, resulting in faster simulations. Finally, we derive an analytical expression for calculating the permeability in all three models; a critical component for accurately matching simulation parameters to physical permeabilities. / Download |

11. | S.D.C. Walsh, M.O. Saar Magma yield stress and permeability: Insights from multiphase percolation theory, Journal of Volcanology and Geothermal Research, 177, pp. 1011-1019, 2008. Download |

10. | S.D.C. Walsh, M.O. Saar Numerical Models of Stiffness and Yield Stress Growth in Crystal-Melt Suspensions, Earth and Planetary Science Letters, 267/1-2, pp. 32-44, 2008. Abstract Magmas and other suspensions that develop sample-spanning crystal networks undergo a change in rheology from Newtonian to Bingham flow due to the onset of a finite yield stress in the crystal network. Although percolation theory provides a prediction of the crystal volume fraction at which this transition occurs, the manner in which yield stress grows with increasing crystal number densities is less-well understood. This paper discusses a simple numerical approach that models yield stress in magmatic crystalline assemblies. In this approach, the crystal network is represented by an assembly of soft-core interpenetrating cuboid (rectangular prism) particles, whose mechanical properties are simulated in a network model. The model is used to investigate the influence of particle shape and alignment anisotropy on the yield stress of crystal networks with particle volume fractions above the percolation threshold. In keeping with previous studies, the simulation predicts a local minimum in the onset of yield stress for assemblies of cubic particles, compared to those with more anisotropic shapes. The new model also predicts the growth of yield stress above (and close to) the percolation threshold. The predictions of the model are compared with results obtained from a critical path analysis. Good agreement is found between a characteristic stiffness obtained from critical path analysis, the growth in assembly stiffness predicted by the model (both of which have approximately cubic power-law exponents) and, to a lesser extent, the growth in yield stress (with a power-law exponent of 3.5). The effect of preferred particle alignment on yield stress is also investigated and found to obey similar power-law behavior. / Download |

9. | R.A. Edwards, B. Rodriguez-Brito, L. Wegley, M. Haynes, M. Breitbart, D.M. Petersen, M.O. Saar, S.C. Alexander, E.C. Alexander Jr., F. Rohwer Using pyrosequencing to shed light on deep mine microbial ecology, BMC Genomics, 2006. Abstract Background Contrasting biological, chemical and hydrogeological analyses highlights the fundamental processes that shape different environments. Generating and interpreting the biological sequence data was a costly and time-consuming process in defining an environment. Here we have used pyrosequencing, a rapid and relatively inexpensive sequencing technology, to generate environmental genome sequences from two sites in the Soudan Mine, Minnesota, USA. These sites were adjacent to each other, but differed significantly in chemistry and hydrogeology. Results Comparisons of the microbes and the subsystems identified in the two samples highlighted important differences in metabolic potential in each environment. The microbes were performing distinct biochemistry on the available substrates, and subsystems such as carbon utilization, iron acquisition mechanisms, nitrogen assimilation, and respiratory pathways separated the two communities. Although the correlation between much of the microbial metabolism occurring and the geochemical conditions from which the samples were isolated could be explained, the reason for the presence of many pathways in these environments remains to be determined. Despite being physically close, these two communities were markedly different from each other. In addition, the communities were also completely different from other microbial communities sequenced to date. Conclusion We anticipate that pyrosequencing will be widely used to sequence environmental samples because of the speed, cost, and technical advantages. Furthermore, subsystem comparisons rapidly identify the important metabolisms employed by the microbes in different environments. / Download |

8. | L.B. Christiansen, S. Hurwitz, M.O. Saar, S.E. Ingebritsen, P.A. Hsieh Seasonal seismicity at western United States volcanic centers, Earth and Planetary Science Letters, 240, pp. 307-321, 2005. no_download |

7. | M.O. Saar, M.C. Castro, C.M. Hall, M. Manga, T.P. Rose Quantifying magmatic, crustal, and atmospheric helium contributions to volcanic aquifers using all stable noble gases: Implications for magmatism and groundwater flow, Geochemistry Geophysics Geosystems, 6/3, 2005. Abstract [1] We measure all stable noble gases (He, Ne, Ar, Kr, Xe) in spring waters in the Oregon Cascades volcanic arc and in eastern Oregon, USA. We show that in order to estimate magmatic helium (He) contributions it is critical to simultaneously consider He isotopic ratios, He concentrations, and mixing of He components. Our component mixing analysis requires consideration of all measured noble gases but no other elements and is particularly insightful when strong dilution by air-saturated water has occurred. In addition, this approach can allow distinction between crustal and magmatic He components and facilitates their identification in deep groundwaters that have been diluted by near-surface water. Using this approach, we show that some cold springs on the eastern flanks of the Oregon Cascades exhibit He isotopic ratios that indicate significant magmatic He contributions comparable to those observed in thermal springs on the western flanks. Furthermore, while these magmatic He contributions are largest in deep groundwaters near the Cascades crest, greater magmatic excess He fractions than may be inferred from He isotopic ratios alone are present in all (deep) groundwaters including those at larger distances (>70 km) from the volcanic arc. We also suggest that excess He and heat discharge without dilution by air-saturated water may be restricted to spring discharge along faults. / Download |

6. | M.O. Saar, M. Manga Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modeling constraints, Journal of Geophysical Research, 109/B4, B04204, 2004. Abstract [1] We investigate the decrease in permeability, k, with depth, z, in the Oregon Cascades employing four different methods. Each method provides insight into the average permeability applicable to a different depth scale. Spring discharge models are used to infer shallow (z < 0.1 km) horizontal permeabilities. Coupled heat and groundwater flow simulations provide horizontal and vertical k for z < 1 km. Statistical investigations of the occurrences of earthquakes that are probably triggered by seasonal groundwater recharge yield vertical k for z < 5 km. Finally, considerations of magma intrusion rates and water devolatilization provide estimates of vertical k for z < 15 km. For depths >0.8 km, our results agree with the power law relationship, k = 10−14 m2 (z/1 km)−3.2, suggested by Manning and Ingebritsen [1999] for continental crust in general. However, for shallower depths (typically z ≤ 0.8 km and up to z ≤ 2) we propose an exponential relationship, k = 5 × 10−13 m2 exp (−z/0.25 km), that both fits data better (at least for the Cascades and seemingly for continental crust in general) and allows for a finite near-surface permeability and no singularity at zero depth. In addition, the suggested functions yield a smooth transition at z = 0.8 km, where their permeabilities and their gradients are similar. Permeabilities inferred from the hydroseismicity model at Mount Hood are about one order of magnitude larger than expected from the above power law. However, higher permeabilities in this region may be consistent with advective heat transfer along active faults, causing observed hot springs. Our simulations suggest groundwater recharge rates of 0.5 ≤ uR ≤ 1 m/yr and a mean background heat flow of Hb ≈ 0.080–0.134 W/m2 for the investigated region. / Download |

5. | M.A. Jellinek, M. Manga, M.O. Saar Did melting glaciers cause volcanic eruptions in eastern California? Probing the mechanics of dike formation, Journal of Geophysical Research, 109/B9, B09206, 2004. Abstract [1] A comparison of time series of basaltic and silicic eruptions in eastern California over the last 400 kyr with the contemporaneous global record of glaciation suggests that this volcanism is influenced by the growth and retreat of glaciers occurring over periods of about 40 kyr. Statistically significant cross correlations between changes in eruption frequency and the first derivative of the glacial time series imply that the temporal pattern of volcanism is influenced by the rate of change in ice volume. Moreover, calculated time lags for the effects of glacial unloading on silicic and basaltic volcanism are distinct and are 3.2 ± 4.2 kyr and 11.2 ± 2.3 kyr, respectively. A theoretical model is developed to investigate whether the increases in eruption frequency following periods of glacial unloading are a response ultimately controlled by the dynamics of dike formation. Applying results from the time series analysis leads, in turn, to estimates for the critical magma chamber overpressure required for eruption as well as constraints on the effective viscosity of the wall rocks governing dike propagation. / Download |

4. | M.O. Saar, M. Manga Seismicity induced by seasonal groundwater recharge at Mt. Hood, Oregon, Earth and Planetary Science Letters, 214, pp. 605-618, 2003. no_download |

3. | M.O. Saar, M. Manga Continuum percolation for randomly oriented soft-core prisms, Physical Review E, 65/056131, 2002. no_download |

2. | M.O. Saar, M. Manga, K.V. Cashman, S. Fremouw Numerical models of the onset of yield strength in crystal-melt suspensions, Earth and Planetary Science Letters, 187, pp. 367-379, 2001. Abstract The formation of a continuous crystal network in magmas and lavas can provide finite yield strength, τy, and can thus cause a change from Newtonian to Bingham rheology. The rheology of crystal–melt suspensions affects geological processes, such as ascent of magma through volcanic conduits, flow of lava across the Earth’s surface, melt extraction from crystal mushes under compression, convection in magmatic bodies, and shear wave propagation through partial melting zones. Here, three-dimensional numerical models are used to investigate the onset of ‘static’ yield strength in a zero-shear environment. Crystals are positioned randomly in space and can be approximated as convex polyhedra of any shape, size and orientation. We determine the critical crystal volume fraction, φc, at which a crystal network first forms. The value of φc is a function of object shape and orientation distribution, and decreases with increasing randomness in object orientation and increasing shape anisotropy. For example, while parallel-aligned convex objects yield φc=0.29, randomly oriented cubes exhibit a maximum φc of 0.22. Approximations of plagioclase crystals as randomly oriented elongated and flattened prisms (tablets) with aspect ratios between 1:4:16 and 1:1:2 yield 0.08<φc<0.20, respectively. The dependence of φc on particle orientation implies that the flow regime and resulting particle ordering may affect the onset of yield strength. φc in zero-shear environments is a lower bound for φc. Finally the average total excluded volume is used, within its limitation of being a ‘quasi-invariant’, to develop a scaling relation between τy and φ for suspensions of different particle shapes. / Download |

1. | M.O. Saar, M. Manga Permeability-porosity relationship in vesicular basalts, Geophysical Research Letters, 26/1, pp. 111-114, 1999. no_download |

### PROCEEDINGS REFEREED

8. | E. Rossi, M. Kant, F. Amann, M.O. Saar, P. Rudolf von Rohr The effects of flame-heating on rock strength: Towards a new drilling technology, Proceedings ARMA 2017, (in press). AbstractThe applicability of a combined thermo-mechanical drilling technique is investigated. The working principle of this method is based on the implementation of a heat source as a mean to either provoke thermal spallation on the surface or to weaken the rock material, when spallation is not possible. Thermal spallation drilling has already been proven to work in hard crystalline rocks, however, several difficulties hamper its application for deep resource exploitation. In order to prove the effectiveness of a combined thermo-mechanical drilling method, the forces required to export the treated sandstone material with a polycrystalline diamond compact (PDC) cutter are analyzed. The main differences between oven and flame treatments are studied by comparing the resulting strength after heat-treating the samples up to temperatures of \(650\, ^{\circ}C\) and for heating rates ranging from \(0.17 \,^{\circ}C/s\) to \(20 ^{\circ}C/s\). For moderate temperatures (\(300-450 \,^{\circ}C\)) the unconfined compressive strength after flame treatments monotonously decreased, opposed to the hardening behavior observed after oven treatments. Thermally induced intra-granular cracking and oxidation patterns served as an estimation of the treated depth due to the flame heat treatment. Therefore, conclusions on preferred operating conditions of the drilling system are drawn based on the experimental results. / no_download |

7. | D. Vogler, R.R. Settgast, C.S. Sherman, V.S. Gischig, R. Jalali, J.A. Doetsch, B. Valley, K.F. Evans, F. Amann, M.O. Saar Modeling the Hydraulic Fracture Stimulation performed for Reservoir Permeability Enhancement at the Grimsel Test Site, Switzerland, Proceedings of the 42nd Workshop on Geothermal Reservoir Engineering Stanford University, 2017. Download |

6. | N. Garapati, J. Randolph, S. Finsterle, M.O. Saar Simulating Reinjection of Produced Fluids Into the Reservoir, Proceedings of 41st Workshop on Geothermal Reservoir Engineering, 2016. no_download |

5. | N. Garapati, J.B. Randolph, J.L. Valencia Jr., M.O. Saar Design of CO2-Plume Geothermal (CPG) subsurface system for various geologic parameters, Proceedings of the Fifth International Conference on Coupled Thermo-Hydro-Mechanical-Chemical (THMC) Processes in Geosystems: Petroleum and Geothermal Reservoir Geomechanics and Energy Resource Extraction, 2015. no_download |

4. | N. Garapati, J.B. Randolph, M.O. Saar Superheating Low-Temperature Geothermal Resources to Boost Electricity Production, Proceedings of the 40th Workshop on Geothermal Reservoir Engineering 2015, 2, pp. 1210-1221, 2015. no_download |

3. | M.O. Saar, Th. Buscheck, P. Jenny, N. Garapati, J.B. Randolph, D. Karvounis, M. Chen, Y. Sun, J.M. Bielicki Numerical Study of Multi-Fluid and Multi-Level Geothermal System Performance, Proceedings World Geothermal Congress 2015, 2015. no_download |

2. | T.A. Buscheck, J.M. Bielicki, M. Chen, Y. Sun, Y. Hao, T.A. Edmunds, J.B. Randolph, M.O. Saar Multi-Fluid Sedimentary Geothermal Energy Systems for Dispatchable Renewable Electricity, Proceedings to the World Geothermal Congress, 2015. no_download |

1. | P. Bailey, J. Myre, S.C.D. Walsh, D.J. Lilja, M.O. Saar Accelerating Lattice Boltzmann Fluid Flow Simulations Using Graphics Processors, IEEE, pp. 550-557, 2009. Abstract Lattice Boltzmann Methods (LBM) are used for the computational simulation of Newtonian fluid dynamics. LBM-based simulations are readily parallelizable; they have been implemented on general-purpose processors, field-programmable gate arrays (FPGAs), and graphics processing units (GPUs). Of the three methods, the GPU implementations achieved the highest simulation performance per chip. With memory bandwidth of up to 141 GB/s and a theoretical maximum floating point performance of over 600 GFLOPS, CUDA-ready GPUs from NVIDIA provide an attractive platform for a wide range of scientific simulations, including LBM. This paper improves upon prior single-precision GPU LBM results for the D3Q19 model by increasing GPU multiprocessor occupancy, resulting in an increase in maximum performance by 20%, and by introducing a space-efficient storage method which reduces GPU RAM requirements by 50% at a slight detriment to performance. Both GPU implementations are over 28 times faster than a single-precision quad-core CPU version utilizing OpenMP. / Download |

### THESES

1. | M.O. Saar The Relationship Between Permeability, Porosity, and Microstucture in Vesicular Basalts, Theses MSc (University of Oregon), 1998. Abstract This thesis presents measurements of permeability, \(k\), porosity, \(\Phi\), and microstructural parameters of vesicular basalts. Measurements are compared with theoretical models. A percolation theory and a Kozeny-Carman model are used to interpret the measurements and to investigate relationships between porosity, microstructure, and permeability. Typical permeabilities for vesicular basalts are in the range of \(10^{-14}\) < \(k\) < \(10^{-10} m^2\). Best permeability estimates, following power laws predicted by percolation theory, are obtained when samples are used that show 'impeded aperture widening' due to rapid cooling and no bubble collapse (scoria and some flow samples). However, slowly cooled diktytaxitic samples contain elongated, 'collapsed' bubbles. Measurements indicate that the vesicle pathway network remains connected and preserves high permeabilities. Image-analysis techniques are unsuccessful if used for Kozeny-Carman equation parameter determination for vesicular materials, probably because the average interbubble aperture size that determines \(k\) is not resolved with such a technique. / no_download |

show/hide list of publications

### REFEREED PUBLICATIONS IN JOURNALS

56. | A.M.M. Leal, D.A. Kulik, W.R. Smith, M.O. Saar An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling, Pure and Applied Chemistry, 89/5, pp. 597-643, 2017. Abstract We present an overview of novel numerical methods for chemical equilibrium and kinetic calculations for complex non-ideal multiphase systems. The methods we present for equilibrium calculations are based either on Gibbs energy minimization (GEM) calculations or on solving the system of extended law of mass-action (xLMA) equations. In both methods, no a posteriori phase stability tests, and thus no tentative addition or removal of phases during or at the end of the calculations, are necessary. All potentially stable phases are considered from the beginning of the calculation, and stability indices are immediately available at the end of the computation to determine which phases are actually stable at equilibrium. Both GEM and xLMA equilibrium methods are tailored for computationally demanding applications that require many rapid local equilibrium calculations, such as reactive transport modeling. The numerical method for chemical kinetic calculations we present supports both closed and open systems, and it considers a partial equilibrium simplification for fast reactions. The method employs an implicit integration scheme that improves stability and speed when solving the often stiff differential equations in kinetic calculations. As such, it requires compositional derivatives of the reaction rates to assemble the Jacobian matrix of the resultant implicit algebraic equations that are solved at every time step. We present a detailed procedure to calculate these derivatives, and we show how the partial equilibrium assumption affects their computation. These numerical methods have been implemented in Reaktoro (reaktoro.org), an open-source software for modeling chemically reactive systems. We finish with a discussion on the comparison of these methods with others in the literature. / Download |

55. | A.J. Luhmann, B.M. Tutolo, B.C. Bagley, D.F.R. Mildner, W.E. Seyfried Jr., M.O. Saar Permeability, porosity, and mineral surface area changes in basalt cores induced by reactive transport of CO2-rich brine, Water Resources Research, 53, pp. 1-20, 2017. Abstract Four reactive flow-through laboratory experiments (two each at 0.1 mL/min and 0.01 mL/min flow rates) at 150°C and 150 bar (15 MPa) are conducted on intact basalt cores to assess changes in porosity, permeability, and surface area caused by CO2-rich fluid-rock interaction. Permeability decreases slightly during the lower flow rate experiments and increases during the higher flow rate experiments. At the higher flow rate, core permeability increases by more than one order of magnitude in one experiment and less than a factor of two in the other due to differences in preexisting flow path structure. X-ray computed tomography (XRCT) scans of pre- and post-experiment cores identify both mineral dissolution and secondary mineralization, with a net decrease in XRCT porosity of ∼0.7%–0.8% for the larger pores in all four cores. (Ultra) small-angle neutron scattering ((U)SANS) data sets indicate an increase in both (U)SANS porosity and specific surface area (SSA) over the ∼1 nm to 10 µm scale range in post-experiment basalt samples, with differences due to flow rate and reaction time. Net porosity increases from summing porosity changes from XRCT and (U)SANS analyses are consistent with core mass decreases. (U)SANS data suggest an overall preservation of the pore structure with no change in mineral surface roughness from reaction, and the pore structure is unique in comparison to previously published basalt analyses. Together, these data sets illustrate changes in physical parameters that arise due to fluid-basalt interaction in relatively low pH environments with elevated CO2 concentration, with significant implications for flow, transport, and reaction through geologic formations. / Download |

54. | A.J. Luhmann, B.M. Tutolo, C. Tan, B.M. Moskowitz, M.O. Saar, W.E. Seyfried, Jr. Whole rock basalt alteration from CO2-rich brine during flow-through experiments at 150°C and 150 bar, Chemical Geology, 453, pp. 92-110, 2017. Abstract Four flow-through experiments at 150 °C were conducted on intact cores of basalt to assess alteration and mass transfer during reaction with CO2-rich fluid. Two experiments used a flow rate of 0.1 ml/min, and two used a flow rate of 0.01 ml/min. Permeability increased for both experiments at the higher flow rate, but decreased for the lower flow rate experiments. The experimental fluid (initial pH of 3.3) became enriched in Si, Mg, and Fe upon passing through the cores, primarily from olivine and titanomagnetite dissolution and possibly pyroxene dissolution. Secondary minerals enriched in Al and Si were present on post-experimental cores, and an Fe2O3-rich phase was identified on the downstream ends of the cores from the experiments at the lower flow rate. While we could not specifically identify if siderite (FeCO3) was present in the post-experimental basalt cores, siderite was generally saturated or supersaturated in outlet fluid samples, suggesting a thermodynamic drive for Fe carbonation from basalt-H2O-CO2 reaction. Reaction path models that employ dissolution kinetics of olivine, labradorite, and enstatite also suggest siderite formation at low pH. Furthermore, fluid-rock interaction caused a relatively high mobility of the alkali metals; up to 29% and 99% of the K and Cs present in the core, respectively, were preferentially dissolved from the cores, likely due to fractional crystallization effects that made alkali metals highly accessible. Together, these datasets illustrate changes in chemical parameters that arise due to fluid-basalt interaction in relatively low pH environments with elevated CO2. / Download |

53. | J.M. Myre, E. Frahm, D.J. Lilja, M.O. Saar TNT-NN: A Fast Active Set Method for Solving Large Non-Negative Least Squares Problems, (in press). no_download |

52. | A.M.M. Leal, D. Kulik, G. Kosakowski, M.O. Saar Computational methods for reactive transport modeling: An extended law of mass-action, xLMA, method for multiphase equilibrium calculations, Advances in Water Resources, 96, pp. 405-422, 2016. Abstract We present a numerical method for multiphase chemical equilibrium calculations based on a Gibbs energy minimization approach. The method can accurately and efficiently determine the stable phase assemblage at equilibrium independently of the type of phases and species that constitute the chemical system. We have successfully applied our chemical equilibrium algorithm in reactive transport simulations to demonstrate its effective use in computationally intensive applications. We used FEniCS to solve the governing partial differential equations of mass transport in porous media using finite element methods in unstructured meshes. Our equilibrium calculations were benchmarked with GEMS3K, the numerical kernel of the geochemical package GEMS. This allowed us to compare our results with a well-established Gibbs energy minimization algorithm, as well as their performance on every mesh node, at every time step of the transport simulation. The benchmark shows that our novel chemical equilibrium algorithm is accurate, robust, and efficient for reactive transport applications, and it is an improvement over the Gibbs energy minimization algorithm used in GEMS3K. The proposed chemical equilibrium method has been implemented in Reaktoro, a unified framework for modeling chemically reactive systems, which is now used as an alternative numerical kernel of GEMS. / Download |

51. | B.M. Tutolo, D.F. Mildner, C.V. Gagnon, M.O. Saar, W.E. Seyfried Nanoscale constraints on porosity generation and fluid flow during serpentinization, Geology, 44/2, pp. 103-106, 2016. Abstract Field samples of olivine-rich rocks are nearly always serpentinized—commonly to completion—but, paradoxically, their intrinsic porosity and permeability are diminishingly low. Serpentinization reactions occur through a coupled process of fluid infiltration, volumetric expansion, and reaction-driven fracturing. Pores and reactive surface area generated during this process are the primary pathways for fluid infiltration into and reaction with serpentinizing rocks, but the size and distribution of these pores and surface area have not yet been described. Here, we utilize neutron scattering techniques to present the first measurements of the evolution of pore size and specific surface area distribution in partially serpentinized rocks. Samples were obtained from the ca. 2 Ma Atlantis Massif oceanic core complex located off-axis of the Mid-Atlantic Ridge and an olivine-rich outcrop of the ca. 1.1 Ga Duluth Complex of the North American Mid-Continent Rift. Our measurements and analyses demonstrate that serpentine and accessory phases form with their own, inherent porosity, which accommodates the bulk of diffusive fluid flow during serpentinization and thereby permits continued serpentinization after voluminous serpentine minerals fill reaction-generated porosity. / Download |

50. | A.M.M. Leal, D.A. Kulik, M.O. Saar Enabling Gibbs energy minimization algorithms to use equilibrium constants of reactions in multiphase equilibrium calculations, Chemical Geology, 437, pp. 170-181, 2016. Abstract The geochemical literature provides numerous thermodynamic databases compiled from equilibrium constants of reactions. These databases are typically used in speciation calculations based on the law of mass action (LMA) approach. Unfortunately, such LMA databases cannot be directly used in equilibrium speciation methods based on the Gibbs energy minimization (GEM) approach because of their lack of standard chemical potentials of species. Therefore, we present in this work a simple conversion approach that calculates apparent standard chemical potentials of species from equilibrium constants of reactions. We assess the consistency and accuracy of the use of apparent standard chemical potentials in GEM algorithms by benchmarking equilibrium speciation calculations using GEM and LMA methods with the same LMA database. In all cases, we use PHREEQC to perform the LMA calculations, and we use its LMA databases to calculate the equilibrium constants of reactions. GEM calculations are performed using a Gibbs energy minimization method of Reaktoro — a unified open-source framework for numerical modeling of chemically reactive systems. By comparing the GEM and LMA results, we show that the use of apparent standard chemical potentials in GEM methods produces consistent and accurate equilibrium speciation results, thus validating our new, practical conversion technique that enables GEM algorithms to take advantage of many existing LMA databases, consequently extending and diversifying their range of applicability. / Download |

49. | T.A. Buscheck, J.M. Bielicki, T.A. Edmunds, Y. Hao, Y. Sun, J.B. Randolph, M.O. Saar Multifluid geo-energy systems: Using geologic CO2 storage for geothermal energy production and grid-scale energy storage in sedimentary basins, Geosphere, 12/3, pp. 1-19, 2016. Download |

48. | B.M. Tutolo, X.-Z. Kong, W.E. Seyfried Jr., M.O. Saar High performance reactive transport simulations examining the effects of thermal, hydraulic, and chemical (THC) gradients on fluid injectivity at carbonate CCUS reservoir scales, International Journal of Greenhouse Gas Control, 39, pp. 285-301, 2015. Abstract Carbonate minerals and CO2 are both considerably more soluble at low temperatures than they are at elevated temperatures. This inverse solubility has led a number of researchers to hypothesize that injecting low-temperature (i.e., less than the background reservoir temperature) CO2 into deep, saline reservoirs for CO2 Capture, Utilization, and Storage (CCUS) will dissolve CO2 and carbonate minerals near the injection well and subsequently exsolve and re-precipitate these phases as the fluids flow into the geothermally warm portion of the reservoir. In this study, we utilize high performance computing to examine the coupled effects of cool CO2 injection and background hydraulic head gradients on reservoir-scale mineral volume changes. We employ the fully coupled reactive transport simulator PFLOTRAN with calculations distributed over up to 800 processors to test 21 scenarios designed to represent a range of reservoir depths, hydraulic head gradients, and CO2 injection rates and temperatures. In the default simulations, 50 °C CO2 is injected at a rate of 50 kg/s into a 200 bar, 100 °C calcite or dolomite reservoir. By comparing these simulations with others run at varying conditions, we show that the effect of cool CO2 injection on reservoir-scale mineral volume changes tends to be relatively minor. We conclude that the low heat capacity of CO2 effectively prevents low-temperature CO2 injection from decreasing the temperature across large portions of the simulated carbonate reservoirs. This small thermal perturbation, combined with the low relative permeability of brine within the supercritical CO2 plume, yields limited dissolution and precipitation effects directly attributable to cool CO2 injection. Finally, we calculate that relatively high water-to-rock ratios, which may occur over much longer CCUS reservoir lifetimes or in materials with sufficiently high brine relative permeability within the supercritical CO2 plume, would be required to substantially affect injectivity through thermally-induced mineral dissolution and precipitation. Importantly, this study shows the utility of reservoir scale-reactive transport simulators for testing hypotheses and placing laboratory-scale observations into a CCUS reservoir-scale context. / Download |

47. | B.M. Tutolo, A.J. Luhmann, X.-Z. Kong, M.O. Saar, W.E. Seyfried Jr. CO2 sequestration in feldspar-rich sandstone: Coupled evolution of fluid chemistry, mineral reaction rates, and hydrogeochemical properties, Geochimica Et Cosmochimica Acta, 160, pp. 132-154, 2015. Abstract To investigate CO / Download_{2} Capture, Utilization, and Storage (CCUS) in sandstones, we performed three 150 °C flow-through experiments on K-feldspar-rich cores from the Eau Claire formation. By characterizing fluid and solid samples from these experiments using a suite of analytical techniques, we explored the coupled evolution of fluid chemistry, mineral reaction rates, and hydrogeochemical properties during CO_{2} sequestration in feldspar-rich sandstone. Overall, our results confirm predictions that the heightened acidity resulting from supercritical CO_{2} injection into feldspar-rich sandstone will dissolve primary feldspars and precipitate secondary aluminum minerals. A core through which CO_{2}-rich deionized water was recycled for 52 days decreased in bulk permeability, exhibited generally low porosity associated with high surface area in post-experiment core sub-samples, and produced an Al hydroxide secondary mineral, such as boehmite. However, two samples subjected to ?3 day single-pass experiments run with CO_{2}-rich, 0.94 mol/kg NaCl brines decreased in bulk permeability, showed generally elevated porosity associated with elevated surface area in post-experiment core sub-samples, and produced a phase with kaolinite-like stoichiometry. CO_{2}-induced metal mobilization during the experiments was relatively minor and likely related to Ca mineral dissolution. Based on the relatively rapid approach to equilibrium, the relatively slow near-equilibrium reaction rates, and the minor magnitudes of permeability changes in these experiments, we conclude that CCUS systems with projected lifetimes of several decades are geochemically feasible in the feldspar-rich sandstone end-member examined here. Additionally, the observation that K-feldspar dissolution rates calculated from our whole-rock experiments are in good agreement with literature parameterizations suggests that the latter can be utilized to model CCUS in K-feldspar-rich sandstone. Finally, by performing a number of reactive transport modeling experiments to explore processes occurring during the flow-through experiments, we have found that the overall progress of feldspar hydrolysis is negligibly affected by quartz dissolution, but significantly impacted by the rates of secondary mineral precipitation and their effect on feldspar saturation state. The observations produced here are critical to the development of models of CCUS operations, yet more work, particularly in the quantification of coupled dissolution and precipitation processes, will be required in order to produce models that can accurately predict the behavior of these systems. |

46. | N. Garapati, J.B. Randolph, M.O. Saar Brine displacement by CO2, energy extraction rates, and lifespan of a CO2-limited CO2-Plume Geothermal (CPG) system with a horizontal production well, Geothermics, 55, pp. 182-194, 2015. Abstract Several studies suggest that CO2-based geothermal energy systems may be operated economically when added to ongoing geologic CO2 sequestration. Alternatively, we demonstrate here that CO2-Plume Geothermal (CPG) systems may be operated long-term with a finite amount of CO2. We analyze the performance of such CO2-limited CPG systems as a function of various geologic and operational parameters. We find that the amount of CO2 required increases with reservoir depth, permeability, and well spacing and decreases with larger geothermal gradients. Furthermore, the onset of reservoir heat depletion decreases for increasing geothermal gradients and for both particularly shallow and deep reservoirs. / Download |

45. | B.M. Tutolo, A.T. Schaen, M.O. Saar, W.E. Seyfried Jr. Implications of the redissociation phenomenon for mineral-buffered fluids and aqueous species transport at elevated temperatures and pressures, Applied Geochemistry, 55, pp. 119-127, 2015. Abstract Aqueous species equilibrium constants and activity models form the foundation of the complex speciation codes used to model the geochemistry of geothermal energy production, extremophilic ecosystems, ore deposition, and a variety of other processes. Researchers have shown that a simple three species model (i.e., Na+, Cl?, and NaCl(aq)) can accurately describe conductivity measurements of concentrated NaCl and KCl solutions at elevated temperatures and pressures (Sharygin et al., 2002). In this model, activity coefficients of the charged species (e.g., Na+, K+, Cl?) become sufficiently low that the complexes must redisocciate with increasing salt concentration in order to meet equilibrium constant constraints. Redissociation decreases the proportion of the elements bound up as neutral complexes, and thereby increases the true ionic strength of the solution. In this contribution, we explore the consequences of the redissociation phenomenon in albite–paragonite–quartz (APQ) buffered systems. We focus on the implications of the redissociation phenomenon for mineral solubilities, particularly the observation that, at certain temperatures and pressures, calculated activities of charged ions in solution remain practically constant even as element concentrations increase from <1 molal to 4.5 molal. Finally, we note that redissociation has a similar effect on pH, and therefore aqueous speciation, in APQ-hosted systems. The calculations and discussion presented here are not limited to APQ-hosted systems, but additionally apply to many others in which the dominant cations and anions can form neutral complexes. / Download |

44. | B.M. Adams, T.H. Kuehn, J.M. Bielicki, J.B. Randolph, M.O. Saar A comparison of electric power output of CO2 Plume Geothermal (CPG) and brine geothermal systems for varying reservoir conditions, Applied Energy, 140, pp. 365-377, 2015. Abstract In contrast to conventional hydrothermal systems or enhanced geothermal systems, CO2 Plume Geothermal (CPG) systems generate electricity by using CO2 that has been geothermally heated due to sequestration in a sedimentary basin. Four CPG and two brine-based geothermal systems are modeled to estimate their power production for sedimentary basin reservoir depths between 1 and 5km, geothermal temperature gradients from 20 to 50°Ckm-1, reservoir permeabilities from 1×10-15 to 1×10-12m2 and well casing inner diameters from 0.14m to 0.41m. Results show that CPG direct-type systems produce more electricity than brine-based geothermal systems at depths between 2 and 3km, and at permeabilities between 10-14 and 10-13m2, often by a factor of two. This better performance of CPG is due to the low kinematic viscosity of CO2, relative to brine at those depths, and the strong thermosiphon effect generated by CO2. When CO2 is used instead of R245fa as the secondary working fluid in an organic Rankine cycle (ORC), the power production of both the CPG and the brine-reservoir system increases substantially; for example, by 22% and 20% for subsurface brine and CO2 systems, respectively, with a 35°Ckm-1 thermal gradient, 0.27m production and 0.41m injection well diameters, and 5×10-14m2 reservoir permeability. / Download |

43. | A.J. Luhmann, M. Covington, J. Myre, M. Perne, S.W. Jones, C.E. Alexander Jr., M.O. Saar Thermal damping and retardation in karst conduits, Hydrology and Earth System Sciences, 19/1, pp. 137-157, 2015. Abstract Water temperature is a non-conservative tracer in the environment. Variations in recharge temperature are damped and retarded as water moves through an aquifer due to heat exchange between water and rock. However, within karst aquifers, seasonal and short-term fluctuations in recharge temperature are often transmitted over long distances before they are fully damped. Using analytical solutions and numerical simulations, we develop relationships that describe the effect of flow path properties, flow-through time, recharge characteristics, and water and rock physical properties on the damping and retardation of thermal peaks/troughs in karst conduits. Using these relationships, one can estimate the thermal retardation and damping that would occur under given conditions with a given conduit geometry. Ultimately, these relationships can be used with thermal damping and retardation field data to estimate parameters such as conduit diameter. We also examine sets of numerical simulations where we relax some of the assumptions used to develop these relationships, testing the effects of variable diameter, variable velocity, open channels, and recharge shape on thermal damping and retardation to provide some constraints on uncertainty. Finally, we discuss a multitracer experiment that provides some field confirmation of our relationships. High temporal resolution water temperature data are required to obtain sufficient constraints on the magnitude and timing of thermal peaks and troughs in order to take full advantage of water temperature as a tracer. / Download |

42. | B.M. Tutolo, A.J. Luhmann, X.-Z. Kong, M.O. Saar, W.E. Seyfried Jr. Experimental observation of permeability changes in dolomite at CO2 sequestration conditions, Environmental Science and Technology, 48/4, pp. 2445-2452, 2014. Abstract Injection of cool CO2 into geothermally warm carbonate reservoirs for storage or geothermal energy production may lower near-well temperature and lead to mass transfer along flow paths leading away from the well. To investigate this process, a dolomite core was subjected to a 650 h, high pressure, CO2 saturated, flow-through experiment. Permeability increased from 10–15.9 to 10–15.2 m2 over the initial 216 h at 21 °C, decreased to 10–16.2 m2 over 289 h at 50 °C, largely due to thermally driven CO2 exsolution, and reached a final value of 10–16.4 m2 after 145 h at 100 °C due to continued exsolution and the onset of dolomite precipitation. Theoretical calculations show that CO2 exsolution results in a maximum pore space CO2 saturation of 0.5, and steady state relative permeabilities of CO2 and water on the order of 0.0065 and 0.1, respectively. Post-experiment imagery reveals matrix dissolution at low temperatures, and subsequent filling-in of flow passages at elevated temperature. Geochemical calculations indicate that reservoir fluids subjected to a thermal gradient may exsolve and precipitate up to 200 cm3 CO2 and 1.5 cm3 dolomite per kg of water, respectively, resulting in substantial porosity and permeability redistribution. / Download |

41. | A.J. Luhmann, X.-Z. Kong, B.M. Tutolo, N. Garapati, B.C. Bagley, M.O. Saar, W.E. Seyfried Jr. Experimental dissolution of dolomite by CO2-charged brine at 100oC and 150 bar: Evolution of porosity, permeability, and reactive surface area, Chemical Geology, 380, pp. 145-160, 2014. Abstract Hydrothermal flow experiments of single-pass injection of CO2-charged brine were conducted on nine dolomite cores to examine fluid–rock reactions in dolomite reservoirs under geologic carbon sequestration conditions. Post-experimental X-ray computed tomography (XRCT) analysis illustrates a range of dissolution patterns, and significant increases in core bulk permeability were measured as the dolomite dissolved. Outflow fluids were below dolomite saturation, and cation concentrations decreased with time due to reductions in reactive surface area with reaction progress. To determine changes in reactive surface area, we employ a power-law relationship between reactive surface area and porosity (Luquot and Gouze, 2009). The exponent in this relationship is interpreted to be a geometrical parameter that controls the degree of surface area change per change in core porosity. Combined with XRCT reconstructions of dissolution patterns, we demonstrate that this exponent is inversely related to both the flow path diameter and tortuosity of the dissolution channel. Even though XRCT reconstructions illustrate dissolution at selected regions within each core, relatively high Ba and Mn recoveries in fluid samples suggest that dissolution occurred along the core’s entire length and width. Analysis of porosity–permeability data indicates an increase in the rate of permeability enhancement per increase in porosity with reaction progress as dissolution channels lengthen along the core. Finally, we incorporate the surface area–porosity model of Luquot and Gouze (2009) with our experimentally fit parameters into TOUGHREACT to simulate experimental observations. / Download |

40. | B.M. Adams, T.H. Kuehn, J.M. Bielicki, J.B. Randolph, M.O. Saar On the importance of the thermosiphon effect in CPG (CO2-Plume geothermal) power systems, Energy, 69, pp. 409-418, 2014. Abstract CPG (CO2 Plume Geothermal) energy systems use CO2 to extract thermal energy from naturally permeable geologic formations at depth. CO2 has advantages over brine: high mobility, low solubility of amorphous silica, and higher density sensitivity to temperature. The density of CO2 changes substantially between geothermal reservoir and surface plant, resulting in a buoyancy-driven convective current – a thermosiphon – that reduces or eliminates pumping requirements. We estimated and compared the strength of this thermosiphon for CO2 and for 20 weight percent NaCl brine for reservoir depths up to 5 km and geothermal gradients of 20, 35, and 50 °C/km. We found that through the reservoir, CO2 has a pressure drop approximately 3–12 times less than brine at the same mass flowrate, making the CO2 thermosiphon sufficient to produce power using reservoirs as shallow as 0.5 km. At 2.5 km depth with a 35 °C/km gradient – the approximate western U.S. continental mean – the CO2 thermosiphon converted approximately 10% of the energy extracted from the reservoir to fluid circulation, compared to less than 1% with brine, where additional mechanical pumping is necessary. We found CO2 is a particularly advantageous working fluid at depths between 0.5 km and 3 km. / Download |

39. | B.M. Tutolo, X.-Z. Kong, W.E. Seyfried Jr., M.O. Saar Internal consistency in aqueous geochemical data revisited: Applications to the aluminum system, Geochimica et Cosmochimica Acta, 133, pp. 216-234, 2014. Abstract Internal consistency of thermodynamic data has long been considered vital for confident calculations of aqueous geochemical processes. However, an internally consistent mineral thermodynamic data set is not necessarily consistent with calculations of aqueous species thermodynamic properties due, potentially, to improper or inconsistent constraints used in the derivation process. In this study, we attempt to accommodate the need for a mineral thermodynamic data set that is internally consistent with respect to aqueous species thermodynamic properties by adapting the least squares optimization methods of Powell and Holland (1985). This adapted method allows for both the derivation of mineral thermodynamic properties from fluid chemistry measurements of solutions in equilibrium with mineral assemblages, as well as estimates of the uncertainty on the derived results. Using a large number of phase equilibria, solubility, and calorimetric measurements, we have developed a thermodynamic data set of 12 key aluminum-bearing mineral phases. These data are derived to be consistent with Na+ and K+ speciation data presented by Shock and Helgeson (1988), H4SiO4(aq) data presented by Stefánsson (2001), and the Al speciation data set presented by Tagirov and Schott (2001). Many of the constraining phase equilibrium measurements are exactly the same as those used to develop other thermodynamic data, yet our derived values tend to be quite different than some of the others’ due to our choices of reference data. The differing values of mineral thermodynamic properties have implications for calculations of Al mineral solubilities; specifically, kaolinite solubilities calculated with the developed data set are as much as 6.75 times lower and 73% greater than those calculated with Helgeson et al. (1978) and Holland and Powell (2011) data, respectively. Where possible, calculations and experimental data are compared at low T, and the disagreement between the two sources reiterates the common assertion that low-T measurements of phase equilibria and mineral solubilities in the aluminum system rarely represent equilibrium between water and well-crystallized, aluminum-bearing minerals. As an ancillary benefit of the derived data, we show that it may be combined with high precision measurements of aqueous complex association constants to derive neutral species activity coefficients in supercritical fluids. Although this contribution is specific to the aluminum system, the methods and concepts developed here can help to improve the calculation of water–rock interactions in a broad range of earth systems. / Download |

38. | N. Garapati, J.B. Randolph, J.L. Valencia Jr., M.O. Saar CO2-Plume Geothermal (CPG) Heat Extraction in Multi-layered Geologic Reservoirs, Energy Procedia, 63, pp. 7631-7643, 2014. Abstract CO2-Plume Geothermal (CPG) technology involves injecting CO2 into natural, highly permeable geologic units to extract energy. The subsurface CO2 absorbs heat from the reservoir, buoyantly rises to the surface, and drives a power generation system. The CO2 is then cooled and reinjected underground. Here, we analyze the effects of multi-layered geologic reservoirs on CPG system performance by examining the CO2 mass fraction in the produced fluid, pore-fluid pressure buildup during operation, and heat energy extraction rates. The produced CO2 mass fraction depends on the stratigraphic positions of highly permeable layers which also affect the pore-fluid pressure drop across the reservoir. / Download |

37. | T.A. Buscheck, J.M. Bielicki, M. Chen, Y. Sun, Y. Hao, T.A. Edmunds, M.O. Saar, J.B. Randolph Integrating CO2 Storage with Geothermal Resources for Dispatchable Renewable Electricity, Energy procedia, 63, pp. 7619-7630, 2014. Abstract We present an approach that uses the huge fluid and thermal storage capacity of the subsurface, together with geologic CO2 storage, to harvest, store, and dispatch energy from subsurface (geothermal) and surface (solar, nuclear, fossil) thermal resources, as well as energy from electrical grids. Captured CO2 is injected into saline aquifers to store pressure, generate artesian flow of brine, and provide an additional working fluid for efficient heat extraction and power conversion. Concentric rings of injection and production wells are used to create a hydraulic divide to store pressure, CO2, and thermal energy. Such storage can take excess power from the grid and excess/waste thermal energy, and dispatch that energy when it is demanded, enabling increased penetration of variable renewables. Stored CO2 functions as a cushion gas to provide enormous pressure-storage capacity and displaces large quantities of brine, which can be desalinated and/or treated for a variety of beneficial uses. Geothermal power and energy-storage applications may generate enough revenues to justify CO2 capture costs. / Download |

36. | B.M. Adams, T.H. Kuehn, J.B. Randolph, M.O. Saar The reduced pumping power requirements from increasing the injection well fluid density, Transactions – Geothermal Resources Council, 37, pp. 667-672, 2013. no_download |

35. | J.B. Randolph, M.O. Saar, J.M. Bielicki Geothermal energy production at geologic CO2 sequestration sites: Impact of thermal drawdown on reservoir pressure, Energy Procedia, 37, pp. 6625-6635, 2013. Abstract Recent geotechnical research shows that geothermal heat can be efficiently mined by circulating carbon dioxide through naturally permeable rock formations — a method called CO2 Plume Geothermal — the same geologic reservoirs that are suitable for deep saline aquifer CO2 sequestration or enhanced oil recovery. This paper describes the effect of thermal drawdown on reservoir pressure buildup during sequestration operations, revealing that geothermal heat mining can decrease overpressurization by 10% or more. / Download |

34. | R. Gottardi, P.-H. Kao, M.O. Saar, Ch. Teyssier Effects of permeability fields on fluid, heat, and oxygen isotope transport in extensional detachment systems, Geochemistry, Geophysics, Geosystems, 14/5, pp. 1493-1522, 2013. Abstract [1] Field studies of Cordilleran metamorphic core complexes indicate that meteoric fluids permeated the upper crust down to the detachment shear zone and interacted with highly deformed and recrystallized (mylonitic) rocks. The presence of fluids in the brittle/ductile transition zone is recorded in the oxygen and hydrogen stable isotope compositions of the mylonites and may play an important role in the thermomechanical evolution of the detachment shear zone. Geochemical data show that fluid flow in the brittle upper crust is primarily controlled by the large-scale fault-zone architecture. We conduct continuum-scale (i.e., large-scale, partial bounce-back) lattice-Boltzman fluid, heat, and oxygen isotope transport simulations of an idealized cross section of a metamorphic core complex. The simulations investigate the effects of crust and fault permeability fields as well as buoyancy-driven flow on two-way coupled fluid and heat transfer and resultant exchange of oxygen isotopes between meteoric fluid and rock. Results show that fluid migration to middle to lower crustal levels is fault controlled and depends primarily on the permeability contrast between the fault zone and the crustal rocks. High fault/crust permeability ratios lead to channelized flow in the fault and shear zones, while lower ratios allow leakage of the fluids from the fault into the crust. Buoyancy affects mainly flow patterns (more upward directed) and, to a lesser extent, temperature distributions (disturbance of the geothermal field by ~25°C). Channelized fluid flow in the shear zone leads to strong vertical and horizontal thermal gradients, comparable to field observations. The oxygen isotope results show δ18O depletion concentrated along the fault and shear zones, similar to field data. / Download |

33. | S.D.C. Walsh, M.O. Saar Developing extensible lattice-Boltzmann simulators for general-purpose graphics-processing units, Communications in Computational Physics, 13/3, pp. 867-879, 2013. Abstract Lattice-Boltzmann methods are versatile numerical modeling techniques capable of reproducing a wide variety of fluid-mechanical behavior. These methods are well suited to parallel implementation, particularly on the single-instruction multiple data (SIMD) parallel processing environments found in computer graphics processing units (GPUs). Although recent programming tools dramatically improve the ease with which GPUbased applications can be written, the programming environment still lacks the flexibility available to more traditional CPU programs. In particular, it may be difficult to develop modular and extensible programs that require variable on-device functionality with current GPU architectures. This paper describes a process of automatic code generation that overcomes these difficulties for lattice-Boltzmann simulations. It details the development of GPU-based modules for an extensible lattice-Boltzmann simulation package – LBHydra. The performance of the automatically generated code is compared to equivalent purposewritten codes for both single-phase,multiphase, andmulticomponent flows. The flexibility of the new method is demonstrated by simulating a rising, dissolving droplet moving through a porous medium with user generated lattice-Boltzmann models and subroutines. / Download |

32. | X.-Z. Kong, M.O. Saar DBCreate: A SUPCRT92-based program for producing EQ3/6, TOUGHREACT, and GWB thermodynamic databases at user-defined T and P, Computers and Geosciences, 51, pp. 415-417, 2013. Abstract SUPCRT92 is a widely used software package for calculating the standard thermodynamic properties of minerals, gases, aqueous species, and reactions. However, it is labor-intensive and error-prone to use it directly to produce databases for geochemical modeling programs such as EQ3/6, the Geochemist’s Workbench, and TOUGHREACT. DBCreate is a SUPCRT92-based software program written in FORTRAN90/95 and was developed in order to produce the required databases for these programs in a rapid and convenient way. This paper describes the overall structure of the program and provides detailed usage instructions. / Download |

31. | X.-Z. Kong, M.O. Saar Numerical study of the effects of permeability heterogeneity on density-driven convective mixing during CO2 dissolution storage, Int. J. Greenhouse Gas Control, 19, pp. 160-173, 2013. Abstract Permanence and security of carbon dioxide (CO2) in geologic formations requires dissolution of CO2 into brine, which slightly increases the brine density. Previous studies have shown that this small increase in brine density induces convective currents, which greatly enhances the mixing efficiency and thus CO2 storage capacity and rate in the brine. Density-driven convection, in turn, is known to be largely dominated by permeability heterogeneity. This study explores the relationship between the process of density-driven convection and the permeability heterogeneity of an aquifer during CO2 dissolution storage, using high-resolution numerical simulations. While the porosity is kept constant, the heterogeneity of the aquifer is introduced through a spatially varying permeability field, characterized by the Dykstra-Parsons coefficient and the correlation length. Depending on the concentration profile of dissolved CO2, we classify the convective finger patterns as dispersive, preferential, and unbiased fingering. Our results indicate that the transition between unbiased and both preferential and dispersive fingering is mainly governed by the Dykstra-Parsons coefficient, whereas the transition between preferential and dispersive fingering is controlled by the permeability correlation length. Furthermore, we find that the CO2 dissolution flux at the top boundary will reach a time-independent steady state. Although this flux strongly correlates with permeability distribution, it generally increases with the permeability heterogeneity when the correlation length is less than the system size. / Download |

30. | A.J. Luhmann, X.-Z. Kong, B.M. Tutolo, K. Ding, M.O. Saar, W.E. Seyfried Jr. Permeability reduction produced by grain reorganization and accumulation of exsolved CO2 during geologic carbon sequestration: A new CO2 trapping mechanism, Environmental Science and Technlogy, 47/1, pp. 242-251, 2013. Abstract Carbon sequestration experiments were conducted on uncemented sediment and lithified rock from the Eau Claire Formation, which consisted primarily of K-feldspar and quartz. Cores were heated to accentuate reactivity between fluid and mineral grains and to force CO2 exsolution. Measured permeability of one sediment core ultimately reduced by 4 orders of magnitude as it was incrementally heated from 21 to 150 °C. Water-rock interaction produced some alteration, yielding sub-?m clay precipitation on K-feldspar grains in the core’s upstream end. Experimental results also revealed abundant newly formed pore space in regions of the core, and in some cases pores that were several times larger than the average grain size of the sediment. These large pores likely formed from elevated localized pressure caused by rapid CO2 exsolution within the core and/or an accumulating CO2 phase capable of pushing out surrounding sediment. CO2 filled the pores and blocked flow pathways. Comparison with a similar experiment using a solid arkose core indicates that CO2 accumulation and grain reorganization mainly contributed to permeability reduction during the heated sediment core experiment. This suggests that CO2 injection into sediments may store more CO2 and cause additional permeability reduction than is possible in lithified rock due to grain reorganization. / Download |

29. | M. Covington, A.F. Banwell, J. Gulley, M.O. Saar Quantifying the effects of glacier conduit geometry and recharge on proglacial hydrograph form, Journal of Hydrology, 414-415, pp. 59-71, 2012. Abstract The configuration of glacier hydrological systems is often inferred from proxy data, such as hydrographs, that are collected in proglacial streams. Seasonal changes in the peakedness of hydrographs are thought to reflect changes in the configuration of the subglacial drainage system. However, the amount of information that proglacial hydrographs contain about drainage system configurations depends critically on the degree to which the drainage systems modify recharge hydrographs. If the drainage system does not modify recharge hydrographs, then proglacial hydrographs primarily reflect the recharge conditions produced by supraglacial inputs. Here, we develop a theoretical framework to determine the circumstances under which glacier drainage systems can modify recharge hydrographs and the circumstances under which recharge pulses pass through glaciers unchanged. We address the capability of single conduits, simple arborescent conduit networks, and linked cavity systems to modify diurnal recharge pulses. Simulations of discharge through large sets of such systems demonstrate that, unless large reservoirs or significant constrictions are present, the discharge hydrographs of simple glacial conduit systems are nearly identical to their recharge hydrographs. Conduit systems tend not to modify hydrographs because the changes in storage within englacial and subglacial conduit networks on short time scales are typically small compared to their ability to transmit water. This finding suggests that proglacial hydrographs reflect a variety of factors, including surface melt rate, surface water transfer, and subglacial water transfer. In many cases the influence of suglacial processes may be relatively minor. As a result, the evolution of proglacial hydrographs cannot be used unambiguously to infer changes in the structure or efficiency of englacial or subglacial hydrological systems, without accurate knowledge of the nature of the recharge hydrograph driving the flow. / Download |

28. | J.B. Randolph, B.M. Adams, T.H. Kuehn, M.O. Saar Wellbore heat transfer in CO2-based geothermal systems, Geothermal Resources Council (GRC) Transactions, 36, pp. 549-554, 2012. no_download |

27. | M.D. Covington, A.J. Luhmann, C. Wicks, M.O. Saar Process length scales and longitudinal damping in karst conduits, Journal of Geophysical Research – Earth Surface, 117, F01025, 2012. Abstract [1] Simple mathematical models often allow an intuitive grasp of the function of physical systems. We develop a mathematical framework to investigate reactive or dissipative transport processes within karst conduits. Specifically, we note that for processes that occur within a characteristic timescale, advection along the conduit produces a characteristic process length scale. We calculate characteristic length scales for the propagation of thermal and electrical conductivity signals along karst conduits. These process lengths provide a quantitative connection between karst conduit geometry and the signals observed at a karst spring. We show that water input from the porous/fractured matrix is also characterized by a length scale and derive an approximation that accounts for the influence of matrix flow on the transmission of signals through the aquifer. The single conduit model is then extended to account for conduits with changing geometries and conduit flow networks, demonstrating how these concepts can be applied in more realistic conduit geometries. We introduce a recharge density function, ϕR, which determines the capability of an aquifer to damp a given signal, and cast previous explanations of spring variability within this framework. Process lengths are a general feature of karst conduits and surface streams, and we conclude with a discussion of other potential applications of this conceptual and mathematical framework. / Download |

26. | S.C. Alexander, M.O. Saar Improved characterization of small u for Jacob pumping test analysis methods, Ground Water, 50/2, pp. 256-265, 2012. Abstract Numerous refinements have been proposed to traditional pumping test analyses, yet many hydrogeologists continue to use the Jacob method due to its simplicity. Recent research favors hydraulic tomography and inverse numerical modeling of pumping test data. However, at sites with few wells, or relatively short screens, the data requirements of these methods may be impractical within physical and fiscal constraints. Alternatively, an improved understanding of the assumptions and limitations of Theis and, due to their widespread usage, Jacob analyses, leads to improved interpretations in data-poor environments. A fundamental requirement of Jacob is a “small” value of u = f(r2/t), with radial distance, r, and pumping time, t. However, selection of a too stringent (i.e., too low) maximum permissible u-value, umax, results in rejection of usable data from wells beyond a maximum radius, rmax. Conversely, data from small radii, less than rmin, where turbulent- and vertical-flow components arise, can result in acceptance of inappropriate data. Usage of drawdown data from wells too close to the pumping well, and exclusion of data from wells deemed too far, can cause unrealistic aquifer transmissivity, permeability, and storativity determinations. Here, data from an extensive well field in a glacial-outwash aquifer in north-central Minnesota, USA, are used to develop a new estimate for umax. Traditionally quoted values for umax range from 0.01 to 0.05. Our proposed value for Jacob distance-drawdown analyses is significantly higher with umax up to 0.2, resulting in larger allowable rmax-values and a higher likelihood of inclusion of additional wells in such pumping test analyses. / Download |

25. | M.D. Covington, A.J. Luhmann, F. Gabrovsek, M.O. Saar, I. Willis, C.M. Wicks Mechanisms of heat exchange between water and rock in karst conduits, Water Resources Research, 47, W10514/10, 2011. Abstract [1] Previous studies, motivated by understanding water quality, have explored the mechanisms for heat transport and heat exchange in surface streams. In karst aquifers, temperature signals play an additional important role since they carry information about internal aquifer structures. Models for heat transport in karst conduits have previously been developed; however, these models make different, sometimes contradictory, assumptions. Additionally, previous models of heat transport in karst conduits have not been validated using field data from conduits with known geometries. Here we use analytical solutions of heat transfer to examine the relative importance of heat exchange mechanisms and the validity of the assumptions made by previous models. The relative importance of convection, conduction, and radiation is a function of time. Using a characteristic timescale, we show that models neglecting rock conduction produce spurious results in realistic cases. In contrast to the behavior of surface streams, where conduction is often negligible, conduction through the rock surrounding a conduit determines heat flux at timescales of weeks and longer. In open channel conduits, radiative heat flux can be significant. In contrast, convective heat exchange through the conduit air is often negligible. Using the rules derived from our analytical analysis, we develop a numerical model for heat transport in a karst conduit. Our model compares favorably to thermal responses observed in two different karst settings: a cave stream fed via autogenic recharge during a snowmelt event, and an allogenically recharged cave stream that experiences continuous temperature fluctuations on many timescales. / Download |

24. | J.B. Randolph, M.O. Saar Impact of reservoir permeability on the choice of subsurface geothermal heat exchange fluid: CO2 versus water and native brine, Geothermal Resources Council (GRC) Transactions, 35, pp. 521-526, 2011. no_download |

23. | J.B. Randolph, M.O. Saar Combining geothermal energy capture with geologic carbon dioxide sequestration, Geophysical Research Letters, 38, L10401, 2011. Abstract [1] Geothermal energy offers clean, renewable, reliable electric power with no need for grid-scale energy storage, yet its use has been constrained to the few locations worldwide with naturally high geothermal heat resources and groundwater availability. We present a novel approach with the potential to permit expansion of geothermal energy utilization: heat extraction from naturally porous, permeable formations with CO2 as the injected subsurface working fluid. Fluid-mechanical simulations reveal that the significantly higher mobility of CO2, compared to water, at the temperature/pressure conditions of interest makes CO2 an attractive heat exchange fluid. We show numerically that, compared to conventional water-based and engineered geothermal systems, the proposed approach provides up to factors of 2.9 and 5.0, respectively, higher geothermal heat energy extraction rates. Consequently, more regions worldwide could be economically used for geothermal electricity production. Furthermore, as the injected CO2 is eventually geologically sequestered, such power plants would have negative carbon footprints. / Download |

22. | M.A. Davis, S.D.C. Walsh, M.O. Saar Statistically reconstructing continuous isotropic and anisotropic two-phase media while preserving macroscopic material properties, Physical Review E, 83, 026706, 2011. Abstract We propose a method to generate statistically similar reconstructions of two-phase media. As with previous work, we initially characterize the microstructure of the material using two-point correlation functions (a subset of spatial correlation functions) and then generate numerical reconstructions using a simulated annealing method that preserves the geometric relationships of the material’s phase of interest. However, in contrast to earlier contributions that consider reconstructions composed of discrete arrays of pixels or voxels alone, we generate reconstructions based on assemblies of continuous, three-dimensional, interpenetrating objects. The result is a continuum description of the material microstructure (as opposed to a discretized or pixelated description), capable of efficiently representing large disparities in scale. Different reconstruction methods are considered based on distinct combinations of two-point correlation functions of varying degrees of complexity. The quality of the reconstruction methods are evaluated by comparing the total pore fraction, specific surface area of the percolating cluster, pore fraction of the percolating cluster, tortuosity, and permeability of the reconstructions to those of a set of reference assemblies. Elsewhere it has been proposed that two-phase media could be statistically reproduced with only two spatial correlation functions: the two-point probability function (the probability that two points lie within the same phase) and the lineal path function (the probability that a line between two points lies entirely within the same phase). We find that methods employing the two-point probability function and lineal path function are improved if the percolating cluster volume is also considered in the reconstruction. However, to reproduce more complicated geometric assemblies, we find it necessary to employ the two-point probability, two-point cluster, and lineal path function in addition to the percolating cluster volume to produce a generally accurate statistical reconstruction. / Download |

21. | J.B. Randolph, M.O. Saar Coupling carbon dioxide sequestration with geothermal energy capture in naturally permeable, porous geologic formations: Implications for CO2 sequestration, Energy Procedia, 4, pp. 2206-2213, 2011. Abstract Carbon dioxide (CO2) sequestration in deep saline aquifers and exhausted oil and natural gas fields has been widely considered as a means for reducing CO2 emissions to the atmosphere as a counter-measure to global warming. However, rather than treating CO2 merely as a waste fluid in need of permanent disposal, we propose that it could also be used as a working fluid in geothermal energy capture, as its thermodynamic and fluid mechanical properties suggest it transfers geothermal heat more efficiently than water. Energy production and sales in conjunction with sequestration would improve the economic viability of CO2 sequestration, a critical challenge for large-scale implementation of the technology. In addition, using CO2 as the working fluid in geothermal power systems may permit utilization of lower temperature geologic formations than those that are currently deemed economically viable, leading to more widespread utilization of geothermal energy. Here, we present the results of early-stage calculations demonstrating the geothermal energy capture potential of CO2-based geothermal systems and implications of such energy capture for the economic viability of geologic CO2 sequestration. / Download |

20. | M.O. Saar Review: Geothermal heat as a tracer of large-scale groundwater flow and as a means to determine permeability fields, special theme issue on Environmental Tracers and Groundwater Flow, editor-invited peer-reviewed contribution, Hydrogeology Journal, 19, pp. 31-52, 2011. Abstract A review of coupled groundwater and heat transfer theory is followed by an introduction to geothermal measurement techniques. Thereafter, temperature-depth profiles (geotherms) and heat discharge at springs to infer hydraulic parameters and processes are discussed. Several studies included in this review state that minimum permeabilities of approximately 5 × 10−17 < kmin <10−15 m2 are required to observe advective heat transfer and resultant geotherm perturbations. Permeabilities below kmin tend to cause heat-conduction-dominated systems, precluding inversion of temperature fields for groundwater flow patterns and constraint of permeabilities other than being / Download |

19. | J. Myre, S.D.C. Walsh, D.J. Lilja, M.O. Saar Performance analysis of single-phase multiphase, and multicomponent lattice-Boltzmann fluid flow simulations on GPU clusters, Concurrency and Computation: Practice and Experience, 23, pp. 332-350, 2010. Abstract Abstract The lattice-Boltzmann method is well suited for implementation in single-instruction multiple-data (SIMD) environments provided by general purpose graphics processing units (GPGPUs). This paper discusses the integration of these GPGPU programs with OpenMP to create lattice-Boltzmann applications for multi-GPU clusters. In addition to the standard single-phase single-component lattice-Boltzmann method, the performances of more complex multiphase, multicomponent models are also examined. The contributions of various GPU lattice-Boltzmann parameters to the performance are examined and quantified with a statistical model of the performance using Analysis of Variance (ANOVA). By examining single- and multi-GPU lattice-Boltzmann simulations with ANOVA, we show that all the lattice-Boltzmann simulations primarily depend on effects corresponding to simulation geometry and decomposition, and not on the architectural aspects of GPU. Additionally, using ANOVA we confirm that the metrics of Efficiency and Utilization are not suitable for memory-bandwidth-dependent codes. Copyright © 2010 John Wiley & Sons, Ltd. / Download |

18. | S. Dasgupta, M.O. Saar, R.L. Edwards, C.-C. Shen, H. Cheng, C.E. Alexander Jr. Three thousand years of extreme rainfall events recorded in stalagmites from Spring Valley Caverns, Minnesota, Earth and Planetary Science Letters, 300, pp. 46-54, 2010. Abstract Annual layer analysis in two stalagmites collected from Spring Valley Caverns, southeastern Minnesota, reveals hydrological response of the cave to extreme rainfall events in the Midwest, USA. Cave-flooding events are identified within the two samples by the presence of detrital layers composed of clay sized particles. Comparison with instrumental records of precipitation demonstrates a strong correlation between these cave-flood events and extreme rainfall observed in the Upper Mississippi Valley. A simple model is developed to assess the nature of rainfall capable of flooding the cave. The model is first calibrated to the last 50-yr (1950–1998 A.D.) instrumental record of daily precipitation data for the town of Spring Valley and verified with the first 50 yr of record from 1900 to 1949 A.D. Frequency analysis shows that these extreme flood events have increased from the last half of the nineteenth century. Comparison with other paleohydrological records shows increased occurrence of extreme rain events during periods of higher moisture availability. Our study implies that increased moisture availability in the Midwestern region, due to rise in temperature from global warming could lead to an increase in the occurrence of extreme rainfall events. / Download |

17. | S.D.C. Walsh, M.O. Saar Interpolated lattice-Boltzmann boundary conditions for surface reaction kinetics, Physical Review E, 82, 066703, 2010. Abstract This paper describes a method for implementing surface reaction kinetics in lattice Boltzmann simulations. The interpolated boundary conditions are capable of simulating surface reactions and dissolution at both stationary and moving solid-fluid and fluid-fluid interfaces. Results obtained with the boundary conditions are compared to analytical solutions for first-order and constant-flux kinetic surface reactions in a one-dimensional half space, as well as to the analytical solution for evaporation from the surface of a cylinder. Excellent agreement between analytical and simulated results is obtained for a wide range of diffusivities, lattice velocities, and surface reaction rates. The boundary model’s ability to represent dissolution in binary fluid mixtures is demonstrated by modeling diffusion from a rising bubble and dissolution of a droplet near a flat plate. / Download |

16. | J.B. Randolph, M.O. Saar Coupling geothermal energy capture with carbon dioxide sequestration in naturally permeable, porous geologic formations: A comparison with enhanced geothermal systems, Geothermal Resources Council (GRC) Transactions, 34, pp. 433-438, 2010. no_download |

15. | S.D.C. Walsh, M.O. Saar Macroscale lattice-Boltzmann methods for low-Peclet-number solute and heat transport in heterogeneous porous media, Water Resources Research, 46, W07517, 2010. Abstract [1] This paper introduces new methods for simulating subsurface solute and heat transport in heterogeneous media using large-scale lattice-Boltzmann models capable of representing both macroscopically averaged porous media and open channel flows. Previous examples of macroscopically averaged lattice-Boltzmann models for solute and heat transport are only applicable to homogeneous media. Here, we extend these models to properly account for heterogeneous pore-space distributions. For simplicity, in the majority of this paper we assume low Peclet number flows with an isotropic dispersion tensor. Nevertheless, this approach may also be extended to include anisotropic-dispersion by using multiple relaxation time lattice-Boltzmann methods. We describe two methods for introducing heterogeneity into macroscopically averaged lattice-Boltzmann models. The first model delivers the desired behavior by introducing an additional time-derivative term to the collision rule; the second model by separately weighting symmetric and anti-symmetric components of the fluid packet densities. Chapman-Enskog expansions are conducted on the governing equations of the two models, demonstrating that the correct constitutive behavior is obtained in both cases. In addition, methods for improving model stability at low porosities are also discussed: (1) an implicit formulation of the model; and (2) a local transformation that normalizes the lattice-Boltzmann model by the local porosity. The model performances are evaluated through comparisons of simulated results with analytical solutions for one- and two-dimensional flows, and by comparing model predictions to finite element simulations of advection isotropic-dispersion in heterogeneous porous media. We conclude by presenting an example application, demonstrating the ability of the new models to couple with simulations of reactive flow and changing flow geometry: a simulation of groundwater flow through a carbonate system. / Download |

14. | M. Covington, C.M. Wicks, M.O. Saar A dimensionless number describing the effects of recharge and geometry on discharge from simple karst aquifers, Water Resources Research, 45, W11410, 2009. Abstract [1] The responses of karstic aquifers to storms are often used to obtain information about aquifer geometry. In general, spring hydrographs are a function of both system geometry and recharge. However, the majority of prior work on storm pulses through karst has not studied the effect of recharge on spring hydrographs. To examine the relative importance of geometry and recharge, we break karstic aquifers into elements according to the manner of their response to transient flow and demonstrate that each element has a characteristic response timescale. These fundamental elements are full pipes, open channels, reservoir/constrictions, and the porous matrix. Taking the ratio of the element timescale with the recharge timescale produces a dimensionless number, γ, that is used to characterize aquifer response to a storm event. Using sets of simulations run with randomly selected element parameters, we demonstrate that each element type has a critical value of γ below which the shape of the spring hydrograph is dominated by the shape of the recharge hydrograph and above which the spring hydrograph is significantly modified by the system geometry. This allows separation of particular element/storm pairs into recharge-dominated and geometry-dominated regimes. While most real karstic aquifers are complex combinations of these elements, we draw examples from several karst systems that can be represented by single elements. These examples demonstrate that for real karstic aquifers full pipe and open channel elements are generally in the recharge-dominated regime, whereas reservoir/constriction elements can fall in either the recharge- or geometry-dominated regimes. / Download |

13. | S.D.C. Walsh, M.O. Saar, P. Bailey, D.J. Lilja Accelerating geoscience and engineering system simulations on graphics hardware, Computers and Geosciences, 35/12, pp. 2353-2364, 2009. Abstract Many complex natural systems studied in the geosciences are characterized by simple local-scale interactions that result in complex emergent behavior. Simulations of these systems, often implemented in parallel using standard central processing unit (CPU) clusters, may be better suited to parallel processing environments with large numbers of simple processors. Such an environment is found in graphics processing units (GPUs) on graphics cards. This paper discusses GPU implementations of three example applications from computational fluid dynamics, seismic wave propagation, and rock magnetism. These candidate applications involve important numerical modeling techniques, widely employed in physical system simulations, that are themselves examples of distinct computing classes identified as fundamental to scientific and engineering computing. The presented numerical methods (and respective computing classes they belong to) are: (1) a lattice-Boltzmann code for geofluid dynamics (structured grid class); (2) a spectral-finite-element code for seismic wave propagation simulations (sparse linear algebra class); and (3) a least-squares minimization code for interpreting magnetic force microscopy data (dense linear algebra class). Significant performance increases (between 10×× and 30×× in most cases) are seen in all three applications, demonstrating the power of GPU implementations for these types of simulations and, more generally, their associated computing classes. / Download |

12. | S.D.C. Walsh, H. Burwinkle, M.O. Saar A new partial-bounceback lattice-Boltzmann method for fluid flow through heterogeneous media, Computers and Geosciences, 35/6, pp. 1186-1193, 2009. Abstract Partial-bounceback lattice-Boltzmann methods employ a probabilistic meso-scale model that varies individual lattice node properties to reflect a material’s local permeability. These types of models have great potential in a range of geofluid, and other science and engineering, simulations of complex fluid flow. However, there are several different possible approaches for formulating partial-bounceback algorithms. This paper introduces a new partial-bounceback algorithm and compares it to two pre-existing partial-bounceback models. Unlike the two other partial-bounceback methods, the new approach conserves mass in heterogeneous media and shows improvements in simulating buoyancy-driven flow as well as diffusive processes. Further, the new model is better-suited for parallel processing implementations, resulting in faster simulations. Finally, we derive an analytical expression for calculating the permeability in all three models; a critical component for accurately matching simulation parameters to physical permeabilities. / Download |

11. | S.D.C. Walsh, M.O. Saar Magma yield stress and permeability: Insights from multiphase percolation theory, Journal of Volcanology and Geothermal Research, 177, pp. 1011-1019, 2008. Download |

10. | S.D.C. Walsh, M.O. Saar Numerical Models of Stiffness and Yield Stress Growth in Crystal-Melt Suspensions, Earth and Planetary Science Letters, 267/1-2, pp. 32-44, 2008. Abstract Magmas and other suspensions that develop sample-spanning crystal networks undergo a change in rheology from Newtonian to Bingham flow due to the onset of a finite yield stress in the crystal network. Although percolation theory provides a prediction of the crystal volume fraction at which this transition occurs, the manner in which yield stress grows with increasing crystal number densities is less-well understood. This paper discusses a simple numerical approach that models yield stress in magmatic crystalline assemblies. In this approach, the crystal network is represented by an assembly of soft-core interpenetrating cuboid (rectangular prism) particles, whose mechanical properties are simulated in a network model. The model is used to investigate the influence of particle shape and alignment anisotropy on the yield stress of crystal networks with particle volume fractions above the percolation threshold. In keeping with previous studies, the simulation predicts a local minimum in the onset of yield stress for assemblies of cubic particles, compared to those with more anisotropic shapes. The new model also predicts the growth of yield stress above (and close to) the percolation threshold. The predictions of the model are compared with results obtained from a critical path analysis. Good agreement is found between a characteristic stiffness obtained from critical path analysis, the growth in assembly stiffness predicted by the model (both of which have approximately cubic power-law exponents) and, to a lesser extent, the growth in yield stress (with a power-law exponent of 3.5). The effect of preferred particle alignment on yield stress is also investigated and found to obey similar power-law behavior. / Download |

9. | R.A. Edwards, B. Rodriguez-Brito, L. Wegley, M. Haynes, M. Breitbart, D.M. Petersen, M.O. Saar, S.C. Alexander, E.C. Alexander Jr., F. Rohwer Using pyrosequencing to shed light on deep mine microbial ecology, BMC Genomics, 2006. Abstract Background Contrasting biological, chemical and hydrogeological analyses highlights the fundamental processes that shape different environments. Generating and interpreting the biological sequence data was a costly and time-consuming process in defining an environment. Here we have used pyrosequencing, a rapid and relatively inexpensive sequencing technology, to generate environmental genome sequences from two sites in the Soudan Mine, Minnesota, USA. These sites were adjacent to each other, but differed significantly in chemistry and hydrogeology. Results Comparisons of the microbes and the subsystems identified in the two samples highlighted important differences in metabolic potential in each environment. The microbes were performing distinct biochemistry on the available substrates, and subsystems such as carbon utilization, iron acquisition mechanisms, nitrogen assimilation, and respiratory pathways separated the two communities. Although the correlation between much of the microbial metabolism occurring and the geochemical conditions from which the samples were isolated could be explained, the reason for the presence of many pathways in these environments remains to be determined. Despite being physically close, these two communities were markedly different from each other. In addition, the communities were also completely different from other microbial communities sequenced to date. Conclusion We anticipate that pyrosequencing will be widely used to sequence environmental samples because of the speed, cost, and technical advantages. Furthermore, subsystem comparisons rapidly identify the important metabolisms employed by the microbes in different environments. / Download |

8. | L.B. Christiansen, S. Hurwitz, M.O. Saar, S.E. Ingebritsen, P.A. Hsieh Seasonal seismicity at western United States volcanic centers, Earth and Planetary Science Letters, 240, pp. 307-321, 2005. no_download |

7. | M.O. Saar, M.C. Castro, C.M. Hall, M. Manga, T.P. Rose Quantifying magmatic, crustal, and atmospheric helium contributions to volcanic aquifers using all stable noble gases: Implications for magmatism and groundwater flow, Geochemistry Geophysics Geosystems, 6/3, 2005. Abstract [1] We measure all stable noble gases (He, Ne, Ar, Kr, Xe) in spring waters in the Oregon Cascades volcanic arc and in eastern Oregon, USA. We show that in order to estimate magmatic helium (He) contributions it is critical to simultaneously consider He isotopic ratios, He concentrations, and mixing of He components. Our component mixing analysis requires consideration of all measured noble gases but no other elements and is particularly insightful when strong dilution by air-saturated water has occurred. In addition, this approach can allow distinction between crustal and magmatic He components and facilitates their identification in deep groundwaters that have been diluted by near-surface water. Using this approach, we show that some cold springs on the eastern flanks of the Oregon Cascades exhibit He isotopic ratios that indicate significant magmatic He contributions comparable to those observed in thermal springs on the western flanks. Furthermore, while these magmatic He contributions are largest in deep groundwaters near the Cascades crest, greater magmatic excess He fractions than may be inferred from He isotopic ratios alone are present in all (deep) groundwaters including those at larger distances (>70 km) from the volcanic arc. We also suggest that excess He and heat discharge without dilution by air-saturated water may be restricted to spring discharge along faults. / Download |

6. | M.O. Saar, M. Manga Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modeling constraints, Journal of Geophysical Research, 109/B4, B04204, 2004. Abstract [1] We investigate the decrease in permeability, k, with depth, z, in the Oregon Cascades employing four different methods. Each method provides insight into the average permeability applicable to a different depth scale. Spring discharge models are used to infer shallow (z < 0.1 km) horizontal permeabilities. Coupled heat and groundwater flow simulations provide horizontal and vertical k for z < 1 km. Statistical investigations of the occurrences of earthquakes that are probably triggered by seasonal groundwater recharge yield vertical k for z < 5 km. Finally, considerations of magma intrusion rates and water devolatilization provide estimates of vertical k for z < 15 km. For depths >0.8 km, our results agree with the power law relationship, k = 10−14 m2 (z/1 km)−3.2, suggested by Manning and Ingebritsen [1999] for continental crust in general. However, for shallower depths (typically z ≤ 0.8 km and up to z ≤ 2) we propose an exponential relationship, k = 5 × 10−13 m2 exp (−z/0.25 km), that both fits data better (at least for the Cascades and seemingly for continental crust in general) and allows for a finite near-surface permeability and no singularity at zero depth. In addition, the suggested functions yield a smooth transition at z = 0.8 km, where their permeabilities and their gradients are similar. Permeabilities inferred from the hydroseismicity model at Mount Hood are about one order of magnitude larger than expected from the above power law. However, higher permeabilities in this region may be consistent with advective heat transfer along active faults, causing observed hot springs. Our simulations suggest groundwater recharge rates of 0.5 ≤ uR ≤ 1 m/yr and a mean background heat flow of Hb ≈ 0.080–0.134 W/m2 for the investigated region. / Download |

5. | M.A. Jellinek, M. Manga, M.O. Saar Did melting glaciers cause volcanic eruptions in eastern California? Probing the mechanics of dike formation, Journal of Geophysical Research, 109/B9, B09206, 2004. Abstract [1] A comparison of time series of basaltic and silicic eruptions in eastern California over the last 400 kyr with the contemporaneous global record of glaciation suggests that this volcanism is influenced by the growth and retreat of glaciers occurring over periods of about 40 kyr. Statistically significant cross correlations between changes in eruption frequency and the first derivative of the glacial time series imply that the temporal pattern of volcanism is influenced by the rate of change in ice volume. Moreover, calculated time lags for the effects of glacial unloading on silicic and basaltic volcanism are distinct and are 3.2 ± 4.2 kyr and 11.2 ± 2.3 kyr, respectively. A theoretical model is developed to investigate whether the increases in eruption frequency following periods of glacial unloading are a response ultimately controlled by the dynamics of dike formation. Applying results from the time series analysis leads, in turn, to estimates for the critical magma chamber overpressure required for eruption as well as constraints on the effective viscosity of the wall rocks governing dike propagation. / Download |

4. | M.O. Saar, M. Manga Seismicity induced by seasonal groundwater recharge at Mt. Hood, Oregon, Earth and Planetary Science Letters, 214, pp. 605-618, 2003. no_download |

3. | M.O. Saar, M. Manga Continuum percolation for randomly oriented soft-core prisms, Physical Review E, 65/056131, 2002. no_download |

2. | M.O. Saar, M. Manga, K.V. Cashman, S. Fremouw Numerical models of the onset of yield strength in crystal-melt suspensions, Earth and Planetary Science Letters, 187, pp. 367-379, 2001. Abstract The formation of a continuous crystal network in magmas and lavas can provide finite yield strength, τy, and can thus cause a change from Newtonian to Bingham rheology. The rheology of crystal–melt suspensions affects geological processes, such as ascent of magma through volcanic conduits, flow of lava across the Earth’s surface, melt extraction from crystal mushes under compression, convection in magmatic bodies, and shear wave propagation through partial melting zones. Here, three-dimensional numerical models are used to investigate the onset of ‘static’ yield strength in a zero-shear environment. Crystals are positioned randomly in space and can be approximated as convex polyhedra of any shape, size and orientation. We determine the critical crystal volume fraction, φc, at which a crystal network first forms. The value of φc is a function of object shape and orientation distribution, and decreases with increasing randomness in object orientation and increasing shape anisotropy. For example, while parallel-aligned convex objects yield φc=0.29, randomly oriented cubes exhibit a maximum φc of 0.22. Approximations of plagioclase crystals as randomly oriented elongated and flattened prisms (tablets) with aspect ratios between 1:4:16 and 1:1:2 yield 0.08<φc<0.20, respectively. The dependence of φc on particle orientation implies that the flow regime and resulting particle ordering may affect the onset of yield strength. φc in zero-shear environments is a lower bound for φc. Finally the average total excluded volume is used, within its limitation of being a ‘quasi-invariant’, to develop a scaling relation between τy and φ for suspensions of different particle shapes. / Download |

1. | M.O. Saar, M. Manga Permeability-porosity relationship in vesicular basalts, Geophysical Research Letters, 26/1, pp. 111-114, 1999. no_download |

### PROCEEDINGS REFEREED

8. | E. Rossi, M. Kant, F. Amann, M.O. Saar, P. Rudolf von Rohr The effects of flame-heating on rock strength: Towards a new drilling technology, Proceedings ARMA 2017, (in press). AbstractThe applicability of a combined thermo-mechanical drilling technique is investigated. The working principle of this method is based on the implementation of a heat source as a mean to either provoke thermal spallation on the surface or to weaken the rock material, when spallation is not possible. Thermal spallation drilling has already been proven to work in hard crystalline rocks, however, several difficulties hamper its application for deep resource exploitation. In order to prove the effectiveness of a combined thermo-mechanical drilling method, the forces required to export the treated sandstone material with a polycrystalline diamond compact (PDC) cutter are analyzed. The main differences between oven and flame treatments are studied by comparing the resulting strength after heat-treating the samples up to temperatures of \(650\, ^{\circ}C\) and for heating rates ranging from \(0.17 \,^{\circ}C/s\) to \(20 ^{\circ}C/s\). For moderate temperatures (\(300-450 \,^{\circ}C\)) the unconfined compressive strength after flame treatments monotonously decreased, opposed to the hardening behavior observed after oven treatments. Thermally induced intra-granular cracking and oxidation patterns served as an estimation of the treated depth due to the flame heat treatment. Therefore, conclusions on preferred operating conditions of the drilling system are drawn based on the experimental results. / no_download |

7. | D. Vogler, R.R. Settgast, C.S. Sherman, V.S. Gischig, R. Jalali, J.A. Doetsch, B. Valley, K.F. Evans, F. Amann, M.O. Saar Modeling the Hydraulic Fracture Stimulation performed for Reservoir Permeability Enhancement at the Grimsel Test Site, Switzerland, Proceedings of the 42nd Workshop on Geothermal Reservoir Engineering Stanford University, 2017. Download |

6. | N. Garapati, J. Randolph, S. Finsterle, M.O. Saar Simulating Reinjection of Produced Fluids Into the Reservoir, Proceedings of 41st Workshop on Geothermal Reservoir Engineering, 2016. no_download |

5. | N. Garapati, J.B. Randolph, J.L. Valencia Jr., M.O. Saar Design of CO2-Plume Geothermal (CPG) subsurface system for various geologic parameters, Proceedings of the Fifth International Conference on Coupled Thermo-Hydro-Mechanical-Chemical (THMC) Processes in Geosystems: Petroleum and Geothermal Reservoir Geomechanics and Energy Resource Extraction, 2015. no_download |

4. | N. Garapati, J.B. Randolph, M.O. Saar Superheating Low-Temperature Geothermal Resources to Boost Electricity Production, Proceedings of the 40th Workshop on Geothermal Reservoir Engineering 2015, 2, pp. 1210-1221, 2015. no_download |

3. | M.O. Saar, Th. Buscheck, P. Jenny, N. Garapati, J.B. Randolph, D. Karvounis, M. Chen, Y. Sun, J.M. Bielicki Numerical Study of Multi-Fluid and Multi-Level Geothermal System Performance, Proceedings World Geothermal Congress 2015, 2015. no_download |

2. | T.A. Buscheck, J.M. Bielicki, M. Chen, Y. Sun, Y. Hao, T.A. Edmunds, J.B. Randolph, M.O. Saar Multi-Fluid Sedimentary Geothermal Energy Systems for Dispatchable Renewable Electricity, Proceedings to the World Geothermal Congress, 2015. no_download |

1. | P. Bailey, J. Myre, S.C.D. Walsh, D.J. Lilja, M.O. Saar Accelerating Lattice Boltzmann Fluid Flow Simulations Using Graphics Processors, IEEE, pp. 550-557, 2009. Abstract Lattice Boltzmann Methods (LBM) are used for the computational simulation of Newtonian fluid dynamics. LBM-based simulations are readily parallelizable; they have been implemented on general-purpose processors, field-programmable gate arrays (FPGAs), and graphics processing units (GPUs). Of the three methods, the GPU implementations achieved the highest simulation performance per chip. With memory bandwidth of up to 141 GB/s and a theoretical maximum floating point performance of over 600 GFLOPS, CUDA-ready GPUs from NVIDIA provide an attractive platform for a wide range of scientific simulations, including LBM. This paper improves upon prior single-precision GPU LBM results for the D3Q19 model by increasing GPU multiprocessor occupancy, resulting in an increase in maximum performance by 20%, and by introducing a space-efficient storage method which reduces GPU RAM requirements by 50% at a slight detriment to performance. Both GPU implementations are over 28 times faster than a single-precision quad-core CPU version utilizing OpenMP. / Download |

### THESES

1. | M.O. Saar The Relationship Between Permeability, Porosity, and Microstucture in Vesicular Basalts, Theses MSc (University of Oregon), 1998. Abstract This thesis presents measurements of permeability, \(k\), porosity, \(\Phi\), and microstructural parameters of vesicular basalts. Measurements are compared with theoretical models. A percolation theory and a Kozeny-Carman model are used to interpret the measurements and to investigate relationships between porosity, microstructure, and permeability. Typical permeabilities for vesicular basalts are in the range of \(10^{-14}\) < \(k\) < \(10^{-10} m^2\). Best permeability estimates, following power laws predicted by percolation theory, are obtained when samples are used that show 'impeded aperture widening' due to rapid cooling and no bubble collapse (scoria and some flow samples). However, slowly cooled diktytaxitic samples contain elongated, 'collapsed' bubbles. Measurements indicate that the vesicle pathway network remains connected and preserves high permeabilities. Image-analysis techniques are unsuccessful if used for Kozeny-Carman equation parameter determination for vesicular materials, probably because the average interbubble aperture size that determines \(k\) is not resolved with such a technique. / no_download |