Abstract

It has long been recognized that quartz precipitation from circulating hydrothermal fluids may reduce porosity and permeability near intrusions. However, the magnitude of permeability changes and potential feedbacks between flow, heat transfer, and quartz precipitation/dissolution remain largely unquantified. Here, we present numerical simulations of fluid convection around upper crustal intrusions which explicitly incorporate the feedback between quartz solubility and rock permeability. As groundwater is heated to ~350°C, silica dissolves from the host rock, increasing porosity and permeability. Further heating to supercritical conditions leads to intensive quartz precipitation and consequent permeability reduction. The initial host rock permeability and porosity are found to be main controls on the magnitude and timescales of permeability changes. While the permeability changes induced by quartz precipitation are moderate in host rocks with a primary porosity ≥ 0.05, quartz precipitation may reduce rock permeability by more than an order of magnitude in host rocks with a primary porosity of 0.025. Zones of quartz precipitation transiently change locations as the intrusion cools, thereby limiting the clogging effect, except for host rocks with low initial porosity. This permeability reduction occurs in timescales of hundreds of years in host rocks with initial high permeability and thousands of years in host rocks with intermediate permeability.

1. Introduction

Rock permeability in hydrothermal systems is affected by an interplay between fluid flow, rock alteration, and the stress state. The dissolution of primary host rock minerals and precipitation of secondary minerals by aqueous fluids change the hydraulic and mechanical properties of the rock media [1, 2]. Depending on the spatial configuration and temporal evolution of zones of precipitation and dissolution, fluid-rock interaction may have contrasting effects on porosity and permeability. Whereas extensive secondary mineral precipitation reduces permeability and thereby causes more diffuse fluid flow, dissolution increases permeability and may aid in focusing fluid flow into high permeability zones due to positive feedbacks between flow and reaction [35].

Several authors have proposed that quartz precipitation close to an intrusive heat source reduces permeability over timescales of decades to hundreds of years and thereby restricts fluid flow at high temperatures [612]. Silica is typically the most abundant dissolved solid in dilute geothermal fluids of meteoric origin, with concentrations closely approaching equilibrium with respect to quartz at temperatures above ca. 180°C [13, 14]. The possible paths of precipitation and/or dissolution can be understood from Figure 1. The concentration of silica in pure water assuming quartz equilibrium increases with temperature and reaches a maximum solubility between 350 and 450°C (depending on fluid pressure), above which it strongly decreases. This region of retrograde quartz solubility overlaps with the temperature range that develops during hydrothermal convection around intrusions [15, 16]. In addition, quartz solubility is strongly reduced during boiling [1719].

As liquid is heated to supercritical conditions, quartz deposition and consequent permeability reduction may decrease fluid mass fluxes and reduce the rate of heat transfer from the intrusion. The quantitative effects of these potential feedback mechanisms have not yet been studied, although rock permeability exerts a strong control on the spatial extent and enthalpy of supercritical zones near intrusions [20] and the style of the hydrothermal system [16]. Below a permeability of ~10−16 m2, fluid advection is limited and heat transfer is dominated by thermal conduction through the rock [15, 21, 22]. While host rocks with permeability near 10−15 m2 may potentially develop supercritical zones extending hundreds of meters above the intrusion, supercritical fluid flow is confined to a thin layer (10 s of m) enveloping the intrusion in host rocks with a higher permeability of 10−14 m2 [20]. Rock permeability controls supercritical geothermal resource formation by determining the temperature to which circulating fluids are heated during passage near an intrusion and the extent to which supercritical fluids mix with cooler circulating waters upon ascent [20]. This points to a complex interaction between quartz deposition and permeability changes accompanying supercritical resource formation. Although a larger share of the fluid circulating near an intrusion emplaced in intermediate permeability host rocks is heated to the high temperatures of quartz retrograde solubility, the much higher fluid fluxes and steeper temperature gradients in high permeability systems suggest a greater overall potential for porosity and permeability reduction.

The magnitude of permeability reduction induced by mineral deposition strongly depends on the initial or “primary” porosity of the host rock [10, 23, 24]. A given volume of secondary mineral precipitates will reduce the permeability of initially low porosity rocks more than in high porosity rocks. While the clogging of pore space by secondary mineral precipitates in low porosity rocks limits the maximum extent of alteration, alteration in high porosity rocks may proceed to completion, with a total replacement of the host rock by secondary minerals [24]. The porosity of rocks hosting geothermal reservoirs can vary over a wide range, from <0.05 to >0.6 [25, 26]. However, evidence from deep boreholes in active geothermal systems indicates that porosity decreases with depth, and so-called “tight rocks” with low porosity (<0.03) are common at depths ≥ 2 km [2, 27, 28]. On top of these variations, the relation between permeability and porosity is complex, depending on the interconnectedness of the pore space, the relation between fracture and matrix porosity, and the stress state of the rock [2931].

In this study, we investigate whether and under what conditions quartz precipitation affects permeability near shallow intrusions. We present numerical models simulating quartz precipitation and silica dissolution during fluid convection around transiently cooling intrusions and notably do not consider the effect of significant fluid expulsion from hydrous magmas [32]. We explore the role of initial host rock permeability, primary porosity, and intrusion depth on the dynamics of silica transport and permeability changes induced by the dissolution of silica from the host rock and precipitation of secondary quartz. We find that quartz precipitation does not strongly reduce permeability unless the initial host rock porosity is low (~0.025). Depending on the initial permeability and porosity of the host rocks, quartz precipitation has contrasting and unexpected effects on supercritical zones, reducing fluxes of supercritical vapor and rates of heat transport near the intrusion if initial host rock permeability is 10−15 m2, but promoting the development of larger supercritical resources if the initial permeability of the host rocks is high (10−14 m2).

2. Methods

To study the mutual feedbacks between fluid flow, quartz precipitation/dissolution, and rock permeability, a numerical model incorporating three main components was implemented: (i) a numerical method for multiphase flow of compressible, variably miscible H2O; (ii) a formulation for quartz solubility and silica transport; and (iii) a dynamic permeability model describing the effect of quartz precipitation and dissolution on pore space.

