Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Earth Sci., 12 January 2026

Sec. Georeservoirs

Volume 13 - 2025 | https://doi.org/10.3389/feart.2025.1706071

Physicochemical kinetic reactive transport modeling in quartz arenite: implications on geothermal potential

  • Department of Geology and Geography, West Virginia University, Morgantown, WV, United States

Kinetic reactive transport modeling has been performed on the quartz-dominant Tuscarora Sandstone. Physical parameters including temperature, porosity, and permeability, as well as chemical parameters including pH, dissolved ions, and mineral precipitation and dissolution have been explored as a function of flow rate. This approach has shown strong spatial relationships between fluid-rock interactions and flow, particularly with respect to extensive reservoir cooling. At 10 L/s, the reservoir is cooled to ∼75 °C after 50 years and is nearly entirely equilibrated with the working fluid at 30 L/s. There is a slight increase in porosity immediately near the injection well likely due to reservoir dissolution. Simulation results in the reservoir block model show negligible secondary mineral precipitation at all modeled injection rates. We suggest that the proposed model can be utilized in evaluating the geothermal potential of low porosity quartz arenite reservoirs by predicting spatial temperature variations as well as observed fluid-mineral assemblages in field settings.

1 Introduction

To address the ever-increasing demand for clean energy, unconventional hydro(geo)thermal reservoirs must be sought out. A first-principle approach in assessing the potential viability of a prospective geothermal reservoir is to develop a simplified approximation of the likely fluid-rock interactions by implementing matrix-transport simulations. A similar note has been made by Wanner et al. (2014), who mention that many geothermal reactive transport models (RTM) consider only fluid flow and heat transfer. More recently, a series of RTM studies have evaluated fluid-rock interactions, including mineralization in geothermal systems (Erol et al., 2022; Holmslykke et al., 2023; Benjakul et al., 2020). Rock permeability is a critical parameter in assessing a potential geothermal reservoir as it plays a key role in controlling the subsurface flow rates, thus affecting the quantity of heat extracted and returned for operational use at the surface (Donselaar et al., 2015). Permeability, in turn, can be impacted by mineral dissolution and precipitation reactions (which are themselves controlled by the composition and flow rate via kinetics caused by fluid-rock residence time) (Figure 1). Therefore, a more complete RTM for assessing a potential geothermal reservoir should incorporate additional aspects such as mineral precipitation or scaling and changes to reservoir porosity and permeability.

Figure 1
Flowchart illustrating the general hierarchy of mineral processes. It begins with

Figure 1. Schematic flow map of major parameters in geothermal exploration and RTM as they relate to pressure (P), temperature (T), and composition (X). Note that flow rate is the variable being altered and that reservoir T and brine X are affected.

In the United States, geothermal energy is traditionally restricted to the country’s western portion where the resource is more readily accessible (i.e., located at shallower depths) (Tester et al., 2021). In the Eastern United States, significantly less implementation and reservoir identification have been performed due to lower heat flux and geothermal gradients. A few studies have been conducted to evaluate the geothermal potential of sedimentary reservoirs in the Appalachian Basin (Camp et al., 2018). Despite the elevated geothermal gradient found in the northern Appalachian Basin (Blackwell et al., 2011; McCleery et al., 2018), the prospective reservoirs in this region are largely understudied. Recently, West Virginia University, in collaboration with industry partners and the National Energy Technology Laboratory of the Department of Energy, has drilled an exploratory geothermal well near the university campus in Morgantown, WV (Monongahela County, see Figure 2). Referred to as the Geothermal Science Energy and Environment Laboratory (GSEEL) well site, efforts of critical downhole data collection aid in reservoir selection for possible future deep direct-use (DDU) at depths above 10,000 ft (>3 km). One of the proposed DDU reservoirs in the Appalachian Basin is the Tuscarora Sandstone (TS) (Garapati and Hause, 2021). This formation was deposited in a nearshore (‘submerged deltaic-estuarine’) environment and is characterized by two defining layers; the lower segment is a graywacke, while the upper segment is a quartz arenite (Folk, 1960). Sediments of the Lower Silurian TS were transported by fluvial processes as evidenced by cross-bedding (Yeakel, 1962). Near the GSEEL site, this formation is not faulted, although it appears along the eastern fold limb of the Chestnut Ridge Anticline (Figure 2). Fractures increase to the east due to steep forelimb folds in comparison to west-dipping fold limbs near the GSEEL site (Krame, 2013). The overlying Rose Hill Formation is comprised of shale with intermixed layers of siltstone and hematite-bearing sandstone (Folk, 1960). Noted for its robust utility as a carbon sequestration caprock, the Rose Hill Formation serves to constrain potential injected TS geofluids as well (Krame, 2013).

Figure 2
Geological map and cross-sections depicting the Chestnut Ridge and Hiram anticlines in West Virginia. Panel A shows locations in Marion and Preston Counties with key features marked. Panel B presents a cross-section showing stratigraphic layers, including the McKenzie Limestone and Tuscarora Sandstone, with vertical exaggeration noted. Panel C details the Tuscarora Sandstone flow direction and RTM block characteristics. The GSEEL site is indicated with a star, and distance scales are provided in miles and kilometers.

Figure 2. Geologic and geospatial map of the proposed Tuscarora Sandstone. The spatial map and cross-sections are modified from Ryder et al. (2009). (A) Location of the GSEEL site in relation to the cross-section. (B) Cross section along the transect (orange semi-diagonal line). (C) Magnified portion outlining the relative location of the Tuscarora Sandstone at the ∼3.2 km depth (10,385 ft). The horizontal length of the red box is to scale and represents the 1 km RTM block model visualized in Figure 3. The flow arrows are also in agreement with those in Figure 3 and are not to scale.

Although the TS (McCleery et al., 2018; McDowell et al., 2020; Garapati et al., 2019) has been identified as a possible geothermal reservoir spatiotemporal simulations approximating fluid-rock interactions, and indirectly transmissivity, are undetermined. Compared to hydrothermal systems which recycle substantial quantities of water efficiently with production and injection rates often to 103 tons/hr (Kaya et al., 2011), the efficiency and throughput of manmade reservoirs are typically significantly lower (Lu, 2018). To this end, a combined physicochemical RTM is utilized to estimate changes to the physical reservoir and its chemical composition with respect to time and lateral distance from the injection well in the Tuscarora Sandstone. RTM is noted for its utility in exploring relevant subjects such as porosity development (e.g., Schott et al., 2009) and the effect of mineral dissolution or precipitation on porosity evolution (e.g., Abarca et al., 2019). In addition, the interrelated nature of porosity and permeability are affected by chemical reactions; namely, those involving mineral precipitation or dissolution, which are dependent upon (among other variables) ionic strength and fluid total dissolved solids (TDS). These considerations make RTM ideal for evolving systems, especially those over relatively large spatial scales in which metastable equilibrium is disrupted by human-induced physicochemical changes (e.g., pumping, antiscalant additives).

This modeling work will primarily focus on fluid-rock interactions at varied injection flow rates in a quartz arenite. The model setup is based on the recent interest in the effect of injection flow rate in fluid-rock interaction (Erol et al., 2022). The quartz arenites contain little Al-bearing minerals, therefore the likelihood of clay swelling and secondary pore-clogging clay deposits are lessened. In addition, scaling plays a significant role in determining the lifespan of the geothermal reservoir (Moeck, 2014). Sandstone is often identified and evaluated as a potential geothermal reservoir medium, with specific examples including the Frio Formation in Texas, United States (Bebout et al., 1978), the Triassic/Jurassic sandstones in Denmark (Weibel et al., 2020), the Delft Sandstone in the Netherlands (Donselaar et al., 2015), and the Ripon Formation in South Africa (Campbell et al., 2016). Although the TS has been selected as a case study, general conclusions made from RTM observations can be applied to other quartz-dominated, low-permeability sandstone reservoirs. In addition, the RTM results below (where fluid is only transferred through existing pores) may identify general constraints applicable to prospective quartz-rich sandstone reservoirs slated for hydraulic fracturing. If controls of specific fluid-rock interactions are documented, these baseline observations have application to geothermal exploration.

2 Methodology