2.1. Multiphase Fluid Flow

The governing equations of multiphase mass and energy conservation are solved using a continuum porous media approach with a pressure-enthalpy-based formulation in a control volume finite element method numerical scheme using the complex systems modeling platform (CSMP++), which has been described in detail by Weis et al. [33] and is thus only described briefly. Phase velocities () were obtained using an extended two-phase (i = liquid, vapor) form of Darcy’s law: where is the rock permeability, the relative permeability of the fluid phase, dynamic viscosity, the pressure gradient, phase density, and the gravitational acceleration vector. A linear relative permeability model with a liquid residual saturation of 0.3 and vapor residual saturation of zero is adopted [15, 33]. The pore fluid consists of pure water, and all fluid properties correspond to the formulation of Haar et al. [34].

The conservation of fluid mass is given by where is rock porosity, is the volumetric saturation of each phase, and is a fluid source term. Energy conservation accounts for the conduction of heat in the rock and advection of enthalpy by fluid, given by where is the thermal conductivity of the rock and is an energy source term. The energy source term results from the assumption that fluid and rock are in local thermal equilibrium, and total enthalpy is distributed over the fluid and rock contained in a control volume such that they are at the same temperature [33].

2.2. Silica Solubility, Precipitation, and Dissolution

The solubility of quartz in aqueous fluids has been measured over a wide range of temperatures, pressures, and fluid compositions (e.g., [3543]). At temperatures above approximately ~180°C, the concentration of dissolved silica in natural waters closely approaches equilibrium with respect to quartz according to the reaction: where refers to solvated aqueous silica independent of hydration state [37, 40]. In this study, we adopt the formulation of Manning [40], in which silica solubility is calculated as a function of temperature and fluid density. Where both liquid and vapor phases are present, silica solubility is calculated separately for each phase, and the bulk fluid solubility reflects a weighted mass balance. Depending on the bulk fluid enthalpy, the solubility of quartz in boiling fluids can vary strongly at a given temperature (Figure 1(b)).

Silica transport occurs via fluid advection according to

The mass of quartz contained in the rock within a given control volume, , is adjusted at each time step such that the mass fraction of dissolved silica in each fluid phase, , corresponds to equilibrium with respect to quartz. The mass of quartz precipitated or dissolved in each time step is given by the difference between the mass concentration of silica in the fluid after advection of liquid and/or vapor () and the bulk fluid equilibrium concentration of quartz (): with mass concentrations specified in kg SiO2 per m3 fluid. The volume-normalized rate of quartz precipitation or dissolution for the fluid-rock medium is calculated as a postprocessing step, according to where is the dimensions of mol SiO2 per unit volume [m−3] of rock-fluid medium and unit time [s−1], is the length of the time step in seconds, and is the molar mass of silica (0.06008 kg mol−1). For a given time step, if the bulk fluid concentration of silica exceeds the equilibrium concentration, is positive and quartz precipitates; if bulk fluid silica concentration is less than the equilibrium concentration, is negative and silica dissolves from the host rock. The rock is assumed to be always quartz-bearing; that is, silica can freely dissolve from the rock even if an equivalent mass of secondary quartz has not yet been precipitated. This assumption is valid for quartz-bearing intrusive host rocks, such as granite, andesite, or rhyolite or for basaltic host rocks that already contain abundant secondary quartz because of previous episodes of heating and fluid circulation. However, the results may not be fully representative for systems with rocks that lack free quartz.

The assumption of quartz equilibrium may not be valid at all conditions in natural systems. For example, quartz solubility below ~180°C is controlled by an amorphous silica polymorph such as chalcedony [6]. Moreover, the slow kinetics of quartz precipitation at temperatures < 200°C [44] as well as silica polymerization [45] may allow silica concentrations to exceed quartz equilibrium. Equilibrium is reached much more rapidly at >300°C [3, 44]. However, the rate at which primary host rock minerals dissolve depends on the temperature, crystallinity, and pH [46]. In this study, the calculated rate of quartz precipitation or dissolution is not related to the kinetic rate of quartz-fluid reaction and instead represents the rate at which quartz is precipitated or silica is dissolved to maintain equilibrium with respect to quartz.

2.3. Porosity-Permeability Feedback

Rock permeability reflects changes in the pore volume resulting from quartz precipitation and the dissolution of silica from the host rock, as well as the effect of the brittle-ductile transition around the intrusion (e.g. [47, 48]). Assuming a unit volume of the porous rock, the porosity change resulting from quartz precipitation or dissolution is calculated at each time step according to where represents the density of quartz (2650 kg m−3). This is converted into a permeability change according to the relation of Weir and White [49]: where is the initial permeability (set to 10−14 m2 or 10−15 m2), is the initial porosity (set to 0.1, 0.05, or 0.025), and is a critical porosity set to half of the initial porosity, following White and Mroczek [10]. In our simulations, if porosity is reduced below , permeability is set to 10−16 m2, below which permeability fluid advection is minor and heat transport is controlled by conduction [15, 21, 22, 50]. Furthermore, porosity is not allowed to be reduced below 0.0005. As described below, porosity is reduced below the critical porosity and approaches 0.0005 only in systems with low initial porosity (0.025). However, these limitations are necessary to avoid extremely short time steps, which would lead to impractically long run times.

While permeability is discretized on the elements, quartz solubility and other fluid properties are nodal variables. Therefore, porosity changes were calculated at the element barycenter through averaging of the porosity changes calculated on the nodes belonging to that element. Permeability changes are quantified using the permeability reduction factor, the ratio of current permeability to initial permeability [10].

Temperature-dependent permeability is also included to simulate the transition from advection- to conduction-dominated heat transfer around the intrusion. We adopt the formulation of Hayba and Ingebritsen [15] where permeability decreases logarithmically with increasing temperature above a select brittle-ductile transition temperature, . Here, we set to 550°C, a value appropriate for basaltic rocks [51] that has been shown previously to be favorable for the development of supercritical geothermal resources [20]. Above , rock permeability is set to the temperature-dependent value only if it is lower than the corresponding quartz precipitation-dependent value. Below , permeability only depends on the porosity changes resulting from quartz precipitation and dissolution (equation 9). It should be noted that the initial porosity values are also applied to the intrusion and should be interpreted as primary porosity after full solidification.

2.4. Model Set-Up
2.4.1. One-Dimensional Simulation

To illustrate the effect of fluid flow and fluid density changes on quartz dissolution and precipitation, we conducted one-dimensional (1-D) fluid flow simulations of heating of liquid to supercritical conditions. The model domain consists of a 1000 m long vertical pipe divided into 101 nodes, resulting in an element size of 10 × 1 m. The rock properties in the one-dimensional test are described in Table 1. The heating path simulation assumed constant pressure and temperature boundary conditions of 30 MPa/300°C at the bottom of the pipe and 20 MPa/450°C at the top of the pipe. The top and side boundaries allowed no flow, while the bottom and top boundaries allowed flow into and out of the pipe, respectively. Initially, the temperature was set to 300°C throughout the domain, and a nodal heat source (term in (3)) of 7 W m−2 is applied to all nodes in the model domain, simulating the heat addition across the conductive boundary layer surrounding a magmatic heat source. This value for the nodal heat flux was selected because it produced the desired temperature profiles over a relatively short simulation time (e.g., 20 years). Note that porosity and permeability changes resulting from quartz precipitation and silica dissolution were not considered in the 1-D simulations.

2.4.2. Two-Dimensional Simulations of Cooling Intrusions

The two-dimensional models are intended to study processes during the main lifetime of a hydrothermal system after intrusion of a magmatic body. Thus, the intrusion process itself is not simulated and simplified as an instantaneous emplacement in host rocks with homogenous, isotropic permeability. Three main model configurations were considered (Figure 2): (i)High host rock permeability (10−14 m2) and a shallow intrusion emplaced at 2 km depth(ii)Intermediate host rock permeability (10−15 m2) and a shallow intrusion(iii)Intermediate host rock permeability and a deep intrusion emplaced at 5 km depth

A system configuration with high host rock permeability and a deep intrusion is not included as it was thought unlikely that permeability as high as 10−14 m2 would be found at depths of ~5 km [48, 52]. The shallow intrusion is centered at a 2.5 km depth with horizontal and vertical axis lengths of 2 km and 1 km, in a computational domain of 5 and 15 km in vertical and horizontal extents. The lengths of the axes are doubled in simulations with the deep intrusion, and the size of the domain is increased to 7 × 20 km.

Initially, the porous medium is saturated with water and thermally equilibrates with a basal heat flux of 0.15 W m−2, thought to be a reasonable approximation of a high heat flux in volcanic areas. The initial pressure distribution is hydrostatic, with a water table coinciding with the ground surface with flat topography, and lithostatic within the intrusion. The top boundary is fixed at atmospheric pressure. To allow hot fluids to vent, a mixed energy boundary condition was applied to the top [53]. In elements along the top boundary where the total volume flux is upwards, fluids vent at the temperature of ascending hydrothermal fluids; elements with a downward volume flux take in the fluid at a fixed temperature of 10°C. This simulates the effect of the recharge of cold meteoric water, sufficient to maintain a stable elevation of the water table. The other boundaries are no-flow boundaries, placed sufficiently far from the heat source that they do not affect fluid convection near the intrusion.

The primary porosity of the host rock is initially homogenous and set to 0.05 (hereafter referred to as “intermediate” porosity), 0.1 (“high” porosity), or 0.025 (“low” porosity). At the onset of simulation time, fluid within the intrusion is set to a temperature of 900°C and lithostatic pressure, describing an instantaneous intrusion of magma into the upper crust. The release of latent heat during crystallization is taken into consideration with a temperature-dependent rock heat capacity [15, 54], which increases linearly from 880 J kg−1°C−1 at temperatures below 750°C to 1760 J kg−1°C−1 at temperatures greater than 800°C.

2.5. Limitations of the Simulations

The objectives of this study are to assess the first-order control of host rock permeability and porosity on silica transport and quartz precipitation and to investigate the magnitude of permeability changes resulting from quartz precipitation. However, there are numerous limiting assumptions used in the model calculations that limit their bearing on natural systems. Fluid flow velocities are calculated based on a Darcian approach in a porous media, and the assumed porosity-permeability relationship (9) is based on theoretical studies of the effects of surface deposition and dissolution in an initial rhombohedral array of uniform spheres [49]. However, (9) is likely not fully appropriate for crystalline rocks with pore structures consisting of an interconnected network of pores and microfractures [55]. Rather, we assume this porosity-permeability relationship because it captures the essential dynamic of the initial rapid decrease in permeability for small amounts of quartz precipitation, followed by slower changes close to the critical porosity, as has been seen in experimental studies [9, 56], and moreover allows comparison with previous studies [10]. Furthermore, we assume that permeability is 10−16 m2 once porosity is reduced to the critical porosity (set to half of the primary porosity) and that further porosity reduction below the critical porosity (to a minimum value of 0.0005) does not lead to further permeability decrease.

One feature of our models is that the magmatic intrusion gradually transforms into host rock with progressive cooling by hydrothermal circulation. Thus, the permeability of a node initially located within the intrusion gradually assumes the initial host rock value as it cools from temperatures > 750°C to <550°C. However, thermal and mechanical stresses generated during the cooling of magmatic intrusions may result in pervasive microfracturing, porosity increase, and permeability enhancement [57]. Our models do not account for the effects of thermal cracking. The timescales of intrusion cooling, system development, and permeability changes resulting from quartz precipitation/dissolution are specific to the intrusion size of our model setup, the initial values of permeability and porosity, and the assumed permeability-porosity relationship. Thus, the transient evolution may be different in systems with smaller or larger intrusions, intrusions that undergo replenishment, or intrusions that expel large quantities of magmatic fluids.