There are multiple programs capable of performing RTM, but we have chosen to invoke Geochemist’s Workbench v.17 (GWB) for equilibrium analysis of the TS, given its strength in addressing geochemical problems. In the RTM simulations, the unmodified LLNL thermo. tdat database was used for thermodynamic and kinetic calculations by GWB and X2t. Although GWB does not incorporate an equation of state (Shevalier et al., 2014), phase diagram analysis of the pressure and temperature conditions associated with the TS does not suggest that this will present adverse effects on the RTM results. At the pressure and temperature of the model, the fluid plots in the liquid domain below both the critical and boiling points-suggesting that both volumetric changes (for water at 90.56 after interpolating between 25 and 100 there is a volume increase of ∼3.5% [Parr. Safety in the Operation (2012); and reference therein] and phase mixing are minimal and a single liquid phase is expected. Geophysical properties of the TS are poor for gas accumulation, thus permitting the exclusion of an in situ CO2 and hydrocarbon presence in modeling (Ryder, 2004). The choice of GWB is reasonable for RTM of the TS because, in previous work, GWB, PHREEQC, and TOUGHREACT softwares have yielded similar results at temperatures below 100 (Klunk et al., 2021) with the ensuing model being simulated at a maximum temperature of ∼91 . Two fundamental aspects which are incorporated in RTM codes are the reduced time-dependent conservation of mass (Equation 1) and heat transfer (Equations 2 and 3, e.g., as utilized in GWB). For the i th species in a fluid;

tφsπiCi+·Ji=Ri1(1)

where φ, sπi, Ci, (e.g., mol/m3) ,Ji, (mol/m2s), and Ri are porosity, fractional pore volume saturated by fluid containing i, the concentration, and transport flux (combination of advective and diffusive transport), and the ‘source/sink’ term respectively (Lichtner et al., 1996). The Ri term is also sometimes referred to as the ‘total reaction rate’ and is in units of mol/m3 porous media/time (SteefelC. I., 2009). Two-dimensional heat transfer in GWB (Bethke et al., 2024) is given as

qHx=KHTxforx,and(2)
qHy=KHTyfory(3)

where qH (cal/cm2/s), KH (cal/cm/s oC) and the differential (oC/cm) are the conductive heat flux, thermal conductivity, and temperature-distance slope respectively.

2.1 Model input parameters and setup

2.1.1 Tuscarora Sandstone rock composition and physical properties

The mineralogy and physical parameters of the TS used in the RTM simulations, are given in Table 1. TS is comprised predominately of quartz [95% Borhara and Onasch (2020)] and has an average volumetric quartz abundance of 91% (Houseknecht, 1988) along with minor lithic clasts. The GSEEL depth of 3.2 km indicates that the location within the TS would be the quartz-rich upper segment (Figure 2). Using four temperature logs over a time period of approximately 2 months, the estimated temperature of the TS at ∼3,200 m (local depth in north central West Virginia) is 90.56 (Carr, T., 20231); similarly, the depth corresponds to a pressure of 63.01 MPa (pers. comm. Northeast Natural Energy, 2023). The permeability of the Tuscarora Sandstone is low (Heald and Anderegg, 1960), with varied porosity depending on the depositional environment (Castle and Byrnes, 2005). Locally, the Medina-Tuscarora Sandstone is located in the fluvial/estuarine depositional setting, with Klinkenberg fluvial permeability ranging from 0.002 to 450 mD with an average fluvial porosity of 8.1% (range of 1.9%–18.1%) (Castle and Byrnes, 2005). This permeability is in good agreement with the matrix permeabilities obtained by McDowell et al. (2020), who provide a range between 1 and 2 mD (McDowell et al., 2020). As noted in the GWB program, the permeability is calculated in accordance with these basic input data (Bethke, 2008). Permeability is not static because fluid-rock interactions can enhance or degrade permeability (Yuan et al., 2017) as the fluid travels laterally and over time, indicating that this dynamic approach taken by GWB is robust. Finally, the thickness of the TS ranges from 30 to 100 m (Borhara and Onasch, 2020) and is in good agreement with our average calculated thickness of 61 m (1σ = 32 m) using data from McCleery et al. (2018). The TS is a good potential candidate for geothermal exploration because it is vertically bound between two shale units, the older Juniata Shale below and the younger Rose Hill Shale and Siltstone above (Youse, 1970). Due to the prevalence of quartz in the TS, the heat capacity of the bulk reservoir is well approximated by obtaining the heat capacity of quartz at reservoir temperature. Using empirically obtained heat capacity values for 300 K and 400 K found in Robie and Hemingway (1995) and interpolating to 363.71 K (90.56 ), we find a value of 835.3 J/kgK.

Table 1
www.frontiersin.org

Table 1. Mineralogy and accompanying physical properties for the Tuscarora Sandstone used in the RTM simulations.

Equilibrium modeling is usually considered for determining the stability of different mineral phases in any geochemically evolving system. However, this approach is less relevant when the temporal extent for mineral-brine interaction is too short for complete equilibration of all the phases. Additional constraints such as rates of mineral dissolution/precipitation found in experiments or field observations do not often conform with complete equilibrium approaches. To avoid any overestimation for the geochemical perturbation of the studied system, we have incorporated known kinetic parameters for dissolution of the primary minerals (Table 2). The reprecipitation (if any) of minerals are considered in equilibrium with the brine as the available kinetic data for precipitation comes with several uncertainties. Reactions are described in Supplementary Table S1 (Supplementary Material). The dissolution rate (r, in mol/cm2/s) of any mineral phase is calculated using

r=S*kx*expEaxRTα*aij*1QK.(4)

Table 2
www.frontiersin.org

Table 2. Kinetic parameters (calculated for 90.56 °C from Equations 1A and 1B) associated with different mineral phases. Parameters at 25 °C are taken from literature.

S is the available reactive surface area (cm2/g) of the mineral, k (x) is the dissolution rate constant (mol/cm2/s) for different mechanisms (i.e., k (a) for acid, k (n) for neutral, and k (b) for base). The activation energy is denoted as Ea and is in units of kJ/mol. Note that k (x), Ea(x) also pertains to the three pH range mechanisms (i.e., Ea(a), Ea(n), and Ea(b). Moreover, a(i) is the activity of H+ or OH species depending upon the pH mechanism and reaction order of j. The saturation index of the mineral component in question is given by Q/K and indicates the deviance from equilibrium (where Q is the ion activity product of the reaction, and capital K is the equilibrium solubility at the corresponding temperature of the reaction). Lastly, the units of temperature (T) are Kelvin. These parameters are summarized in Table 2. For input into GWB, the average of the kinetic rate constant is used as the neutral-basic mechanism rate constants are several orders of magnitude lower than acid mechanism and therefore have negligible influence in this study.

2.1.2 Brine and working fluid composition

Tuscarora brine chemistry was not retrieved from the GSEEL well. As such brine data from a depth of ∼3,340 m in Centre Co. in central Pennsylvania (Dresel and Rose, 2010) was substituted to simulate undiluted formation water. The brine and working fluid compositions used in the RTM simulations are provided in Table 3. For cations without defined redox states (Table 3), the charges ascribed for the RTM are Cu+, Fe2+, Mn2+, and Zn2+. Since there are tremendous discrepancies between the number of elements, aqueous species, and minerals available in thermodynamic databases and the desired activity models, we have chosen to perform RTM using a 2× dilution of the TS formation water, which will serve as the initial fluid component in the TS reservoir. This approach ensures that the ionic strength is low enough (SpecE8 in GWB calculated as 2.3 mol/kg for the 2× diluted solution) for the B-dot activity model used by GWB while still also ensuring that the LLNL database can be used, which contains the most extensive thermodynamic availability of species. Unfortunately, the Pitzer activity model and database are inappropriate because this model does not incorporate either Al or Si and since the TS contains >90% quartz and minor Al-bearing phases, it is a poor choice in this instance. Also, because injected water is highly diluted brine anyway [approximately 1:100 (Pilewski et al., 2019)], it is reasonable that the relative abundance of the brine to diluted brine is minimal. Therefore, the 2× diluted TS formation water was used to serve as the initial solution basis for RTM in the X2t module of GWB. After converting from w:v TDS (mg/L) (Dresel and Rose, 2010) to w:w TDS (mg/kg), we find that the 2× TS brine contains 131,132 mg/kg. The working fluid is simply a 1:100 dilution of the TS brine, and thus, the density was calculated as a two-member mixture using

ρWF=0.01ρbrine+0.99ρpw,(5)

where ρWF ρbrine, and ρpw are the densities for the diluted working fluid, TS brine, and pure water, respectively. Equation 5 yields a working fluid density of 1.00181 g/mL. The initial silica concentration to be used at the start of the RTM simulations was found by calculating the temperature-dependent equilibrium of amorphous silica (H4SiO4) following the equation given by Gunnarsson and Arnórsson (2000) as

logKH4SiO4=8.476485.24T2.268E6T2+3.068*logT.(6)

Table 3
www.frontiersin.org

Table 3. Undiluted brine alongside the 2× initial solution and 100× diluted working fluid compositions used in the RTM simulations. Note that the pH of the working fluid was unchanged. Also, the densities of the diluted solutions are weighted mixtures of the undiluted brine with pure water. The absolute uncertainties found in (Dresel and Rose, 2010) for the undiluted TS brine are simply added to the diluted fluid compositions (except for Cl).

The above equation with temperature (T) in K is the equilibrium solubility for SiO2 + 2H2O = H4SiO4 (Gunnarsson and Arnórsson, 2000). Considering that the TS contains clay minerals, Al3+ must be accounted for in the brine. Using oil and gas brine data for the Third Sand from the Venango Group [(Chapman et al., 2013) and references therein], the Al3+ for the TS was found by simple ratio determination using Equation 7 below. For the Venango Group (VG), the concentration of Al3+ of the TS was estimated using

Al3+TDSVG*TDSTS=AlTS3+,(7)

where all units are in mg/L. It is important to note that Al+3 is sensitive to pH, and the above assumption relies on both brines having the same pH. For the initial basis fluid, Equation 7 gives an Al+3 concentration of ∼0.06 mg/L.

While the dissolved oxygen (DO) concentration is given for the TS brine, the DO of the working fluid was calculated using Henry’s Law at 1 atm (101,325 Pa) and at reinjection temperature of the working fluid after mixing with formation water (43.1 ). Assuming equilibration with standard atmospheric composition (Sander, 2015; Merkel et al., 2005) and yields a dissolved O2 concentration of 5.88 mg/L. Since working fluids are typically sourced from nearby freshwater in contact with the atmosphere, we believe our calculation to be a reasonable estimate. The injected freshwater will have considerably greater DO than the static formation water of the TS (Table 3), and it would also be expected that the Eh of the recirculating water would be commensurately high. Due to the relatively low temperatures and continual input of oxidized working fluid from the surface, reactions of pertinent redox species are decoupled as electric equilibrium is not expected.

2.1.3 Fluid flow and block model setup

The final tabulation of data used to input for the RTM is given in Table 4 and pertains to flow dynamics and reservoir block model geometry. The temperature and pressure-dependent self-diffusion coefficient of H2O was calculated using the equation in Harris and Woolf (1980), which has been slightly rewritten as

DS109m2s=expA0+i=13xiA2i1+A2i*1000T+Ci*1000Ti.(8)

Table 4
www.frontiersin.org

Table 4. Fluid flow and block model parameters for the RTM simulations.

The fitting terms for A0, A1, A2, A3, A4, A5, and A6 are 3.425150, −0.6275*10–3, 0.202474*10–3, 0.114172*10–6, −0.447466*10–7, 0.450105*10–11, and 0, respectively. The fitting terms for C1, C2, and C3 are 0.623898, −0.416757, and 0, respectively. Temperature (T) is in Kelvin. At the TS reservoir conditions, Equation 8 yields DS of 6.98*10–9 m2/s. The uncertainty associated with Equation 8 is % but has only been validated up to 333 K (Harris and Woolf, 1980), although experimental H2O self-diffusion results up to 498 K are reported to be <5% (Krynicki et al., 1978) and appear to be in good agreement with those in (Harris and Woolf, 1980). Therefore, we take the extrapolated uncertainty of Ds at 363 K to be <5%. Another important mass transport parameter is dispersivity, where previous work has shown that the dispersivity for a sandy stratified reservoir is 0.7 cm and is irrespective of distance (Pickens and Grisak, 1981).

Water-rock ratios (W:R) for sedimentary brines are low, especially for fluids with high TDS. However, regressing the TDS vs. W:R relationship curve to TDS of the 100× diluted working fluid gives W:R of 3.06 (Bryndzia and Fay, 2016). Therefore, for each 1 L of the working fluid at 1.00181 kg, there would be an accompanying 0.328 kg of sandstone reservoir. Thus, the constituent masses for TS minerals may be determined as they are titrated with the working fluid over the simulation time period (Table 4).

The temperature of the injection well fluid was calculated to be 43.10 . This calculation was performed by using the average injection temperature (reported range of 45 to 148 ) and the difference between the injection fluid temperature and the reservoir well temperatures for hot water geothermal systems (Kamila et al., 2021). For example, for hot water systems, the average injected temperature is 75 (Tinj¯) and the average difference between the reservoir temperature and the injected temperature is 82.6 (Tresinj¯) yielding an average reservoir temperature of 157.6 (Tres¯) (using values for Tinj¯ and Tresinj¯) from Kamila et al. (2021). Then, for the TS with an assumed reservoir temperature of 90.56 (TTS res), we find

TTSinj=Tinj¯Tres¯*TTSres=75C157.6*90.56=43.10.(9)

While the calculated temperature is lower than the lower bound on the range, this is directly the result of the low temperature of the proposed TS, not a methodological error. Simulations were also run considering the effect of 2.14 μW/m3 heat flow. This value is the radiogenic heat generated from gneiss (Khutorskoi and Polyak, 2016). Basement rocks underlying the TS in northern West Virginia include gneiss (Grenville Province) and granite (Chilhowee Group) (Ryder et al., 2009).

Conducting separate RTM simulations at multiple injection rates allows for a basis for comparison. A conservative injection flow rate of 20 L/s (20.04 kg/s or 16 tons/hr at 1.00181 kg/L) to simulate injection into the TS reservoir. Additional RTM simulations at alternate injection flow rates of 10 L/s and 30 L/s were also performed using otherwise identical input parameters The final parameter found in Table 4 is the reservoir impedance, and to maintain economic viability, a reservoir impedance pressure of 0.1 MPa is recommended (Murphy et al., 1999). The RTM calculations were run over a simulated time period of 50 years, as this corresponds to the time lengths noted for operational lifespan (Sun et al., 2023). Fifty years has been used in fractured simulations of low porosity reservoirs (Gong et al., 2019), suggesting that this is a suitable time estimate. Nonetheless, operational time lengths vary widely from 30 years (Augustine, 2016) to 135 years (Lowry et al., 2018). Given the wide range of operating times, a range of injection flow rates used in the RTM simulations is appropriate to establish general trends and reduce model simulation time bias.

The X2t module within GWB performs 2-D RTM simulations and visualizations of 3-D structures and flow pathways (Figure 3). The fluid travel presented in a map or overhead view represents a single plane in space and that fluid flow would be above and below the illustrated view. Only the reservoir height zn of 61 m, and well location of at x = 0 m, y = 0.5 years m was used consistently throughout RTM simulations as the x and y dimensions are adjusted to explore the changes to the physicochemical system as a function of distance.

Figure 3
Diagram of a rectangular prism with labeled axes: Xn in orange, Yn in green, and Zn in purple. A well marked as

Figure 3. Cartoon box model of a hypothetical subsurface reservoir (not to scale). The reservoir length (xn), width (yn), and height (zn) along with the injection well location (0.5 years) and the number of nodes (nn) are input into X2t of GWB for RTM. The flow of the water (denoted by blue arrows) is forced to travel from the left of box to the right side of the box which is open to continuous flow.

3 Results and discussion

Results presented below and the ensuing discussion pertain to data obtained after the 50-year simulation time. Results for intermediate timepoints at 10, 20, 30, and 40 years of elapsed simulation time for temperature, pH, porosity, permeability, mineralization (as mol/m3), and saturation index (as Q/K) are provided for each injection rate. In addition, the saturation index for amorphous silica is provided as Supplementary Material is also provided at 0- and 50-year simulation endpoints. It is of course acknowledged, that the observed results after 50 simulated years represent the final evolved state of the system after the entire time interval. The effect of simulations considering the incorporated heat flow is negligible (for example, the average absolute difference between temperature with and without heat flow is 1.1 × 10−3 °C with 1 σ of 3.9 × 10−4 °C for all nodes at 10 L/s, and 5.6*10–4 °C with 1σ of 3.0*10–4 °C for all nodes at 30 L/s) and comparative data are presented for each spatial node of the RTM block at each x,y coordinate for both 10 L/s and 30 L/s injections. These data are provided as Supplementary Material. Additionally, the maximum difference between albite mineralization with and without the external heat source is 8.00*10–8 mol/m3 for all nodes at 10 L/s, and 2.8*10–8 mol/m3 for all nodes at 30 L/s. These data are presented as Supplementary Material.

3.1 Reservoir temperature evolution

The effect of injection flow rate on TS reservoir temperature after 50 years across the 1 km (Erol et al., 2022) simulation area is illustrated in Figure 4. When the dilute working fluid is introduced to the reservoir at a rate of 10 L/s (Figure 4A), the ∼91 °C reservoir is cooled to temperatures <75 °C at a distance of 750 m from the injection well. Thermal breakthrough (when reservoir T is 90% of initial) is not achieved at 1 km from the injection well after 50 years, although it has been achieved at a distance between 750 m and 850 m (T ≈ 82 °C). Distances beyond 750 m are warmer (up to 86 °C at 1 km, and corners of RTM block ∼90 °C at 1 km) but all grid nodes experience some varied amount of cooling. As the injection flow rate is increased to 20 L/s and beyond (Figures 4B,C), the reservoir temperature becomes progressively cooled, and outward migration of the cooled margin occurs. At 20 L/s injection, thermal breakthrough (∼81.5 °C) achieved after 30 years at a distance of 950 m from the injection well, while only 20 years is required for 30 L/s—also at 950 m (Supplementary Material). At an injection rate of 30 L/s, the reservoir temperature is nearly entirely equilibrated with the working fluid, even at distances ∼1 km from the injection well site. At this injection rate, the reservoir is cooled more rapidly than the working fluid is heated. To evaluate the lateral extent of the cooling at 30 L/s, the RTM simulation was rerun using a 2 km × 2 km x 61 m block with commensurate updated masses for both brine and primary wallrock minerals. After a period of 50 years, the temperature at 1.9 km from the injection well had decreased by only 1.1 °C, to 89.5 °C (averaging data at y = 0.9 km (89.2 °C) and 1.1 km (89.7 °C) for x = 1.9 km). This RTM temperature data is available as Supplementary Material. However, due to the low spatial resolution and associated uncertainties, this distance should not be overinterpreted, but serve as an overall estimate. The extensive reservoir cooling (for flow rates >30 L/s) predicted by the RTM here is similar to those found in (Gong et al., 2019) who observed cooling at the injection well after a 50 years model simulation in a low porosity, fractured reservoir.

Figure 4
Three heat maps depict temperature variations with fluid flow rates of 10, 20, and 30 liters per second. Each map shows temperature gradients from yellow to red, with arrows indicating flow direction and magnitude. The top map represents 10 liters, the middle 20 liters, and the bottom 30 liters per second. All maps cover a similar spatial grid with temperature ranges from 40 to 95 degrees Celsius.

Figure 4. Effect of TS reservoir temperature as a function of injection rate after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

3.2 Mineral precipitation evolution

After a simulation time period of 50 years, X2t results indicate that the maximum concentration of each mineral in the bulk assemblage consists of quartz (4*104 mol/m3), albite (3*102 mol/m3), illite (2*102 mol/m3), muscovite (2*102 mol/m3), trydimite (0.25 mol/m3), kaolinite (0.2 mol/m3), witherite (8*10–3 mol/m3), and daphnite−14A (2*10–3 mol/m3). Given that quartz, albite, illite, and muscovite occupy the vast majority of the solid phases present, these will be the focus of the ensuing discussion. Tridymite in particular would not be stable at these P-T conditions.

It is undoubtedly true that the precipitation and subsequent dissolution of other secondary minerals will have affected the fluid-rock interaction over the course of the RTM simulation. However, the focus of this work is to evaluate the state of the TS after 50 RTM simulation years and to evaluate the ‘end product’ results at that time. Clay minerals such as kaolinite are well known for their ability to adversely affect sandstone permeability [e.g. (Howard, 1992; Zhou et al., 2022)], and silica (as amorphous silica) is a common scaling solid in geothermal systems (Gallup et al., 2015). Another scaling mineral predicted in the sandstone-based Mezőberény site in Hungary is the mineral goethite [e.g. (Markó et al., 2024)]. As noted by Markó et al. (2024), goethite was predicted in the screened well sections, and was observed in on-site sampling and from filtered water at the surface. One possible explanation for the absence of goethite in the TS simulations is that the Újfalu and Zagyva formations contains both magnetite and chlorite as sources of Fe (the measured Fe concentration at the production well at Mezőberény site was approximately double that of the 100× brine for the TS). One complicating factor to further consider is the redox equilibrium state of both systems. The temperature at Mezőberény is approximately double that of the TS (Markó et al., 2024) which would increase the rate of Fe2+/Fe3+ speciation reactions trending towards equilibrium. A GWB Act2 model calculation was performed at 90.56 C calculating the H2O stability field boundary (such calculations are minimally affected by ionic strength (Bowman et al., 2023) so the results should be reasonable estimate for undilute or dilute brine) and indicates for a constant pH of 5.49 that the upper and lower brine Eh (V) boundaries are 0.87 V and −0.45 V, respectively. Although the true Eh of the brine is unknown at every node point in the RTM simulation, constraints may be found which help elucidate the proportion of reduced to oxidized iron. Thus, at 90.56 as per the Nerst equation, these electric potentials correspond to [Fe2+]/[Fe3+] of 0.04 (favoring Fe3+ abundance) and 7*1016 (favoring Fe2+), respectively. This iron ratio is similar to the 0.049 ratio as measured from recovered 60 geothermal water from a sandstone reservoir in the Decheng District, China (Feng et al., 2025). Despite the oxidizing nature and Fe3+ abundance in the brine, these researchers do not report secondary goethite. It is reasonable that the TS brine is oxidizing and possible that the true Fe3+ concentration is higher than is predicted by the RTM simulation. Still, it is worth mentioning that these calculations are for static closed systems, while true geothermal operations are performed in open systems in 3D where dissolved species are transported from one location to another by advection. Under this consideration, the activity of Fe3+—and by extension, goethite—could remain undersaturated except at specific locations such as the screened wellbore [e.g. (Markó et al., 2024)].

3.2.1 SiO2/quartz

Very little change is observed for quartz irrespective of injection volume flow rate (Figure 5). The lowest amount of quartz is located within the same grid node as the injection well but is otherwise effectively homogenous (or that spatial differences are below detection of the RTM simulation). At 10 L/s injection, quartz is at equilibrium throughout the RTM block with variation in Q/K at the <0.0036 level. This observation is also found for injection rates of 20 L/s and for 30 L/s, as quartz remains practically at equilibrium (0.99 < Q/K < 1) throughout the RTM block. Even though only a single kinetic reaction pathway may be considered by the X2t module (dissolution rate or precipitation rate) at a time for each mineral, the thermodynamic equilibrium prediction is supported by the kinetic results. Secondary quartz precipitation does not appear to occur. The variation in quartz mol/m3 between that at the injection well (40,100 mol/m3) and elsewhere in the RTM block (40,110 mol/m3) is very small—a difference of 0.025%. Given that the TS is predominately comprised of quartz, along with the slow rate of dissolution, it is not surprising that little change occurs. Such observations lend validity to the TS as a prospective geothermal reservoir as the mineralogy is comparatively refractive to change and is relatively simple with the major constituent mineral being quartz prior to and after 50 years of simulation.

Figure 5
Three rows of graphs labeled A, B, and C illustrate quartz concentration at different flow rates: 10, 20, and 30 liters per second. Each row has two plots, one showing quartz concentration in mol/m³ with a red to yellow gradient, and the other displaying the quartz Q/K ratio with a yellow to red gradient. Arrows indicate flow direction, and scales are provided for each graph.

Figure 5. Effect of injection flow rate on quartz mineralization (as mol/m3) and saturation after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

Interestingly, neither chalcedony nor amorphous silica precipitated as secondary minerals. To confirm the absence of both phases, the RTM simulation was re-run for all injection rates with quartz suppression (Supplementary Material). The results again showed that no amorphous silica or chalcedony had precipitated—despite an instantaneous precipitation—and is attributed to its undersaturation and to the decreased solubility of quartz in comparison to amorphous silica at simulation temperatures (Gunnarsson and Arnórsson, 2000). For all spatial nodes of the RTM block, the greatest saturation index for amorphous silica and chalcedony under ideal conditions (quartz suppression applied) is 0.12 and 0.6, respectively. Both instances occur at an injection rate of 10 L/s and decrease progressively as injection rate is increased. At 30 L/s injection rate, the greatest saturation index throughout the RTM block for amorphous silica and chalcedony are 0.09 and 0.57, respectively. Although the dissolution rate for amorphous silica in ‘pure water’ (reasonable for this work because the dilute working fluid dominates the TS formation brine) at 90.56 is 1.32 × 10−10 mol/m2s (Icenhow et al., 2000) is faster than for quartz (6.94 × 10−12 mol/m2s (Schwartzentruber et al., 1987) at 90 in 1 molal NaOH solution), we did not observe amorphous silica as a solid phase. In agreement with these results, RTM simulations performed on a sandstone reservoir containing 25% quartz did not report chalcedony nor amorphous silica precipitation, although K-feldspar and Na-feldspar were reported to have dissolved; note that the dissolution of each mineral generates H4SiO4 (Feng et al., 2025). However, it is reasonable that given the unsaturated state of albite in the TS results, H4SiO4 is not being released as a Na-feldspar dissolution product and therefore the activity of this species is not sufficiently high for SiO2 mineralization.

3.2.2 Albite

Similar to quartz, albite mineralization is insensitive to changes in injection rate (Figure 6). While the spatial maps show a color gradation change, this change, like that of quartz, is below detection. Although albite has the lowest saturation immediately near the injection well site (where the continuously injected fluid is relatively low in Na+ and Al3+), it is still present at 3.05 volume % at all nodes in the RTM block. Although albite is thermodynamically favored to dissolve, the suspected reason that measurable dissolution is not observed is because the prescribed reaction rate is simply too slow for sufficient change over a 50-year period. This limitation is the result of the X2t module which does not allow for transition state theory to be incorporated into kinetic RTM simulations. Instead, only one forward or reverse reaction rate is allowed. RTM simulation performed by Feng et al. (2025) in a sandstone reservoir consisting 13% Na-feldspar observed a volume reduction of 1.474*10–3% after a simulation period of 30 days (Feng et al., 2025). In consideration of the TS RTM results, this does not mean that albite experienced no volume change, rather that the change (if it occurred) was below the numerical detection limit of the software program. The smectite mineral montmorillonite is a byproduct of albite dissolution, and has been reported in experimental water-sandstone experiments at 45 C with Mg2+ concentrations >78.5 mg/L (70). Although the undiluted and 2× diluted brines contain Mg2+ in concentrations much greater than 78.5 mg/L, the injected 100× dilute working fluid (which quickly overcomes the undiluted native formation brine) contains at >5.7× less Mg2+. This discrepancy lends support for the absence of secondary montmorillonite in the GWB X2t simulation.

Figure 6
Three sets of diagrams showing albite concentration and quality at different flow rates: 10, 20, and 30 liters per second. Each set contains two grids; the left grid shows albite concentration in mol/m³ and the right grid shows albite Q/K ratio. The color gradient indicates concentration levels, with arrows representing flow direction.

Figure 6. Effect of injection flow rate on albite mineralization (as mol/m3) and saturation after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

3.2.3 Illite

Illite mineralization is homogenous at 2.88 volume % throughout the RTM block at both 10 L/s and at 20 L/s injection rate (Figure 7). Illite is the third most abundant mineral present in the TS after 50 years, despite having an initial starting volume percentage that is equal to albite and muscovite (i.e., 3% as in Table 1). When the saturation state and the volume % of illite are considered, RT M results show that illite is in disequilibrium and dissolution is favored, although any changes to reservoir volume % are below software quantification. In addition, the degree of saturation decreases in quantity as well as in its spatial extent over time. Like albite, the greatest amount of saturation (albeit still undersaturated) is located away from the injection wellsite where comparatively dilute working fluid is introduced to the TS reservoir. Although smectite clays are not identified either as primary reservoir minerals in the TS [Table 1 (Houseknecht, 1988)], it is not explicitly excluded. In addition, Houseknecht (1988) states that feldspar is “mostly albite” (Houseknecht, 1988), leaving open the potential that potassium feldspar could be present in low abundance. Smectite clays such as montmorillonite are problematic in subsurface reservoirs because of their swelling capability which can adversely affect permeability [e.g. (Zhou et al., 1996)]. At the base of the TS, the matrix is noted for its muddy composition and directly overlies the graywacke and shale-bearing Juniata Formation (Folk, 1960). Therefore, due to possible smectite in this zone it could be rendered less desirable for geothermal development. However, if smectite cannot be avoided, then transformation to illite could mitigate swelling concerns. Illitization occurs where sufficient K+ (RTM simulations of a shale unit of the Da’anzhai Member, China showed that the greatest transformation of smectite to illite occurred for K+ brine concentrations ranging from 3.4425*10–3 mol/L to 1.0318*10–2 mol/L) and primary smectite are available (Yang et al., 2023). The initial K+ concentrations of the undiluted, 2× diluted, and 100× diluted brines fall outside of this K+ concentration zone, although intermediate dilutions ranging from, for example, 10% undiluted and 90% 100× diluted, to 15% undiluted and 85% 100× diluted brine yield favorable K+ concentrations for illitization.

Figure 7
Three panels labeled A, B, and C display data on illite concentration and flow vectors for different flow rates: 10, 20, and 30 liters per second. Each panel has two graphs comparing illite in mol/m³ and illite Q/K, colored from yellow to red indicating concentration levels, with arrows showing flow direction. The concentration is shown over grid coordinates ranging from 0 to 10e4 cm.

Figure 7. Effect of injection flow rate on illite mineralization (as mol/m3) and saturation after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

3.2.4 Muscovite

Across all injection rates, muscovite mineralization is practically homogenous like that for quartz, albite, and illite. There does exist a very small spatial variation in the volume %, with a value of 2.813% in the same node as the injection well, and ∼2.815% throughout the remainder of the block—translating to the 199.9 mol/m3 at the injection well to 200 mol/m3 elsewhere (Figure 8). Despite the homogeneity in mineral volume concentration, there is a strong spatial relationship in equilibrium saturation. At an injection rate of 10 L/s muscovite is undersaturated immediately near the injection well; a region which continues to increase in spatial area commensurate with the injection flow rate. Despite this, muscovite is supersaturated in a radial band at a distance of approximately 650 m from the well site. This supersaturation band is not perfectly symmetrical (taking on an oblate shape) with the magnitude of supersaturation decreasing along the perimeters of the band. This radial band is also observed at 20 L/s, but not at 30 L/s, where the corresponding ion activities have presumably been transported out of the region of interest of the RTM block. There does appear to be a remnant present in the top right of the RTM block (Figure 6C). The RTM results appear to show a spatial relationship between temperature-muscovite saturation and pH-muscovite saturation which are in general agreement with the sandstone-based experimental results noted by Gan et al. (2022). RTM grid points containing saturated to supersaturated muscovite are also in close proximity to points which are on the migration front of the cooler injected working fluid. At an injection rate of 30 L/s the supersaturated muscovite grid points have been transported outside the bounds of the RTM block save for one in the upper right corner (probably a trailing limb of the migrating plume). In addition, grid points in proximity to the injection site where injected fluid is acidic, muscovite is the strongly undersaturated. In contrast, increased muscovite saturation tends to follow the leading edge of the migrating plume which also has an increased pH.

Figure 8
Three pairs of contour plots labeled A, B, and C depict muscovite concentration changes at different flow rates: 10, 20, and 30 liters per second. Each pair consists of a plot for molarity (left) and muscovite reaction quotient (Q/K, right). Overlapping arrow vectors indicate flow direction and magnitude. The molarity plots use a red-yellow scale, and the Q/K plots use an orange-yellow scale, both indicating higher concentrations in yellow areas.

Figure 8. Effect of injection flow rate on muscovite mineralization (as mol/m3) and saturation after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

3.3 pH evolution

Throughout the RTM block and regardless of injection flow rate, the pH of the aqueous fluid in the TS is acidic (Figure 9). At 10 L/s, circumneutral pH is found in an oblate radial zone emanating to a distance of approximately 550 m (pH = 6.6), while at 20 L/s the pH at 550 m is ∼6.9. At 30 L/s and a distance of 550 from the injection well, the pH is ∼6.6. The thickness of this radial band increases outward in the direction of fluid flow in proportion to elapsed time (see Supplementary Material). It is important to note that while this radial region has the highest pH of the RTM block, the pH immediately adjacent to the injection well is the lowest at 5.9. This proportional relationship between brine pH and distance from injection site has been observed by Ma et al. (2023) in another RTM study where acidic fluid was injected at a rate of 30 kg/s (similar to the upper injection rate used here) (Ma et al., 2023). In their study, the RTM formation media was fractured, and the spatio-temporal pH change occurred much more rapidly than in the unfractured TS example. If the TS were fractured either hydraulically or through liquid nitrogen (Wang et al., 2025), evolution in the brine pH should be more rapid. If flow of the working fluid is dominated through the induced fracture network, pH-sensitive dissolution or precipitation reactions may be more localized in proximity to the fracture.

Figure 9
Three contour diagrams labeled A, B, and C show pH variations at depths over time for water flow rates of 10, 20, and 30 liters per second. Each grid uses a color gradient from red to yellow representing pH levels from 5.4 to 7. Arrows indicate flow direction and magnitude. The background color intensity increases with higher flow rates, suggesting changing pH distribution and flow dynamics with increased water movement.

Figure 9. Effect of injection flow rate on aqueous fluid pH after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

3.4 Porosity and permeability evolution

After a simulation time of 50 years, reservoir porosity is constant (Figure 10). Across all injection rates, spatial changes in porosity values are minimal (0.01% at 10 L/s, 0.02% at 20 L/s, and 0.03% at 30 L/s). In addition, the magnitude of the range of porosities relative to the initial starting value of 8.1% is also negligible as shown in Table 5. These results show that although both intra-and inter-porosity ranges are small across all injection rates, the rate of change in porosity increases with injection rate. In other words, porosity (localized changes both decrease and increase) is altered more significantly as the injection rate is increased. Although these small differences have little value in this context, changes to matrix porosity could be statistically significant at much higher injection rates.

Figure 10
Three sets of side-by-side contour plots labeled A, B, and C show porosity and x-permeability at different discharge rates: 10, 20, and 30 liters per second. Each set features a left plot for porosity and a right plot for x-permeability, both with arrow vectors indicating flow direction. The porosity scale ranges from yellow to red, while the x-permeability scale is in shades of red. Measurements are marked for year fifty.

Figure 10. Effect of injection flow rate on porosity (panels in left column) and permeability (panels in right column) after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows. (A) 10 L/s. (B) 20 L/s. (C) 30 L/s.

Table 5
www.frontiersin.org

Table 5. Changes in matrix TS reservoir porosity at varied injection flow rates.

The effect of injection flow rate on permeability is nearly zero (Figure 10). Observed permeabilites of ∼230 mD are well within the stated range of the TS (Section 2.1.1). Even at the injection site, the observed change is apparent, as the value here and throughout the rest of the RTM block is 229.5 mD. Deposition (where mol/m3 > 0) of pore-clogging minerals such as kaolinite (Howard, 1992; Zhou et al., 2022) are in too low of volume concentration to contribute to adverse changes in permeability (Figure 11). This does not mean that these minerals could not contribute to localized reductions in permeability, but rather that any highly localized adverse permeability changes are smoothed out by the large node. Experimental core-flooding results obtained in sandstone performed by Chai et al. (2022) showed that mineral dissolution reduced pore throat size, pore connectivity, and formation mechanical strength leading to decreased water transmissivity (Chai et al., 2022). In comparison to the often-encountered discussion where high salinity leads to secondary mineralization, these researchers interestingly found that adverse water transmissivity was due to a reduction in salinity. Perhaps then, from a geochemical standpoint alone, geothermal reservoirs consisting of a high abundance of quartz (such as the TS) are advantageous relative to feldspar or clay-rich sandstones, as they; (1) experience a simple brine chemistry leading to somewhat predictable secondary mineralization; (2) contain less reactive minerals leading to slower or less dissolution thereby retaining formation structural integrity. Identifying the ideal target brine salinity remains one of the larger unknown concerns not only for the TS, but presumably for other geothermal reservoirs. A high salinity brine favors precipitation of scale-forming minerals which clog pore spaces, however, a low salinity brine has been shown (Chai et al., 2022) to result in formation damage which also reduces water throughput.

Figure 11
Grid chart showing changes in kaolinite, with green squares indicating kaolinite dissolution and pink squares indicating precipitation over time. Arrows represent directional flow, and a color gradient at the bottom shows X-permeability in millidarcies for year 50, ranging from yellow (low) to red (high).

Figure 11. Overlay of the permeability spatial map from Figure 8B at 20 L/s and kaolinite precipitation (mol/m3 is positive)-“dissolution” (kaolinite is not present in these locations, change after 50 years of continual flow. Flow enters from the injection well and is denoted by the crosshairs. Fluid transport direction and magnitude are indicated by arrows.

There is a negative correlation between sandstone pore complexity and its permeability [e.g. (Teng et al., 2025)]. Although the 6 MPa triaxial pressure and mineralogy (54% potassium feldspar and 17% quartz) utilized by these workers are different from the parameters of the TS at the study site (Figure 2), their observation nonetheless provides an interesting insight regarding fluid flow; the high pressure and low permeability of the TS in north central West Virginia suggests that the pore distribution is likely complex and heterogeneous with flow controlled primarily through mesoscopic and macroscopic connected pore space. It is important to emphasize that the mineralogical reservoir matrix in Teng et al. (2025) varies from the TS and that primary mineral reactivity may exert control on fluid flow unequally as one mineral preferentially dissolves over another. For the TS, pore heterogeneity might be further increased with flow favoring dissolution-induced channels. Although the 50-year simulation results shown here do not predict secondary mineralization of swelling clays (illite and kaolinite are non-swelling Farrokhpay et al., 2016), these clays could be present in other geothermal formations elsewhere and could contribute to volume reduction and exacerbate selective channeling. Such insight to other prospective geothermal reservoirs and reservoir engineering and management.

3.5 Model limitations

This modeling approach does not consider a two well (injection and production well) model, which would be present in a real-world operation. The second well was excluded in this model as it would present significant complications and many assumptions. For example, the propensity and direction of induced channel connectivity generated by hydraulic fracturing and the production well flow rate; which can vary enormously, presents a large unknown (Kamila et al., 2021). In addition, the distance between the production well and the injection well looks to optimize production output, which is outside the scope of this study (Lu, 2018). Another limitation of the model is the assumption of well efficiency (i.e., the percentage of injected fluid that is returned to the surface) related to the pore connectivity. As the TS is characterized by low permeability, hydraulic fracture stimulation would be expected for true operation. Since fracture identification from resistivity logs is prone to human error or imprecision (Fathi et al., 2022) thereby introducing inaccurate RTM results, this model considers unfractured rock. To this end, we have opted for a single well to model changes in the TS reservoir over time and distance.

The utility of activity models built into RTM, such as GWB, makes simulations of brines difficult. The actual TS reservoir consists of a 4+ mol/kg solution; ideally, the RTM would incorporate the raw, undiluted brine. We were able to simulate a 2× dilution for the brine, but the results from activity models are not suited for high ionic strength fluids Still, manmade reservoirs are operated using highly diluted brine indicating that 2× diluted brine is unexpected to present significant inaccuracies. The X2t module does not allow for both a forward dissolution rate and a reverse precipitation rate to be considered. Instead, only a single reaction pathway rate may be prescribed. In this modeling work, we have opted to include kinetic dissolution rates for the TS host rock minerals. Secondary precipitation, should it occur, is considered in equilibrium with the evolving brine phase. As GWB utilizes the Law of Mass Action, the formation of solid solutions is not considered.

4 Conclusion

• From a geochemical perspective, the TS appears to be well-suited for geothermal potential. RTM simulation results suggest that this geologic unit experiences very minor chemical change over a typical operational lifespan of 50 years. The simple quartz-dominant mineralogy is retained while secondary minerals are present only at very low abundances (<<< 1 volume %). For quartz-dominated reservoirs such as the TS which experience the injection of a dilute working fluid, geochemical changes are less important in comparison to reservoir cooling.

• RTM simulation results show that the rate at which injected fluid is introduced plays a significant role in the reservoir temperature and heat distribution. Injection flow rates above 10 L/s are likely to substantially lower the reservoir temperature even at distances greater than 1,000 m from the injection site. The temperature threshold required to sustain the requisite work at the surface (e.g., electricity production, direct heating, etc.) is beyond the scope of this work. However, it is an important consideration that should be made in conjunction with other parameters presented here.

• Neither flow nor kaolinite appears to affect permeability. Although there is considerable variation in the initial TS permeability, hydraulic fracturing might not be necessary, perhaps as preferential dissolution fissures or fluid channels are established in situ [e.g. (Worthington and Ford, 2009)]. In north-central WV, the high-pressure and low permeability TS, could indicate pore heterogeneity favoring dissolution-based channels which control fluid flow.

• For fluid transported through matrix porosity alone and without any secondary porosity induced by hydraulic fracturing, reservoir temperature is highly sensitive to injection rate. At 10 L/s, reservoir temperature decreases to ∼75 °C at a distance of 750 m from the injection well. As the injection flow rate is increased, the cooling is more profound and widespread. The presence of successful fluid transmissivity without secondary porosity, as modeled, has implications for mitigating thermal drawdown-induced flow channeling–a common challenge towards regulating thermal longevity in reservoirs. Optimization studies on working fluid injection temperature, production well placement, and operation timescales may lead to an increased capacity for maintaining reservoir temperatures. Thermal breakthrough occurs or may be constrained as follows:

a. 10 L/s injection rate, >750 m after 50 years

b. 20 L/s injection rate, 950 m after 30 years

c. 30 L/s injection rate, reservoir is effectively entirely at 1 km, but after 50 years still retains a temperature of 89.5 °C at 1.9 km distance. Geothermal systems that necessitate a higher injection rate into reservoirs such as the TS would benefit by placing the return/production well farther away from the injection well. This study reiterates the importance of performing elementary model assessments of prospective reservoirs prior to field implementation.

• Identifying suitable production well placement using data from a matrix-only porosity RTM is difficult. Improved well location identification could be established if future TS RTM simulations incorporated secondary porosity.

• Since significant quantities of dilute exogenous water are required, adverse scaling is expected to be significantly minimized compared to natural hydrothermal systems. The results here indicate that both precipitation and dissolution in a quartz arenite-based reservoir are minimal with minor alteration to the framework mineralogy. Indeed, this simple mineralogy suggests that if the TS were selected for geothermal implementation, scale prevention and treatment should also be straightforward.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

SB: Conceptualization, Methodology, Formal analysis, Investigation, Visualization, Project administration, Writing – original draft, Writing – review and Editing. JB: Investigation, Writing – original draft, Writing – review and editing. AP: Investigation, Writing – original draft, Writing – review and editing. SS: Funding acquisition, Resources, Writing – original draft, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. We acknowledge the Department of Energy Office of Energy Efficiency and Renewable Energy (DOE EERE) Geothermal Technologies Program Award # DE-EE0009597 to SS.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2025.1706071/full#supplementary-material

Footnotes

1Personal communication, DOE report 2023

References

Abarca, E., Idiart, A., Grabdia, F., Rodríguez-Morillas, N., Pellan, C., Zen, M., et al. (2019). 3D reactive transport modeling of porosity evolution in a Carbonate reservoir through dolomitization. Chem. Geol. 513, 184–199. doi:10.1016/j.chemgeo.2019.03.017

CrossRef Full Text | Google Scholar

Augustine, C. (2016). Update to enhanced geothermal system resource potential estimate. GRC Trans. 40, 673–677. Available online at: https://www.geothermal-library.org/index.php?mode=pubs&action=view&record=1032382.

Google Scholar

Bebout, D. G., Loucks, R. G., and Gregory, A. R. (1978). Frio sandstone reservoirs in the deep subsurface along the Texas Gulf Coast - their potential for production of geopressured geothermal energy; report of investigations no. 91. Austin, TX: Bureau of Economic Geology.

Google Scholar

Benjakul, R., Hollis, C., Robertson, H. A., Sonnenthal, E. L., and Whitaker, F. F. (2020). Understanding controls on hydrothermal dolomitisation: insights from 3D reactive transport modelling of geothermal convection. Solid Earth 11 (6), 2439–2461. doi:10.5194/se-11-2439-2020

CrossRef Full Text | Google Scholar

Bethke, C. M. (2008). Geochemical and biogeochemical reaction modeling. 3rd ed. Cambridge University Press.

Google Scholar

Bethke, C. M., Farrell, B., and Sharifi, M. (2024). GWB reactive transport modeling guide; the geochemist’s workbench; release 17; aqueous solutions. Champaign: LLC, 191.

Google Scholar

Blackwell, D., Richards, M., Frone, Z., Ruzo, A., Dingwall, R., and Williams, M. (2011). Temperature-at-depth maps for the conterminous US and geothermal resource estimates. GRC Trans. 35, 1545–1550. Available online at: https://www.geothermal-library.org/index.php?mode=pubs&action=view&record=1029452.

Google Scholar

Borhara, K., and Onasch, C. M. (2020). Evidence for silica gel and its role in faulting in the Tuscarora sandstone. J. Struct. Geol. 139, 104140. doi:10.1016/j.jsg.2020.104140

CrossRef Full Text | Google Scholar

Bowman, S., Pathak, A., Agrawal, V., and Sharma, S. (2023). Effect of ionic strength on H2O and Si-species stability field geometry in pH-Eh space. Aquat. Geochem. 29 (4), 207–218. doi:10.1007/s10498-023-09417-0

CrossRef Full Text | Google Scholar

Bryndzia, L. T., and Fay, M. C. (2016). Geochemical analysis of returned treatment waters (RTW) associated with shale gas production in the appalachian basin (USA) and deep basin (Canada): potential use of total dissolved solids (TDS) and oxygen isotope data forssessing water:rock ratios and stimulated rock volume (SRV). OnePetro. doi:10.15530/URTEC-2016-2455905

CrossRef Full Text | Google Scholar

Camp, E. R., Jordan, T. E., Hornbach, M. J., and Whealton, C. A. (2018). A probabilistic application of oil and gas data for exploration stage geothermal reservoir assessment in the appalachian Basin. Geothermics 71, 187–199. doi:10.1016/j.geothermics.2017.09.001

CrossRef Full Text | Google Scholar

Campbell, S. A., Lenhardt, N., Dippenaar, M. A., and Götz, A. E. (2016). Geothermal energy from the main Karoo Basin (South Africa): an outcrop analogue study of Permian sandstone reservoir formations. Energy Procedia 97, 186–193. doi:10.1016/j.egypro.2016.10.050

CrossRef Full Text | Google Scholar

Castle, J. W., and Byrnes, A. P. (2005). Petrophysics of Lower Silurian sandstones and integration with the tectonic-stratigraphic framework, Appalachian Basin, United States. AAPG Bull. 89, 41–60. doi:10.1306/08170404028

CrossRef Full Text | Google Scholar

Chai, R., Liu, Y., Xue, L., Rui, Z., Zhao, R., and Wang, J. (2022). Formation damage of sandstone geothermal reservoirs: during decreased salinity water injection. Appl. Energy 322, 119465. doi:10.1016/j.apenergy.2022.119465

CrossRef Full Text | Google Scholar

Chapman, E. C., Capo, R. C., Stewart, B. W., Hedin, R. S., Weaver, T. J., and Edenborn, H. M. (2013). Strontium isotope quantification of siderite, brine and acid mine drainage contributions to abandoned gas well discharges in the Appalachian Plateau. Appl. Geochem. 31, 109–118. doi:10.1016/j.apgeochem.2012.12.011

CrossRef Full Text | Google Scholar

Donselaar, M. E., Groenenberg, R. M., and Gilding, D. T. (2015). Reservoir geology and geothermal potential of the delft sandstone member in the West Netherlands Basin; Melbourne.

Google Scholar

Dresel, P. E., and Rose, A. W. (2010). Chemistry and origin of oil and gas well brines in Western Pennsylvania. Pa. Dep. Conserv. Nat. Resour. Open-File Report OFOG 10-10.0.

Google Scholar

Erol, S., Akın, T., Başer, A., Saraçoğlu, Ö., and Akın, S. (2022). Fluid-CO2 injection impact in a geothermal reservoir: evaluation with 3-D reactive transport modeling. Geothermics 98, 102271. doi:10.1016/j.geothermics.2021.102271

CrossRef Full Text | Google Scholar

Farrokhpay, S., Ndlovu, B., and Bradshaw, D. (2016). Behaviour of swelling clays versus non-swelling clays in flotation. Min. Eng. 96–97, 59–66. doi:10.1016/j.mineng.2016.04.011

CrossRef Full Text | Google Scholar

Fathi, E., Carr, T. R., Adenan, M. F., Panetta, B., Kumar, A., and Carney, B. J. (2022). High-quality fracture network mapping using high frequency Logging While Drilling (LWD) data: MSEEL case study. Mach. Learn. Appl. 10, 100421. doi:10.1016/j.mlwa.2022.100421

CrossRef Full Text | Google Scholar

Feng, B., Yang, J., Zhao, J., Yang, Y., Tian, H., Feng, G., et al. (2025). Study on porosity and permeability characteristics of sandstone geothermal reservoir under recharge conditions: a case study of Decheng District, Shandong Province. Energies 18 (22), 6060. doi:10.3390/en18226060

CrossRef Full Text | Google Scholar

Folk, R. L. (1960). Petrography and origin of the Tuscarora, Rose Hill, and Keefer formations, lower and Middle Silurian of Eastern West Virginia. J. Sediment. Petrol. 30, 1–58. doi:10.1306/74D709C5-2B21-11D7-8648000102C1865D

CrossRef Full Text | Google Scholar

Gallup, D. L., and von Hirtz, P. (2015). “Chapter 22 - control of Silica-Based scales in cooling and geothermal systems,” in Mineral scales and deposits. Editors Z. Amjad, and K. D. Demadis (Amsterdam: Elsevier), 573–582. doi:10.1016/B978-0-444-63228-9.00022-X

CrossRef Full Text | Google Scholar

Gan, H., Liu, Z., Wang, X., Zhang, Y., Liao, Y., Zhao, G., et al. (2022). Effect of temperature and acidification on reinjection of geothermal water into sandstone geothermal reservoirs: laboratory study. Water 14 (19), 2955. doi:10.3390/w14192955

CrossRef Full Text | Google Scholar

Garapati, N., and Hause, J. (2021). Feasibility of deep direct-use geothermal on the West Virginia University Campus-Morgantown, WV; final technical report. DOE EERE-Geothermal Technologies Program, 143. doi:10.2172/1829981

CrossRef Full Text | Google Scholar

Garapati, N., Alonge, O., Hall, L., Irr, V., Zhang, Y., Smith, J., et al. (2019). “Feasibility of development of geothermal deep direct-use district heating and cooling System at West Virginia University Campus-Morgantown, WV,” in Proceedings. Stanford: SGP-TR-214.

Google Scholar

Gong, F., Guo, T., Sun, W., Li, Z., Yang, B., Chen, Y., et al. (2019). Evaluation of geothermal energy extraction in Enhanced Geothermal System (EGS) with Multiple Fracturing Horizontal Wells (MFHW). Renew. Energy 151, 1339–1351. doi:10.1016/j.renene.2019.11.134

CrossRef Full Text | Google Scholar

Gunnarsson, I., and Arnórsson, S. (2000). Amorphous silica solubility and the thermodynamic properties of H4SiO4 in the range of 0° to 350°C at psat. Geochim. Cosmochim. Acta 64 (13), 2295–2307. doi:10.1016/S0016-7037(99)00426-3

CrossRef Full Text | Google Scholar

Harris, K. R., and Woolf, L. A. (1980). Pressure and temperature dependence of the self diffusion coefficient of water and oxygen-18 water - journal of the Chemical Society, faraday transactions 1: physical chemistry in condensed phases (RSC publishing). J. Chem. Soc. Faraday Trans. 76, 377–385. doi:10.1039/f19807600377

CrossRef Full Text | Google Scholar

Heald, M. T., and Anderegg, R. C. (1960). Differential cementation in the Tuscarora Sandstone [Virginia-West Virginia]. J. Sediment. Res. 30 (4), 568–577. doi:10.1306/74D70AA1-2B21-11D7-8648000102C1865D

CrossRef Full Text | Google Scholar

Holmslykke, H. D., and Kjøller, C. (2023). Reactive transport modelling of potential near-well mineralogical changes during seasonal heat storage (HT-ATES) in Danish geothermal reservoirs. J. Energy Storage 72, 108653. doi:10.1016/j.est.2023.108653

CrossRef Full Text | Google Scholar

Houseknecht, D. W. (1988). Intergranular pressure solution in four quartzose sandstones. J. Sediment. Res. 58 (2), 228–246. doi:10.1306/212F8D64-2B24-11D7-8648000102C1865D

CrossRef Full Text | Google Scholar

Howard, J. J. (1992). Influence of authigenic-clay minerals on permeability, origin, diagenesis, and petrophysics of clay minerals in sandstones. Tulsa, OK: SEPM Society for Sedimentary Geology 47, 257–264.

CrossRef Full Text | Google Scholar

Icenhower, J. P., and Dove, P. M. (2000). The dissolution kinetics of amorphous silica into sodium chloride solutions: effects of temperature and ionic strength. Geochim. Cosmochim. Acta 64 (24), 4193–4203. doi:10.1016/S0016-7037(00)00487-7

CrossRef Full Text | Google Scholar

Kamila, Z., Kaya, E., and Zarrouk, S. J. (2021). Reinjection in geothermal fields: an updated worldwide review 2020. Geothermics 89, 101970. doi:10.1016/j.geothermics.2020.101970

CrossRef Full Text | Google Scholar

Kaya, E., Zarrouk, S. J., and O’Sullivan, M. J. (2011). Reinjection in geothermal fields: a review of worldwide experience. Renew. Sustain. Energy Rev. 15 (1), 47–68. doi:10.1016/j.rser.2010.07.032

CrossRef Full Text | Google Scholar

Khutorskoi, M. D., and Polyak, B. G. (2016). Role of radiogenic heat generation in surface heat flow formation. Geotectonics 50 (2), 179–195. doi:10.1134/S0016852116020047

CrossRef Full Text | Google Scholar

Klunk, M. A., Dasgupta, S., Das, M., Conceição, R. V., Siqueira Xavier, S. J., Chemale, F., et al. (2021). Application of geochemical modelling software as a tool to predict the diagenetic reactions between the marine connate water and the salt dome in a petroleum system. J. South Am. Earth Sci. 109, 103272. doi:10.1016/j.jsames.2021.103272

CrossRef Full Text | Google Scholar

Kramer, C. T. (2013). Regional stratigraphic framework of the Tuscarora sandstone: a model for geologic CO2 storage in West Virginia. MS. Morgantown, WV: West Virginia University. doi:10.33915/etd.331

CrossRef Full Text | Google Scholar

Krynicki, K., Green, C. D., and Sawyer, D. W. (1978). Pressure and temperature dependence of self-diffusion in water. Faraday Discuss. Chem. Soc. 66, 199. doi:10.1039/dc9786600199

CrossRef Full Text | Google Scholar

Lichtner, P. C. (1996). “Continuum formulation of multicomponent-multiphase reactive transport,” in Reactive transport in porous media. Reviews in mineralogy. Editors P. C. Lichtner, C. I. Steefel, and E. H. Oelkers (Washington, DC: Mineralogical Society of America), 34, 1–81.

Google Scholar

Lowry, T. S., Peplinski, W., and Martinez, J. (2018). Examining the tradeoff between sustainability and profitability for enhanced geothermal systems. GRC Trans. 42, 42. Available online at: https://www.geothermal-library.org/index.php?mode=pubs&action=view&record=1034059.

Google Scholar

Lu, S.-M. (2018). A global review of enhanced geothermal system (EGS). Renew. Sustain. Energy Rev. 81, 2902–2921. doi:10.1016/j.rser.2017.06.097

CrossRef Full Text | Google Scholar

Ma, L., Cui, Z., Feng, B., Qi, X., Zhao, Y., and Zhang, C. (2023). Reactive transport modeling of chemical stimulation processes for an Enhanced Geothermal System (EGS). Energies 16 (17), 6229. doi:10.3390/en16176229

CrossRef Full Text | Google Scholar

Markó, Á., Brehme, M., Pedretti, D., Zimmermann, G., and Huenges, E. (2024). Controls of low injectivity caused by interaction of reservoir and clogging processes in a sedimentary geothermal aquifer (Mezőberény, Hungary). Geotherm. Energy 12 (1), 40. doi:10.1186/s40517-024-00317-2

CrossRef Full Text | Google Scholar

McCleery, R. S., McDowell, R. R., Moore, J. P., Garapati, N., Carr, T., and Anderson, B. J. (2018). Development of 3-D geological model of Tuscarora sandstone for feasibility of deep direct-use geothermal at West Virginia university’s main campus, 42.

Google Scholar

McDowell, R. R., Lewis, J. E., Daft, G. W., Dinterman, P. A., Brown, S. R., and Moore, J. P. (2020). Silurian Tuscarora Sandstone in Western West Virginia: will it work as a Geothermal Reservoir Rock? #80726 (2020). AAPG 2019 East. Sect. Meet. Energy Heartl. 12-16 Oct. Columb. 2020.

Google Scholar

Merkel, B. J., and Planer-Friedrich, B. (2005). in Groundwater chemistry: a practical guide to modeling of natural and contaminated aquatic systems. Editor D. K. Nordstrom (Berlin: Springer), 200.

Google Scholar

Moeck, I. S. (2014). Catalog of geothermal play types based on geologic controls. Renew. Sustain. Energy Rev. 37, 867–882. doi:10.1016/j.rser.2014.05.032

CrossRef Full Text | Google Scholar

Murphy, H., Brown, D., Jung, R., Matsunaga, I., and Parker, R. (1999). Hydraulics and well testing of engineered geothermal reservoirs. Geothermics 28 (4), 491–506. doi:10.1016/S0375-6505(99)00025-5

CrossRef Full Text | Google Scholar

Northeast Natural Energy (2023). Pers. Comm.

Google Scholar

Palandri, J., and Kharaka, Y. (2004). Open File Report 2004-1068. A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling. 71. doi:10.3133/ofr20041068

CrossRef Full Text | Google Scholar

Parr. Safety in the operation (2012). Parr. Safety in the operation of laboratory reactors and pressure vessels. Moline: Parr Instrument Company. Available online at: https://www.parrinst.com/wp-content/uploads/downloads/2012/05/230M_Parr_Safety-Lab-Reactors.pdf (Accessed January 02, 2026).

Google Scholar

Pickens, J. F., and Grisak, G. E. (1981). Scale-Dependent dispersion in a stratified granular aquifer. Water Resour. Res. 17 (4), 1191–1211. doi:10.1029/WR017i004p01191

CrossRef Full Text | Google Scholar

Pilewski, J., Sharma, S., Agrawal, V., Hakala, J. A., and Stuckman, M. Y. (2019). Effect of maturity and mineralogy on fluid-rock reactions in the marcellus shale. Environ. Sci. Process. Impacts 21, 845–855. doi:10.1039/c8em00452h

PubMed Abstract | CrossRef Full Text | Google Scholar

Robie, R. A., and Hemingway, B. S. (1995). Thermodynamic properties of minerals and related substances at 298.15 K and 1 Bar (105 Pascals) pressure and at higher temperatures. Washington: United States Government Printing Office. doi:10.3133/sim3067

CrossRef Full Text | Google Scholar

Ryder, R. T. (2004). Stratigraphic framework and depositional sequences in the Lower Silurian regional oil and gas accumulation, appalachian Basin: from Ashland County, Ohio, through Southwestern Pennsylvania, to Preston County, West Virginia; geologic investigations series map; pamphlet I–2810. Reston: U. S. Geological Survey, 11.

Google Scholar

Ryder, R. T., Crangle, R. D., Jr., Trippi, M. H., Swezey, C. S., Lentz, E. E., Rowan, E. L., et al. (2009). Usgs SIM-3067: geologic cross section D–D’ through the Appalachian Basin from the Findlay Arch, Sandusky County, Ohio, to the Valley and Ridge Province, Hardy County, West Virginia; scientific investigations map 3067. Reston: U. S. Geological Survey. doi:10.3133/sim3067

CrossRef Full Text | Google Scholar

Sander, R. (2015). Compilation of Henry’s Law Constants (Version 4.0) for Water as solvent. Atmos. Chem. Phys. 15 (8), 4399–4981. doi:10.5194/acp-15-4399-2015

CrossRef Full Text | Google Scholar

Schott, J., Pokrovsky, O. S., and Oelkers, E. H. (2009). The link between mineral dissolution/precipitation kinetics and solution chemistry. Rev. Mineral. Geochem. 70, 207–258. doi:10.2138/rmg.2009.70.6

CrossRef Full Text | Google Scholar

Schwartzentruber, J., Fürst, W., and Renon, H. (1987). Dissolution of quartz into dilute alkaline solutions at 90°C: a kinetic study. Geochim. Cosmochim. Acta 51 (7), 1867–1874. doi:10.1016/0016-7037(87)90177-3

CrossRef Full Text | Google Scholar

Shevalier, M., Dalkhaa, C., Humez, P., Mayer, B., Becker, V., Nightingale, M., et al. (2014). Coupling of TOUGHREACT-geochemist workbench (GWB) for modeling changes in the isotopic Composition of CO2 leaking from a CCS storage reservoir. Energy Procedia 63, 3751–3760. doi:10.1016/j.egypro.2014.11.404

CrossRef Full Text | Google Scholar

SteefelC. I., (2009). CrunchFlow: software for modeling multicomponent reactive flow and transport; user’s manual. Berkeley: Lawrence Berkeley National Laboratory, 91.

Google Scholar

Sun, Y., Zhang, X., Liu, Y., Zhang, X., Li, X., Sun, C., et al. (2023). Prediction and optimization of productivity and lifespan in multi-well enhanced geothermal system. Appl. Therm. Eng. 234, 121155. doi:10.1016/j.applthermaleng.2023.121155

CrossRef Full Text | Google Scholar

Teng, T., Chen, Y., Wang, Y., and Qiao, X. (2025). In situ nuclear magnetic resonance observation of pore fractures and permeability evolution in rock and coal under triaxial compression. J. Energy Eng. 151 (4), 04025036. doi:10.1061/JLEED9.EYENG-6054

CrossRef Full Text | Google Scholar

Tester, J., W., Beckers, K.; J, F., Hawkins, A., Lukawski, Z., and The, M. (2021). Evolving role of geothermal energy for decarbonizing the United States. Energy Environ. Sci. 14 (12), 6211–6241. doi:10.1039/D1EE02309H

CrossRef Full Text | Google Scholar

Wang, L., Zhu, L., Cao, Z., Liu, J., Xue, Y., Wang, P., et al. (2025). Thermo-mechanical degradation and fracture evolution in low-permeability coal subjected to cyclic heating–cryogenic cooling. Phys. Fluids 37 (8), 086617. doi:10.1063/5.0282266

CrossRef Full Text | Google Scholar

Wanner, C., Peiffer, L., Sonnenthal, E., Spycher, N., Iovenitti, J., and Kennedy, B. M. (2014). Reactive transport modeling of the dixie Valley geothermal area: insights on flow and geothermometry. Geothermics 51, 130–141. doi:10.1016/j.geothermics.2013.12.003

CrossRef Full Text | Google Scholar

Weibel, R., Olivarius, M., Vosgerau, H., Mathiesen, A., Kristensen, L., Nielsen, C. M., et al. (2020). Overview of potential geothermal reservoirs in Denmark. Neth. J. Geosci. 99, e3. doi:10.1017/njg.2020.5

CrossRef Full Text | Google Scholar

Worthington, S. R. H., and Ford, D. C. (2009). Self-organized permeability in carbonate aquifers. Groundwater 47 (3), 326–336. doi:10.1111/j.1745-6584.2009.00551.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, T., Apps, J. A., and Pruess, K. (2005). Mineral sequestration of carbon dioxide in a sandstone–shale system. Chem. Geol. 217 (3), 295–318. doi:10.1016/j.chemgeo.2004.12.015

CrossRef Full Text | Google Scholar

Yang, L., Lu, L., Li, X., Shan, Y., Mo, C., Sun, M., et al. (2023). Clay mineral transformation mechanism modelling of shale reservoir in Da’anzhai member, Sichuan Basin, Southern China. Front. Earth Sci. 11, 1205849. doi:10.3389/feart.2023.1205849

CrossRef Full Text | Google Scholar

Yeakel, L. S., Jr (1962). Tuscarora, juniata, and bald eagle paleocurrents and paleogeography in the central appalachians. GSA Bull. 73 (12), 1515–1539. doi:10.1130/0016-7606(1962)73[1515:tjabep]2.0.co;2

CrossRef Full Text | Google Scholar

Youse, A. C. (1970). “A summary of Tuscarora Sandstone petroleum exploration in West Virginia,” in Silurian stratigraphy, central appalachian Basin-1970 (Charleston: Appalachian Geological Survey), 86–98.

Google Scholar

Yuan, G., Cao, Y., Gluyas, J., and Jia, Z. (2017). Reactive transport modeling of coupled feldspar dissolution and secondary mineral precipitation and its implication for diagenetic interaction in sandstones. Geochim. Cosmochim. Acta 207, 232–255. doi:10.1016/j.gca.2017.03.022

CrossRef Full Text | Google Scholar

Zhou, Z., Gunter, W. D., Kadatz, B., and Cameron, S. (1996). Effect of clay swelling on reservoir quality. J. Can. Pet. Technol. 35 (07), 18–23. doi:10.2118/96-07-02

CrossRef Full Text | Google Scholar

Zhou, Y., Yang, W., and Yin, D. (2022). Experimental investigation on reservoir damage caused by clay minerals after water injection in low permeability sandstone reservoirs. J. Pet. Explor. Prod. Technol. 12, 915–924. doi:10.1007/s13202-021-01356-2

CrossRef Full Text | Google Scholar

Keywords: geothermal reservoir engineering, matrix transport, reactive transport modeling, secondary precipitation, Tuscarora Sandstone

Citation: Bowman S, Barton J, Pathak A and Sharma S (2026) Physicochemical kinetic reactive transport modeling in quartz arenite: implications on geothermal potential. Front. Earth Sci. 13:1706071. doi: 10.3389/feart.2025.1706071

Received: 15 September 2025; Accepted: 11 December 2025;
Published: 12 January 2026.

Edited by:

Giuseppe Oliveto, University of Basilicata, Italy

Reviewed by:

Zhengzheng Cao, Henan Polytechnic University, China
Daniele Pedretti, University of Milan, Italy

Copyright © 2026 Bowman, Barton, Pathak and Sharma. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Samuel Bowman, c2FtaGJvd21hbkBnbWFpbC5jb20=

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.