In this study, we assume that host rock permeability is isotropic and initially homogenous. However, evidence from natural systems generally indicates that porosity and permeability decrease with depth due to increasing confining stress [27, 28, 48, 52, 55]. Moreover, natural systems show substantial permeability heterogeneity (e.g., [48]), controlled mainly by the distribution of major fractures and lithologies. A full investigation of the effect of initially heterogeneous host rock porosity and permeability is beyond the scope of this study.

In addition, we assume a simplified chemical system by neglecting the solubility of alteration minerals other than quartz. While the solubility of quartz is largely pH-independent in fluids with a pH below ~9 (e.g., [37]), a value above that generally measured in primary geothermal fluids [58], quartz solubility depends strongly on salinity [18, 39, 43] and the concentrations of volatile gases such as CO2 [42, 43]. Phase separation of NaCl-rich water at temperatures above the critical temperature of pure water may lead to brine formation and halite precipitation [59], which may significantly reduce porosity and permeability near the intrusion [60]. Thus, these models are mainly appropriate for geothermal systems that are recharged by meteoric fluids and therefore contain dilute groundwater, rather than seawater. Finally, as noted above, we do not simulate kinetically controlled silica dissolution or quartz precipitation. Rather, the assumption of local chemical equilibrium provides a rough estimate of quartz dissolution/precipitation rates.

3. Results

3.1. One-Dimensional Simulation of Quartz Precipitation at Supercritical Conditions

Figure 3 shows the results for the one-dimensional simulation of heating of liquid to supercritical vapor. Heat addition causes the temperature of initially 300°C liquid flowing towards the top boundary to increase with time (Figure 3(a)). After twenty years of simulation time, fluid expands from liquid-like densities (greater than the critical density of pure water, kg m−3) to vapor-like densities (less than ) ~0.5 km into the model domain (Figures 3(b) and 3(c)). The temperature increase and fluid pressure and density decrease cause quartz solubility to decrease from ~14 mmol kg−1 to ~1 mmol kg−1 (Figure 3(d)). The rate of quartz precipitation/silica dissolution, which is controlled by the changing quartz solubility gradient, is negative on the downstream side of the domain where quartz solubility increases, indicating dissolution, and is maximized at 2·10−8 mol m−3 s−1 upstream of the transition from liquid to vapor (Figure 3(e)).

3.2. Quartz Precipitation and Silica Dissolution around Cooling Intrusions

The transient evolution of a geothermal system structure controls the dynamics of quartz precipitation and silica dissolution near a cooling intrusion. Our previous approach [16] considers the temporal evolution of magma-driven geothermal systems in three stages: (i) an incipient stage, when one or more plumes of heated water ascend above the intrusive heat source, (ii) a main stage, when the boiling upflow plume has reached the surface, and (iii) a waning stage, when temperatures in the intrusive heat source have decreased below the brittle-ductile transition temperature and the intrusion is permeable. We focus our analysis on the main stage, when the most extensive supercritical geothermal resources form near the intrusion.

3.2.1. High Permeability System with a Shallow Intrusion

In a system with high initial host rock permeability (10−14 m2), intermediate porosity (0.05), and a shallow intrusion emplacement depth of 2 km, intensive quartz precipitation occurs within a thin zone surrounding the intrusion where circulating liquid is heated and expands to vapor-like densities (from now on called “supercritical vapor”). During the incipient stage (Figure 4(a)), boiling upflow plumes at >300°C with silica solubility > 10 mmol kg−1 ascend above the margins of the intrusion into overlying cold groundwater with low quartz solubility. In the main stage (Figure 4(b)), the system develops a central upflow with boiling conditions extending from the surface to the depth of the intrusion. In the waning stage, after the intrusion has been fully cooled by hydrothermal circulation (Figure 4(c)), boiling conditions remain at the surface and silica solubility in the upper 1 km of the system is like the main stage.

The highest rates of quartz precipitation occur within the supercritical zones surrounding the intrusion, where silica solubility decreases from ~13 mmol kg−1 to ~1 mmol kg−1 over ~0.1 km (Figure 4(d)). The liquid streamlines (blue) shown in Figure 4(e) demonstrate how the supercritical vapor contained within this zone is formed from the heating of liquid. Heated liquid ascends along the margins of the intrusion with vertical fluid mass fluxes exceeding 10−4 kg m−2 s−1 and then flows laterally and slightly downwards towards the top of the intrusion where the highest rates of quartz precipitation are found (~10−8–10−7 mol m−3 s−1). High rates of silica dissolution (10−9 mol m−3 s−1) occur nearby on the margins of the deep boiling zone where liquid is heated to temperatures > 350°C.

3.2.2. Intermediate Permeability System with a Shallow Intrusion

Systems with an intermediate initial host rock permeability (10−15 m2), intermediate porosity (0.05), and a shallow heat source show intensive quartz precipitation on the boundaries of spatially extensive zones of supercritical vapor (Figure 5). During the incipient stage, vapor-rich upflow plumes ascend above the margins of the intrusion and are separated by zones where liquid is heated to supercritical conditions as it flows downwards towards the top of the intrusion (Figure 5(a)). With continued evolution (Figure 5(b)), the system develops a boiling upflow plume centered on the intrusion that is underlain by a supercritical zone extending ~0.3 km above the intrusion. As liquid is heated from boiling to supercritical conditions on the boundaries of a spatially extensive supercritical zone, quartz precipitates at maximum rates of ~10−9 mol m−3 s−1 (Figures 5(d)–5(f)). Silica solubility is highest (~13 mmol kg−1) on the outer edges of the boiling zone and decreases rapidly towards the supercritical zone (Figure 5(d)). The highest vertical fluid fluxes (~2·10−5 kg m−2 s−1) and rates of quartz precipitation occur where this silica-rich liquid is boiled to dryness (Figures 5(e) and 5(f)). The rate of quartz precipitation decreases rapidly within the supercritical zone with increasing proximity to the intrusion. Silica dissolves at rates > 10−10 mol m−3 s−1 along the outer edge of the boiling zone as liquid is heated to ~300°C.

3.2.3. Intermediate Permeability System with a Deep Intrusion

Above a deeper intrusion (~5 km depth), quartz precipitates as fluid transitions from a liquid to a supercritical vapor, without any intermediate two-phase zone, as hydrostatic pressures greater than the critical pressure of pure water (22.06 MPa) prevail above the intrusion. During the incipient stage (Figure 6(a)), the ascending upflow plumes transition directly from supercritical vapor to single-phase liquid at ~3.5 km depth. After the two upflow plumes penetrate the surface (Figure 6(b)), the supercritical vapor zones extend from the depth of the intrusion to ~3 km depth. Intensive quartz precipitation occurs along the edges of the supercritical vapor plumes between 3 and 5 km depth. Silica solubility decreases from ~20 mmol kg−1 within the liquid circulating at 5 km depth to ~2 mmol kg−1 within the vapor plumes (Figure 6(d)). Fluid ascends vertically at rates of 2·10−5 kg m−2 s−1 throughout the supercritical zones (Figure 6(e)), while the highest rates of quartz precipitation (~10−9 mol m−3 s−1) are found at the transition between liquid and supercritical vapor (Figure 6(f)).

3.3. Porosity and Permeability Changes due to Quartz Precipitation and Silica Dissolution

The time-integrated effects of quartz precipitation lead to a maximum porosity decrease of ca. ±0.03 porosity units. Figure 7 shows how the spatial distribution of porosity changes with time in systems with an initial porosity of 0.05. The initial geometry of the intrusion is shown with dashed black lines. In an initial high permeability system (Figures 7(a)–7(e)), porosity increases along the outer edge of the intrusion and decreases within the solidified intrusion. Porosity is decreased more by quartz precipitation in intermediate permeability systems compared to high permeability systems. In an intermediate permeability system with a shallow intrusion (Figures 7(f)–7(j)), porosity is reduced from 0.05 to ~0.02 within a thin zone coinciding with the top of the intrusion. This zone is located near the base of the down-flowing liquid zones between the two upflow plumes during the incipient stage (compare Figures 7(f) and 7(g) with Figure 4(a)) and at later times is situated at the transition between supercritical and boiling conditions. In a system with a deep intrusion, porosity is reduced along the sides of the deep vapor plumes extending from 3 to 5 km depth (Figures 7(k)–7(o)). Porosity increases within the zone of liquid downflow in between the vapor plumes.

The porosity changes induced by quartz precipitation and dissolution have negligible effects (±10%) on the permeability of rocks with an initial porosity of 0.1 and moderate effects (±50%) on the permeability of rocks with an initial porosity of 0.05 but lead to order-of-magnitude permeability changes in rocks with an initial porosity of 0.025. Figure 8 shows the effect of initial porosity on the extent of permeability changes near intrusions by comparing the permeability reduction factor (ratio of current permeability to initial permeability) during the main stage for different system configurations. The snapshots representing a given configuration are compared at the same time.

In a system with high host rock permeability (Figures 8(a)–8(c)), permeability increases on the outer perimeter of the intrusion and decreases within the solidified intrusion, and the magnitude of permeability changes (both increase and decrease) becomes more significant as initial porosity is reduced. In high permeability rocks with initially low porosity (Figure 8(c)), permeability is more than doubled within a ~0.2 km thick zone along the top boundary of the intrusion. The zone of permeability increase is directly underlain by a ~0.1 km thick zone where permeability has been reduced by an order of magnitude.

In intermediate permeability systems, the vertical extent of low permeability zones in low primary porosity rocks increases with increasing intrusion emplacement depth. While an intermediate permeability system with a shallow heat source (Figures 8(d)–8(f)) shows ~0.5 km thick zones of permeability reduction near the initial top of the intrusion, a system with a deep intrusion shows zones of permeability decrease that extend up to 2 km above the intrusion (Figures 8(g)–8(i)). Zones where permeability has been more than doubled border zones of permeability reduction along the sides of the supercritical zones. While the most extensive permeability reduction is found along the sides of the upflow plumes ascending above the margins of the intrusion, the largest area of permeability increase occurs above the intrusion.

3.4. Results Summary

Whereas the maximum rate of quartz precipitation depends mainly on the initial host rock permeability, the magnitude and spatial distribution of permeability changes induced by quartz precipitation and dissolution depend on permeability, primary porosity, and intrusion depth: (i)In initially high permeability host rocks (10−14 m2), quartz precipitates intensively (at maximum rates up to ~10−7 mol m−3 s−1) within a thin zone (0.1–0.2 km thick) around the intrusion where circulating liquid is heated to supercritical vapor. In initially intermediate permeability rocks (10−15 m2), quartz precipitation occurs at lower rates (~10−9 mol m−3 s−1) along the boundaries of larger supercritical zones.(ii)Primary porosity is a key control on the magnitude of permeability change induced by quartz precipitation. While permeability changes are negligible (±10%) in host rocks with a primary porosity of 0.1 and moderate (±50%) in host rocks with a primary porosity of 0.5, quartz precipitation and dissolution may lead to order of magnitude permeability changes in host rocks with a primary porosity of 0.025.(iii)The amount of time necessary to form low permeability zones is controlled by host rock permeability. In initially high permeability host rocks, such zones form within timescales of hundreds of years; in intermediate permeability rocks, such zones form over thousands of years.(iv)In systems with a shallow intrusion (2 km depth), the thickness of low permeability zones along the top boundary of the intrusion is ~0.1 km in initially high permeability host rocks and ~0.3 km in initially intermediate permeability host rocks. However, above a deeper intrusion (5 km depth), zones of intensive permeability reduction are wider (~0.5 km) and may extend ~2 kilometers above the top of the intrusion.

4. Discussion

4.1. Quartz Precipitation during Supercritical Geothermal Resource Formation

Our simulations confirm that quartz precipitation is an inevitable consequence of the heating of circulating groundwater from liquid-like densities to vapor-like densities (supercritical vapor) near the intrusion. Following our previous study [20], we define supercritical resources as zones where fluid has a temperature and enthalpy greater than the critical values of pure water (374.1°C, 2.0183 MJ kg−1) [34] and is found in rocks with a permeability greater than 10−16 m2. Based on this definition, supercritical geothermal resources may develop at pressures both below and above the critical pressure of pure water (22.06 MPa). While the pressure of supercritical geothermal formation is close to ~20 MPa in systems with a shallow intrusion, fluids are heated to supercritical temperatures at pressures up to 50 MPa in a system with a 5 km deep intrusion.

Rates of quartz precipitation and silica dissolution at the transition from liquid to supercritical vapor increase with higher host rock permeability. In a high permeability system (Figure 4), vertical fluid fluxes are maximized (>2·10−4 kg m−2 s−2) within a ~50–100 m thick zone where liquid is heated from ~300°C to supercritical temperatures and quartz solubility drops from ~12 mmol kg−1 to >1 mmol kg−1, resulting in intensive quartz precipitation which occurs at rates up to ~10−7 mol kg−1 s−1. For a system with intermediate permeability host rocks and a shallow intrusion (Figure 5), fluid fluxes and quartz precipitation are optimized close to the critical point at rates of ~2·10−5 kg m−2 s−1 and ~10−9 mol kg−1 s−1, respectively. Circulating fluids above a deep intrusion attain higher fluid densities and maximum silica solubility due to the higher hydrostatic fluid pressures (Figure 6). Therefore, the change in silica solubility at the transition from liquid to supercritical vapor is less extreme compared to the systems with shallow intrusions and the temperature where fluid fluxes and rates of quartz precipitation are maximized shifts to higher values (>400°C) with increasing intrusion depth.

4.2. Effect of Quartz Precipitation on Fluid Flow and Heat Transfer Dynamics Near the Intrusion

Fluid flow dynamics near the intrusion are most strongly affected by quartz precipitation and dissolution in host rocks with low primary porosity. Since permeability changes in a system with a primary porosity of 0.05 are only moderate (Figures 8(b), 8(e), and 8(f)) and nearly negligible in systems with a primary porosity of 0.1 (Figures 8(a), 8(d), and 8(g)), these systems behave similarly to systems that do not include quartz deposition [20]. However, the order of magnitude permeability decrease resulting from quartz precipitation in rocks with a primary porosity of 0.025 significantly affects the spatial extent of supercritical zones and fluid circulation patterns near the intrusion.

In initially high permeability host rocks with low primary porosity, quartz precipitation and permeability reduction decrease fluid fluxes near the intrusion but increase the spatial extent of supercritical zones, compared to systems with primary porosity ≥ 0.05 (Figure 9). If primary porosity is 0.1 or 0.05, vertical fluid mass flux exceeds 2·10−4 kg m−2 s−1 along the margins of and/or above the intrusion (Figures 9(a) and 9(b)). Such high fluid fluxes only allow liquid circulating within ~0.1 km of the intrusion to be heated to supercritical conditions. However, the permeability reduction resulting from quartz precipitation in rocks with a primary porosity of 0.025 slows fluid circulation near the intrusion and thereby allows a greater fraction of circulating fluid to be heated to supercritical temperatures (Figure 9(c)). Vertical fluid mass fluxes > 2·10−4 kg m−2 s−1 are only found in zones ~0.3 km from the margins of the intrusion, where ascending fluid diverts around a low permeability zone above the intrusion (areas where permeability has been reduced by more than an order of magnitude compared to the initial host rock value, denoted by the thick black lines in Figure 9(c)).

The thicker low permeability zones resulting from quartz precipitation in systems with initially intermediate host rock permeability and low primary porosity compartmentalize flow between the supercritical and liquid regimes. For a system with a deep intrusion and host rocks with a primary porosity of 0.1 or 0.05, vertical fluid mass fluxes > 10−5 kg m−2 s−1 are found throughout much of the deep upflow zone above 5 km depth (Figures 10(a) and 10(b)). However, in host rocks with a primary porosity of 0.025 (Figure 10(c)), distinct liquid and supercritical zones with such high fluid fluxes are separated by a ~0.3 km thick low permeability zone (where permeability has been reduced from 10−15 m2 to 2·10−16 m2, denoted by the thick black lines, Figure 10(c)). This low permeability barrier coincides with the outer edges of the supercritical zone and encompasses an area with relatively low fluid mass fluxes (<5·10−6 kg m−2 s−1). Thus, in intermediate permeability host rocks, quartz precipitation may lead to physical segregation of supercritical vapor near the impermeable intrusion from cooler liquid circulating further away. Moreover, while the supercritical zones in Figure 10 extend from the depth of the intrusion to ~3 km depth in host rocks with a primary porosity ≥ 0.05, supercritical zones in host rocks with a primary porosity of 0.025 only extend to ~4 km depth. The reduced spatial extent of the supercritical zone in low primary porosity rocks is a result of significant permeability reduction and a lower rate of heat transfer from the intrusion.

Slower heat transport from intrusions emplaced in low porosity rocks because of quartz precipitation-induced permeability reduction is evident in Figures 810. Since the snapshots of a given system configuration show the same time in system development and the size of the impermeable intrusion increases as primary porosity is reduced, this indicates that lower primary porosity and more significant permeability reduction decrease the cooling rate of the intrusion. This is also seen in Figure 11, which shows the temperature at the center of an intrusion emplaced at 2 km depth (with a center at 2.5 km depth) in systems with intermediate host rock permeability. Quartz precipitation and dissolution delay the time at which temperature rapidly decreases from 700 to 200°C by 0.3 ka in a system with a primary host rock porosity of 0.1 compared to an identical system without quartz precipitation and dissolution. This temperature decrease occurs ~0.5 and ~1 ka later in systems with a primary porosity of 0.05 and 0.025, respectively, indicating that only in rocks with relatively low primary porosity do the permeability changes induced by quartz precipitation and dissolution significantly reduce the rate of heat transfer from the intrusion.

These results suggest that quartz precipitation is a key control on permeability and fluid flow near intrusions in host rocks with low primary porosity. Although quartz precipitation reduces permeability along the outer edges of supercritical zones, attractive supercritical resources with higher rates of mass and heat transport may also form below the low permeability zone. However, the development of economically attractive resources depends on the transition between brittle and plastic deformations around the intrusion. This study assumes that permeability decreases only above 550°C, which is supported by experimental studies showing brittle deformation in basalt to such high temperatures [51]. While it was previously thought that granitic rocks would be impermeable above ~360°C [11], a recent study showed that the temperature conditions of the elastic-plastic transition depend on the effective confining stress [61]. In granitic rocks, potentially exploitable supercritical geothermal resources with temperatures up to ~450°C form in rocks at effective confining stresses of 50 MPa, corresponding to a depth of ~3 km. Higher effective confining stresses at greater depths cause the maximum temperature of supercritical resources at permeable conditions to decrease to 400°C at 4.5 km depth (effective confining stress of 75 MPa). Thus, in granitic rocks with low primary porosity, the low permeability zone formed by quartz precipitation may directly border the intrusion, without an intervening zone with high supercritical fluid fluxes.

4.3. Timescales of Quartz Precipitation, Silica Dissolution, and Permeability Change

While permeability reduction due to quartz precipitation occurs over timescales of hundreds of years in host rocks with initially high permeability, permeability changes are slower in intermediate permeability host rocks and take place over thousands of years. The spatial configuration of zones of precipitation and silica dissolution change with time as the geometry of the impermeable heat source changes. As an intrusion cools, the number and locations of upflow plumes may change [16], and zones of intensive quartz precipitation initially along the sides of the intrusion transform into zones of dissolution. This causes clogging of pore space to remain limited. While zones along the sides of an intrusion experience early precipitation followed by later dissolution, locations above the top of the intrusion are more likely to undergo strong permeability reduction that restricts circulation.

In systems with initially high host rock permeability, intensive quartz dissolution along a flow path with increasing quartz solubility may produce a positive feedback leading to large increases in porosity and permeability. Figures 12(a)–12(g) show how fluid properties change with time at a point along the outer margin of the intrusion emplaced in host rocks with initial high permeability and low primary porosity (yellow star in Figure 9(c)). Following rapid heating from 200°C to ~370°C after emplacement of the intrusion (Figure 12(a)), quartz precipitates at rates near 10−8 mol m−3 s−1 (Figure 12(e)), porosity decreases from 0.025 to ~0.015, and permeability decreases from 10−14 to 5·10−16 m2 (Figure 12(g)). However, between 0.2 and 0.4 ka after intrusion emplacement, the location cools to <350°C (Figure 12(a)), close to the temperature at which quartz solubility is maximized, and vapor saturation decreases to zero (Figure 12(b)). Since quartz solubility increases along the fluid flow pathway, silica dissolves from the host rock at rates up to 4·10−8 mol m−3 s−1 (Figure 12(e)). Intensive dissolution sets up a positive feedback as liquid mass fluxes increase from ~1·10−4 kg m−2 s−1 to ~5·10−4 kg m−2 s−1 (Figure 12(c)), and porosity and permeability reach maximum values of ~0.17 and ~10−13 m2, respectively, by 0.4 ka after intrusion emplacement (Figure 12(d)). As the intrusion retracts inwards, causing the location to cool further to <200°C, quartz solubility decreases, and porosity and permeability remain nearly constant.

In high permeability systems, zones of intensive quartz precipitation may clog pore space within timescales of hundreds of years, particularly where liquid directly above the center of the intrusion flows downwards during heating to supercritical conditions. Figures 12(h)–12(n) show the temporal evolution of fluid properties at a point close to the apex of the intrusion (purple star in Figure 9(c)). Downflow of ~360°C liquid occurs throughout the first 0.1 ka following emplacement in this location at mass flow rates up to 5·10−4 kg m−2 s−1 (Figure 12(j)). Quartz precipitates at rates up to 2·10−8 mol m−3 s−1 (Figure 12(l)), reducing porosity to <0.0125 (the critical porosity for a system with primary porosity of 0.025) and permeability to 10−16 m2 (Figures 12(m) and 12(n)). Continued heating and additional porosity reduction below the critical value do not lead to further permeability reduction, as our model assumes 10−16 m2 is the minimum permeability resulting from quartz precipitation. These two orders of magnitude reduction in host rock permeability restrict supercritical vapor mass fluxes to less than ~10−5 kg m−2 s−1.

Permeability changes occur over longer timescales (thousands of years) in host rocks with initial intermediate permeability due to lower fluid mass fluxes and rates of quartz precipitation (Figure 13). Figures 13(a)–13(f) show the temporal evolution of a location on the perimeter of a deep intrusion (purple star in Figure 10(c)). Vertical liquid mass fluxes near 10−5 kg m−2 s−1 (Figure 13(b)) along a heating path where quartz solubility increases up to ~17 mmol kg−1 leads to quartz dissolution at rates of ~10−9 mol m−3 s−1 (Figure 13(d)), and cause porosity and permeability to increase from their initial values of 0.025 and 10−15 m2 to ~0.06 and 5·10−16 m2, respectively, over ~12 ka of system development (Figures 13(e) and 13(f)). Figures 13(g)–13(l) show the temporal evolution of a point located ~0.3 km closer to the center of the intrusion, where quartz precipitation at rates near 10−9 mol m−3 s−1 from supercritical vapor causes permeability to decrease from 10−15 m2 to 10−16 m2 within ~5.5 ka. Thus, compared to systems with initial high host rock permeability, systems with initial host rock permeability of 10−15 m2 show a lower magnitude of permeability increase within liquid heating zones, and clogging of pore space within supercritical vapor zones occurs over timescales of thousands of years rather than hundreds of years.

4.4. Comparison with Natural Systems

While the simulations have numerous limitations (see Section 2.5), they may help explain the observed spatial distribution of quartz within the deep roots of some exposed fossil geothermal systems. In the simulations, quartz precipitation is the most intensive inside of the host rock above the initial upper boundary of the intrusion (see black lines in Figure 10(c)). This result corresponds well to observations from the Skaergaard complex in Greenland [62], the Isle of Skye [63], and the Geitafell Volcano in Iceland [64], where quartz is abundant at distances of >0.2 km from the intrusive contact while the inner core of the gabbroic intrusion is relatively lacking in quartz. Quartz is found as part of a propylitic alteration assemblage including chlorite and actinolite that overprints higher temperature, pyroxene-bearing assemblages [6466]. Our simulations explain this observation as the consequence of cooling of and gradual permeability increase within the intrusion, which causes zones of quartz precipitation to shift to greater depths with time. The basal parts of sheeted dyke complexes within the Troodos ophiolite contain elongated (up to 1 km wide), pod-shaped, quartz-rich epidosites [67], similar to the quartz-rich, low permeability zones that form in our simulations with initially high permeability host rocks (see Figure 9(c)). However, the dynamics of quartz precipitation in subseafloor hydrothermal systems containing seawater as the pore fluid will differ from systems containing pure water because of phase separation at temperatures above the critical point of pure water, brine formation, and halite precipitation [18, 59, 60].

It should be noted that the dynamics of quartz precipitation from hydrous intrusions that expel large quantities of magmatic fluid are beyond the scope of this study. Observations from porphyry ore deposits suggest that fluid expulsion causes quartz precipitation within stockwork veins (e.g. [68]). As the fracture networks above a fluid-rich magma fill with quartz and permeability decreases, fluid pressure increases to near-lithostatic, leading to subsequent hydrofracturing and permeability increase [32]. However, in our simulations, the convecting fluids are not pressurized to near-lithostatic values even within the low permeability zones. In addition to the absence of a fluid source, this result may also be a consequence of the assumption that the quartz precipitation does not reduce permeability below 10−16 m2. Previous studies suggest that significant fluid pressure build-up in the presence of a fluid source requires lower permeability [69].

5. Conclusions

The simplified reactive transport models of fluid flow around upper crustal intrusions presented in this study show how fluid fluxes, the rates of quartz precipitation and dissolution, and the magnitude and timescales of permeability reduction depend on host rock permeability, porosity, and the intrusion emplacement depth. Relatively high fluid fluxes at the transition from liquid to supercritical vapor above the intrusion in initially high permeability host rocks (10−14 m2) result in quartz precipitation at rates of ~10−8 mol m−3 s−1 (corresponding to volumetric changes of ~100 cm3 quartz per cubic meter of rock per year). Within timescales of hundreds of years, a relatively thin (<0.3 km thick), quartz-rich, low permeability zone develops above the top of the intrusion. Lower mass fluxes in initially intermediate permeability (10−15 m2) host rocks result in lower rates of quartz precipitation (~10−9 mol m−3 s−1) along the edges of larger supercritical zones. However, since the intrusion cools slower in intermediate permeability host rocks, there is more time to develop larger (0.3–0.5 km thick), vertically extensive silicified zones, which may extend ~2 km above intrusions emplaced at 5 km depth.

This study identifies host rock primary porosity as a key control on the magnitude of permeability changes induced by quartz precipitation and dissolution. Quartz precipitation results in an order-of-magnitude permeability decrease in host rocks with primary porosity of 0.025, a realistic value for natural rocks at depths > 2 km. However, permeability changes are only moderate in host rocks with a primary porosity of 0.05 and largely negligible in host rocks with a primary porosity of 0.1. The degree and spatial distribution of permeability changes above the intrusion reflect the transient configuration of upflow plumes as the intrusion cools and changes geometry. Zones of quartz precipitation and permeability reduction along the sides of the supercritical zones are bordered by zones of quartz dissolution and permeability increase, where liquid is heated to temperature ~350°C and quartz solubility is maximized. Thus, our models predict a lower magnitude of permeability reduction compared to previous models that assume a fixed temperature basal boundary condition (e.g., [10]), because zones of precipitation may transform into zones of dissolution as the intrusion cools. This limits the clogging effect, except in rocks with low initial porosity, in which low permeability zones form more rapidly.

The reduced permeability zones above the intrusions may segregate supercritical vapor close to the intrusion from circulating liquid further away, where heated liquid dissolves silica from host rock minerals and permeability is increased. Thus, while quartz precipitation-induced permeability reduction reduces fluxes of supercritical vapor and the rate of heat transfer from the intrusion, the effect is partially counteracted by increased liquid fluxes within the higher permeability zones. This study suggests that quartz precipitation and subsequent permeability reduction does not prevent the development of supercritical geothermal resources, even in rocks with low primary porosity. However, fluid flow dynamics within supercritical zones will also depend on other porosity/permeability creation and destruction processes associated with intrusion emplacement and cooling, such as magmatic fluid release [32], thermal cracking [57], inelastic compaction [31], and the precipitation of secondary minerals other than quartz [66]. The simulations presented in this study do not address the complex interactions between hydrothermal convection, quartz precipitation and dissolution, and these other factors, which should be investigated in future studies.

Data Availability

Unprocessed model output data is provided as a zipped archive of vtk files. The CSMP++ computing platform is co-owned by ETH Zurich, Heriot Watt University, and Montanuniversität Leoben. The source code is not distributed freely; precompiled code libraries may be obtained via one of the owning institutions.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors thank the COTHERM partners, particularly Bruno Thein, Andri Stefánsson, Matylda Hermanska, and Dimitri Kulik. This manuscript benefited greatly from discussions with Matthew Steele-MacInnis and Szandra Fekete, as well as the comments of two anonymous reviewers. This study was funded by the Swiss National Science Foundation [Grant no. CRSII2 141843/1, Sinergia COTHERM].

Supplementary Materials

Zipped archive of model output vtk files from simulations shown in this paper, including temperature, fluid pressure, vapor volumetric saturation, total mass flux (vertical fluid mass flux), porosity, permeability, mass quartz equilibrium total (concentration of dissolved SiO2 in mol/kg), and rate of quartz precipitation. (Supplementary Materials)