Hydrologic windows into the crystalline basement and their controls on groundwater flow patterns across the Paradox Basin, western USA Open Access

Hydrogeologists typically neglect the crystalline basement when developing conceptual and quantitative models of regional groundwater flow systems (Taucare et al., 2020; Meyers et al., 2021). The permeability of siliciclastic and carbonate rocks is typically assumed to be orders of magnitude higher than that of fractured crystalline basement rocks. However, convective heat-flow anomaly data and inferences of solute mass transport gleaned from metamorphic systems suggest that the crystalline basement is permeable to depths of 10 km (Manning and Ingebritsen, 1999; Ingebritsen and Manning, 2010). Fluid residence times, constrained by noble gases, show enhanced permeability within the upper 1 km of Precambrian basement rocks (Ferguson et al., 2023). Continental-scale compilations of pore-fluid stable isotopic data show deeper meteoric water circulation (up to ~5 km depth) in areas of relatively high topographic relief (McIntosh and Ferguson, 2021). Tertiary-age, δ¹⁸O-depleted plutonic rocks in western North America presented by Gregory et al. (1989) argue for meteoric fluid circulation to depths of ~10 km. Deep circulation within the crystalline basement has also been shown to have important implications for mountain-front recharge (Frisbee et al., 2017), deep subsurface microbial activity (Lollar et al., 2019), and near-surface ecosystem health. Deep circulation systems are a global phenomenon (e.g., Stober, 1996; Stober and Bucher, 2004, 2007, 2015a; Bucher et al., 2009; Stober et al., 1999; Bucher et al., 2009; Stober et al., 2016) and explain most thermal anomalies, mineralized springs, and changes in regional hydrochemistry of aquifer systems in contact with basement rocks.

This study focuses on understanding groundwater flow interactions between the crystalline basement and sedimentary units of the Paradox Basin, Utah, USA. The Paradox Basin hydrogeologic system is conceptualized as having an upper and lower aquifer system separated by the Paradox Formation confining unit (Fig. 1C; Thackston et al., 1981). The upper aquifer system includes the Navajo Sandstone, Burro Canyon, Cutler, and Honaker Trail formations. The Redwall Limestone is the principle aquifer of the lower aquifer systems. The Paradox Formation, comprised of evaporites and organic-rich shales (Nuccio and Condon, 1996), acts as a tight confining unit that separates the two aquifer systems. In some regions across the Colorado Plateau, the Paradox Formation acts as a seal, trapping CO₂ and He within underlying reservoir rocks (Heath et al., 2017; Tyne et al., 2022). Across the Paradox Basin, faults act as conduits for hydrocarbons and CO₂ (Shipton et al., 2004), and in the geologic past, also for ore-forming fluids (Jacobs and Kerr, 1965; Chan et al., 2000; Chan et al., 2001; Bailey et al., 2022). Faults are also associated with modern springs (Thackston et al., 1981).

What is perhaps less appreciated is the role of relatively permeable crystalline basement rocks in controlling recharge and deep regional groundwater flow patterns across the Paradox Basin. Geochemical, isotopic, and noble gas data presented by Kim et al. (2022a) and Tyne et al. (2022) indicate that meteoric fluids mixed with Paradox Formation brines within the Redwall Limestone and McCracken Sandstone. This could only be accomplished if the La Sal Mountains were acting as a recharge tower for the lower aquifer system. Here, we define the crystalline basement as including both Precambrian igneous/metamorphic and 28 Ma laccolith units (Condon, 1997). Triggered seismicity related to saline brine injection within the crystalline basement 1.2 km beneath the Redwall Limestone (Ake et al., 2005) suggests that the crystalline basement has non-negligible permeability (Zhang et al., 2013). The injection rates reported by Ake et al. (2005) are too high to be accommodated by low-permeability basement rocks. Crystalline basement-hosted springs discharging along the Colorado River in the Grand Canyon are additional evidence for a permeable crust (Crossey et al., 2009). In this study, we hypothesize that where the crystalline basement outcrops, such as within the La Sal and Abajo Mountains (Fig. 1A), it acts as a hydrologic window, permitting meteoric fluids to descend beneath the Paradox Formation confining unit and recharge the lower aquifer system. The concept of hydrologic windows was first proposed to explain the locations of groundwater up-flow zones associated with crystalline basement-hosted hot springs along the Rio Grand Rift in New Mexico (Barroll and Reiter, 1990; Mailloux et al., 1999; Pepin et al., 2012). However, we argue below that this conceptual model equally applies to groundwater recharge areas.

To assess groundwater flow interactions between the crystalline basement and overlying/adjacent sedimentary units within the Paradox Basin, we constructed a suite of cross-sectional paleo-hydrogeologic models. The models solve for variable-density groundwater flow, heat, and solute transport. We also tracked advective-dispersive isotopic (δ¹⁸O) transport and groundwater residence times (Goode, 1996). Because we compared our model results to geochemical and isotopic data presented by Kim et al. (2022a), our study is focused on the Pleistocene to Modern groundwater flow system originating within the La Sal Mountains and flowing southeast across Lisbon Valley, Utah. We allowed water-table elevations and the oxygen isotopic composition of recharge within the La Sal Mountains to vary between glacial cycles over a 1 m.y. period in an attempt to replicate Pleistocene climate forcing.

We addressed the following questions: How permeable is the crystalline basement within the Paradox Basin? To what extent does the crystalline basement modify the transport of geochemical tracers within the lower aquifer system? How do faults and fault-zone connectivity influence groundwater flow patterns? Have the groundwater flow system and geochemical/isotopic tracers within the Paradox Basin arrived at a dynamic equilibrium with the modern climate?

The sedimentary units of the Paradox Basin discussed in this study are listed in Table 1. Prior to formation of the Paradox Basin during Cambrian and Devonian times, marine shales, sandstones, and carbonates were deposited, including the McCracken Sandstone and the Ouray Limestone, in a continental platform environment (Condon, 1997). The Ordovician and Silurian eras were a period of nondeposition (Condon, 1997). During the Mississippian, a marine transgression resulted in deposition of the Redwall Limestone. Paradox Basin deposition began in earnest during Pennsylvanian–Permian times with the ancestral Rocky Mountains and Uncompahgre uplift. During the Pennsylvanian, the Paradox Formation was deposited to the west of the Uncompahgre uplift. Organic-rich shale, dolomite, and evaporite units were deposited in a marginal marine environment during cycles of marine flooding and regression. Goldhammer et al. (1991) defined 34 Milankovich-driven carbonate cycles during the middle Pennsylvanian. Ductile deformation of the salt beds occurred as sediment was shed off the Uncompahgre uplift (Barbeau, 2003) to create the Honaker Trail and Cutler formations. Salt tectonics led to the development of a series of mini-basins (Rasmussen and Rasmussen, 2009), including Lisbon Valley. This was followed by the deposition of continental units, including the Mesa Verde Group, which includes the aeolian Navajo Sandstone and the Morrison and Burro Canyon formations. During the Cretaceous, western North America was inundated by a shallow sea, which resulted in the deposition of the Mancos Shale. Up to 2 km of marine shales were deposited across the Colorado Plateau during the Cretaceous. The Cenozoic was a period of relative stability and nondeposition (Murray et al., 2016). During the Oligocene, a series of laccoliths were emplaced across the Colorado Plateau, including the La Sal complex (Hunt and Waters, 1958). Rapid erosion over the past 2–5 m.y. (Murray et al., 2016) associated with downcutting of the Colorado River and its tributaries, including the Dolores River, resulted in the formation of the La Sal Mountains, which initiated a regional, topographically driven groundwater flow system.

In La Sal, Utah (elevation 2127 m), a community on the southern side of the La Sal Mountains close to Lisbon Valley, the mean annual temperature is 8.4 °C, and the precipitation is 487 mm (Noyes et al., 2021). Modern evapotranspiration rates exceed precipitation within the lowlands (Table S1¹). Annual precipitation within the La Sal Mountains (elevation 3880 m) is up to 0.83 m (Richmond, 1972). Present-day and paleo-recharge rates were estimated in this study using monthly temperature and precipitation data described in the Supplemental Material. In the desert southwest during the Last Glacial Maximum (LGM), temperatures are estimated to have been at least 5–7 °C cooler, and precipitation is thought to have doubled in New Mexico, USA, which led to the formation of a number of Pleistocene lakes (Benson, 1988; Phillips et al., 1986; Menking et al., 2004; Allen, 2005; Asermon et al., 2010; Reheis et al., 2014). We hypothesize that isotopic/chemical tracers were modified by Pleistocene climatic cycles. Table S1 estimates how temperature, evapotranspiration, recharge, and the stable oxygen isotopic composition of water (δ¹⁸O) in precipitation may have varied between modern times and the LGM. Figure S1 plots changes in simulated isotopic composition at the land surface within the La Sal Mountains and beneath Lisbon Valley during the Pleistocene. Some prior studies indicated that recharge may have been up to three times greater than Holocene levels during the LGM (Zhu et al., 2003). It is likely that some of this available excess water increased runoff during periods of glaciation (Putnam and Broecker, 2017) more than diffuse recharge. In this study, we assume that LGM temperature reduction resulted in ¹⁸O-depleted recharge (Noyes et al., 2021).

Hanshaw and Hill (1969) presented hydrologic and geochemical analyses of the Paradox Formation, Redwall Limestone, and Cutler and Honaker Trail formations (Fig. 2A). They found that, with few exceptions, the salinity (i.e., total dissolved solids) of formation waters in the upper aquifer units ranged from fresh (<1 ppt) to brackish (<10 ppt). However, the Paradox Formation contained brines with up to 400 ppt of salinity. As noted above, Thackston et al. (1981) conceptualized Paradox Basin hydrogeology as having upper and lower aquifer systems separated by the Paradox Formation (evaporites), which serves as the regional confining unit (orange pattern, Fig. 1C). Analysis of head maps (Figs. 1A and 1B) indicates a consistently downward vertical head gradient (blue circles, Fig. 1C). Drill-stem test data (Allis, 2013) reveal fluid underpressure levels of up to 20 MPa (2000 m) within the Paradox Formation (Fig. 1D). This is likely due, in part, to rapid uplift and erosion (Corbet and Bethke, 1992) that began ~2–5 m.y. ago (Murray et al., 2016) or possibly dehydration reactions (Stober and Bucher, 2004). Underpressures within the upper and lower aquifer systems are also likely due to the effects of topography-driven flow (Belitz and Bredehoeft, 1988). Plan-view contour maps of upper (Fig. 1A) and lower (Fig. 1B) piezometric surfaces indicate that the La Sal and Abajo mountains (laccolith intrusions) are acting as recharge towers related to orographic precipitation effects. The lack of upward vertical head gradients in topographically low-lying regions suggests that much of the Paradox Basin acts as a recharge area for the greater Colorado Plateau.

Reitman et al. (2014) developed a three-dimensional model of variable-density groundwater flow and solute transport to quantify salt loading into the Colorado River within Gypsum Valley, which is located ~100 km to the southwest of Lisbon Valley and receives recharge from the Abajo Mountains. They performed a model calibration exercise using water levels and salinity data and estimated that the annual mass of salt dissolved within the upper aquifer system contributed ~2.2 × 10⁵ kg to the Colorado River in Gypsum Valley. Gardner et al. (2020) used a suite of geochemical and isotopic tracers to estimate the location and magnitude of recharge from the La Sal Mountains to Spanish Valley, which is near Moab (Fig. 1A). The principle aquifer in their study area is the Glenn Canyon Group, which includes the Navajo Sandstone. Groundwater flow was from the La Sal Mountains, with discharge into the Colorado River in low-lying areas near Moab, Utah. These authors concluded that there is little depression-focused recharge along arroyos, with the bulk of the recharge occurring at high elevations within the La Sal Mountains. Because the Glenn Canyon Group does not crop out in the uplands, they concluded that much of the recharge is provided via the crystalline basement. Their estimated recharge rate, using a lumped-parameter model, was ~0.09 m/yr. Corrected ¹⁴C ages for the Glenn Canyon Group aquifer ranged between 1700 years and 3700 years, with a mean residence time of 2700 years.

Noyes et al. (2021) used water-well levels and isotopic tracer data to assess the hydrologic connection between the Burro Canyon and Navajo aquifers within Lisbon Valley. Water-level data indicated a relatively high vertical hydraulic gradient (~6.7) between the Navajo and Burrow Canyon aquifers, which are separated by the Morrison Formation. This hydraulic gradient is three times what is reported in Figure 1C. Stable isotopic compositions of water (δ¹⁸O and δ²H) and ¹⁴C ages in these two aquifers are distinctive. The Burro Canyon pore fluids are Holocene in age (11–3.3 ka), while fluids within the Navajo Sandstone are late Pleistocene in age (36–15 ka). Scatter plots of δ¹⁸O, δ²H, and ¹⁴C data indicated that the older groundwater within the Navajo aquifer is isotopically depleted, which is consistent with recharge under cooler conditions during the late Pleistocene. The Navajo aquifer crops out at a higher elevation in the foothills of the La Sal Mountains, which could also partially explain the more depleted isotopic composition. Noyes et al. (2021) concluded that there was little hydrologic communication between the Burro Canyon and Navajo aquifers.

Permeability data for different formations within the Paradox Basin can be found in Woodward-Clyde Consultants (1982), Freethey and Cordy (1991), and Lopes and Hoffmann (1997). We summarize the data in Table 2. For the Paradox Formation evaporites, we relied on measurements from field-pressure tests carried out in situ in bedded salt at the medium-level nuclear waste repository near Carlsbad, New Mexico, as reported in Beauheim and Roberts (2002). For Mancos Shale permeability, we relied on laboratory core measurements reported in Gutierrez et al. (2015). It is worth noting that core measurements can underestimate permeability (Stober and Bucher, 2015a). Neuzil (1994) pointed out that lab measurements were 10⁻²⁰ m² compared to the 10⁻¹⁶ m²Bredehoeft et al. (1983) estimated at the regional scale. Porosity data for the Paradox Basin sedimentary units varied between ~0.05 and 0.2 (Cappa and Rice, 1995; Chidsey et al., 2003; Woodward-Clyde Consultants, 1982; Clem and Brown, 1984). We assigned a relatively high value of porosity for the crystalline basement rocks (0.05; Table 1). However, similar porosity values (0.023) have been inferred from tracer tests within the crystalline basement along the Rhine Graben at depths of 2–4 km (Aquilina et al., 2004).

We constructed NW–SE cross-sectional hydrothermal models (FEMOC; Person et al., 2007) from the La Sal Mountains across Lisbon Valley (red line A–A′ in Fig. 1A). The model domain has a maximum thickness of 7.3 km and is 36 km in lateral extent (Fig. 2A). We included ~2 km of crystalline basement beneath the basal aquifer group. The upper aquifer system includes the Honaker Trail, Cutler, Navajo Sandstone, and Burro Canyon formations. These aquifer units are separated by confining units that include the Chinle-Moenkopi formations, Morrison Formation, and Mancos Shale. The Mancos Shale is absent in portions of Lisbon Valley. The lower aquifer system includes the McCracken Sandstone and Redwall Limestone. These two units were lumped into a single unit in our model. Also included in our model are the La Sal Laccolith and underlying granitic/metamorphic basement rocks (Condon, 1997). We included two fault zones in our study. The crystalline basement varied between 10⁻¹⁶ m² and 10⁻¹⁴ m² in the two fault scenarios. The faults were assigned fault permeabilities two orders of magnitude higher in the z-direction (10⁻¹⁴ m²) than in the x-direction (10⁻¹⁶ m²). These subvertical faults extend through the Paradox Formation to the base of the model domain (solid red lines in Fig. 2A). The GTO and Lisbon Valley fault zone elements have widths of 98 m and 133 m, respectively. The model is comprised of 2846 nodes and 5466 triangular elements. Near the La Sal Mountains, elements have a maximum width of ~770 m. Vertical discretization varied between ~110 m (sedimentary units) and 760 m (crystalline basement elements). We considered additional fault scenarios (not shown) where the faults terminated within the Paradox Formation. We found that they had little effect on the salinity within the crystalline basement.

We solved a variable-density groundwater flow equation (Equation A1 in the Supplemental Material). The dependent variable is the equivalent freshwater head (Fig. 3). We included a sink term in the groundwater flow equation to approximate the development of underpressure due to erosional unloading (Corbet and Bethke, 1992). We did not remove sediments (nodes) along the top surface of the model domain during the 1 m.y. simulation period; for the fluid sink term, we assumed an erosion rate of 0.4 mm/yr, which is consistent with the removal of 2 km of Mancos Shale over a period of 5 m.y. (Table 3). We assigned a relatively high specific storage coefficient of 3 × 10⁻⁵ m⁻¹ to accentuate the development of underpressure within the Paradox Formation (Fig. 1D). This high specific storage coefficient had relatively little effect on computed transient heads within aquifer units.

Equivalent freshwater heads become high if brines are present, but they cannot be used to directly infer directions of vertical groundwater flow because of buoyancy effects (Post et al., 2007). We solved a conductive-convective heat transport equation as well as a series of advective-dispersive transport equations, including transport of solute, isotopic tracers (δ¹⁸O), and mean groundwater residence times (Equations A3–A8 in the Supplemental Material). We neglected fluid-rock isotope exchange reactions due to the relatively shallow depth and associated low temperatures (<160 °C). The equations were solved using the finite element method. We used the modified method of characteristics to approximate advective transport.

We imposed specified-value boundary conditions along the top boundary for heat, solute, and isotopic transport, along with groundwater residence time (Fig. 2A). The bottom boundary conditions were all no flux except for heat transport (Fig. 2A). We used a specified heat flux of 60 mW/m² along the bottom of the model domain. Within the La Sal Mountains, we allowed specified heads to fluctuate by up to 20 m during glacial–interglacial periods of the Pleistocene (Paces et al., 2020). This resulted in only a small increase in recharge, far less than that reported by Zhu et al. (2003). We applied a specified head along the southern edge of our model domain to allow fluids to exit Lisbon Valley (Figs. 2B and 3). Heads decreased from 1830 m to 1400 m between the top and base of the model domain except along the Paradox Formation (Fig. 2B). The imposed decrease in head with depth is consistent with the potentiometric maps of Thackston et al. (1981). By not imposing a specified head along the Paradox Formation, we allowed underpressures to develop within the Paradox Formation that were not influenced by this boundary. A spring boundary (also known as a no-diffusive-flux boundary) was imposed for the other transport equations along this edge. No-flux boundaries were imposed along the northern edge of the model domain.

We assumed local hydrostatic initial conditions for groundwater flow, and imposed a linear increase in salinity, temperature, and mean groundwater residence time with depth. Initial salinity conditions increased with depth from 0 ppt at the land surface to 300 ppt at the base of the model domain. Initial pore fluid residence time increased linearly from 0 m.y. to 4 m.y. between the top and base of the model domain. Initial δ¹⁸O varied from between about −13‰ to −17‰ at the land surface (Table S1) to about +9.3‰ at the base of the model domain. Within the Paradox Formation, we fixed salinity and δ¹⁸O values to be equal to 300 ppt and 5‰, respectively, during the simulation. Initial temperatures increased linearly using a subsurface geothermal gradient of 30 °C/km. We used computed initial conductive temperatures and equilibrium fractionation factors for a mineral assemblage that included quartz, anorthite, muscovite, biotite, hornblende, and calcite to set the initial δ¹⁸O values of the fluids (Bowman et al., 1994). In our model, the modern mean annual land surface temperatures varied with elevation between 11.5 °C and 1.0 °C, given the change of 2100 m in elevation between the Lisbon Valley and the La Sal Mountains. During glacial times, δ¹⁸O values were decreased by 6‰ due to a 6 °C temperature reduction along the top boundary (water table). All simulations were initialized (spun up) and run for 1 m.y. to ensure that the initial salinity, residence times, and initial oxygen isotopic conditions would not have a significant impact on present-day model results. The models were then run at 1.05 m.y. using a time step size of 100 yr; solute and isotopic tracers have established dynamic equilibrium conditions by the end of the simulation. Model runs required up to two weeks of simulation time on our Linux cluster depending on the permeability level of the crystalline basement that was assigned.

We previously ran a number of simulations, varying the permeabilities of the upper and lower aquifer and confining unit. Some of these results can be found in Noyes (2019). In this study, we did not vary the permeability of the aquifer and confining units. Rather, our analysis focused on the effects of crystalline basement permeability and the presence or absence of faults on groundwater flow between the La Sal Mountains and Lisbon Valley. Based on ⁸¹Kr and δ¹⁸O measurements reported by Kim et al. (2022b), relatively young (ca. 1 Ma) meteoric fluids occur within the basal aquifer group. Noble gas results show extensive flushing of remnant basinal brines by meteoric recharge (Tyne et al., 2022). This would only be possible if the La Sal Laccoliths were sufficiently permeable to permit significant volumes of meteoric water to percolate down and mix with basin brines beneath the Paradox Formation.

We used ⁸¹Kr and ¹⁴C age tracers, salinity, and δ¹⁸O data reported by Noyes et al. (2021) and Kim et al. (2022a, 2022b) to constrain and test our model results. Permeability and porosity values assigned to each of the 15 hydrostratigraphic units are presented in Table 1. The thermal conductivity of the Paradox Formation was set about twice as high as that of the clastic and carbonate units. Modeled scenarios and parameters considered in our sensitivity study are listed in Table 4.

Figure 3 presents contour maps of freshwater heads computed for the five scenarios considered in our sensitivity study (Table 4). Due to erosional unloading (0.4 mm/yr), the Paradox Formation had heads below hydrostatic conditions (up to ~970 m below hydrostatic conditions, or ~10 MPa in all modeled scenarios). Topography-driven flow dominates within the upper aquifer system. The northwest to southeast trend of increased hydraulic heads within the crystalline basement beneath the Paradox Formation is due to the increasing permeability of the La Sal Laccolith. As the permeability of the crystalline basement increased from 10⁻¹⁸ m² to 10⁻¹⁴ m² (Figs. 3A–3C), elevated heads propagated southward along the bottom 2 km of the model domain beneath the lower aquifer system. Because of the specified head boundary condition along the right (southern) edge of the model domain (Fig. 1B), groundwater migrated out of Lisbon Valley. Had we chosen a no-flow boundary for the entire right edge of the model domain, Lisbon Valley would have become a groundwater discharge area with upward hydraulic head gradients, which is inconsistent with the water-level measurements of Thackston et al. (1981) and Noyes et al. (2021). The presence of faults that cut the Paradox Formation in scenarios 4 and 5 allowed groundwater from the upper aquifer system to migrate down into the Paradox Formation, creating underpressured cells (Figs. 3D and 3E). Figure 4A presents vertical changes in hydraulic head beneath Lisbon Valley at x = 31 km; (vertical gray line in Fig. 2A indicates location of profile). Within the upper and lower aquifer systems, heads computed for all scenarios compare reasonably well to the estimated range of water levels within Lisbon Valley from the piezometric contour maps (Figs. 1A and 1B) of Thackston et al. (1981; see horizontal black lines in Fig. 4A). Computed heads within the upper aquifer system only agree with the shallowest portion of the hydraulic head data reported in Noyes et al. (2021). Note that only scenario 5 (green dashed line; crystalline basement permeability of 10⁻¹⁴ m²) produced a downward hydraulic gradient near the land surface that is consistent with field observations.

Flow rates and directions were sensitive to crystalline basement permeability (Fig. 5). A relatively small volume of deep recharge from the La Sal Mountains was focused into the lower aquifer system when the crystalline basement permeability was relatively low (10⁻¹⁸ m²; Figs. 4B and 5A); groundwater velocities (Darcy flux [] divided by porosity [ϕ]) within the lower aquifer system were only ~0.0028 m/yr for the low-permeability scenario (Fig. 4B), and 0.0005 m/yr within the underlying tight crystalline basement. Recharge into the upper aquifer system was controlled by the elevation where clastic units such as the Navajo Sandstone outcropped. When the permeability of the crystalline basement was raised to 10⁻¹⁶ m² or 10⁻¹⁴ m² (Table 4, scenarios 2 and 3; Figs. 5B and 5C), the flow rates in the lower aquifer system increased to ~0.1 m/yr (Fig. 4B), as this unit received significant recharge from the La Sal Mountains. When the crystalline basement was assigned a permeability of 10⁻¹⁶ m² or 10⁻¹⁴ m², recharge to the units of the upper aquifer system came not only from the La Sal Mountains across the water table but also from lateral flow below the land surface (Figs. 5B and 5C). The vertical velocity at the water table (top surface) within the La Sal Mountains for the high-permeability scenario was 4 m/yr. Multiplying this by porosity (0.05) yields an estimated diffuse recharge rate (q_z) of 0.19 m/yr. Gardner et al. (2020) estimated a recharge rate of 0.09 m/yr based on lumped-parameter modeling and ¹⁴C groundwater ages. For the intermediate permeability scenario (10⁻¹⁶ m²), the vertical velocity of groundwater within the La Sal Mountains decreased to ~0.12 m/yr, lowering recharge to 0.006 m/yr, which is low in comparison to the rate of Gardner et al. (2020). Groundwater velocities within the crystalline basement increased from ~0.0002 m/yr to 0.05 m/yr to 4 m/yr as the crystalline basement permeability increased from 10⁻¹⁸ m² to 10⁻¹⁶ m² and 10⁻¹⁴ m², respectively (Figs. 4B, 5B, and 5C). Groundwater streak lines (red arrows) in Figure 5 indicate that nearly all flow exits the model domain along the southern boundary. Groundwater flow directions are controlled, in large part, by permeability and the lateral head gradient between the La Sal Mountains (~2700 m) and the southern edge of the model domain (1830–1400 m). Groundwater velocities exceeded 1 m/yr within the permeable units of the upper aquifer system. The inclusion of anisotropic faults (k_z > k_x) did little to change the fluid flux much within the upper aquifer system. For fault scenario 4 (Table 4), fluids moved into the GTO and Lisbon Valley faults within the upper aquifer system and migrated downward, terminating within the Paradox Formation (Fig. 5D) owing to the inward-directed hydraulic gradients within this underpressured formation (Corbet and Bethke, 1992). In fault scenario 5, some of the fluids entering the fault zone above the Paradox Formation migrated into the lower aquifer system (Fig. 5E).

Figure 6 presents computed salinity patterns for all modeled scenarios. The position of the mixing zone within the upper aquifer system was found to be sensitive to both crystalline basement permeability and the presence or absence of faults. Within the shallow units of the upper aquifer system, a topography-driven flow system extends down into the Cutler and Honaker Trail formations, maintaining low salinities. Salinities within the Cutler Formation increased along the flow path to the southeast as laccolith and crystalline basement permeability increased (Figs. 6A–6C). For high-permeability scenario 5, in which faults were added and horizontal permeability was lower (10⁻¹⁶ m²) than that of the aquifers, the freshwater–saline water mixing zone rose into the Cutler Formation (Fig. 6E). The salinity in the lower aquifer system ranged between 100 ppt and 300 ppt, depending on the permeability of the crystalline basement. For the lowest crystalline basement permeability of 10⁻¹⁸ m², the crystalline basement and lower aquifer system is dominated by high salinity (~300 ppt; Figs. 6A and 7A). Simulations that considered higher permeabilities diluted salinity in the lower aquifer system and crystalline basement (Figs. 6B, 6C, and 7A). The intermediate permeability scenario is most consistent with measured salinities within the lower aquifer system across Lisbon Valley (Fig. 7A). The highest permeability scenario has the lowest computed salinity in other regions (e.g., x = 10 km; Fig. 6C). Transient haline-convection cells formed along the flow path within the crystalline basement (Fig. S1). Computed salinities are in reasonably good agreement with conditions observed for the intermediate- and high-permeability scenarios (blue dots, Fig. 7A; Kim et al., 2022a).

Computed mean groundwater ages are presented in Figure 8. Observed groundwater residence times within the upper aquifer system ranged from Pleistocene to Holocene (Fig. 7B). As crystalline basement permeability increased (Figs. 8A–8C), the volumes of relatively young meteoric fluids entering the Cutler and Honaker Trail formations rose. For scenario 4, faults with relatively low horizontal permeability had a barrier effect, increasing simulated groundwater age (Fig. 8D). Groundwater age within the Paradox Formation ranged between 4 Ma and 2 Ma. For the low-permeability crystalline basement scenarios (10⁻¹⁸ m²; Fig. 7B), groundwater age continued to increase in a nearly monotonic trend from the Paradox Formation to the bottom of the model domain. As crystalline basement permeability increased, relatively younger groundwater was introduced beneath the Paradox Formation. Simulated ages for the intermediate- and high-permeability basement scenarios are consistent with ⁸¹Kr groundwater ages measured within Lisbon Valley (Fig. 7B; Kim et al., 2022b). Simulated groundwater ages within the lower aquifer system are in reasonably good agreement with observed conditions for the intermediate- and high-permeability scenarios (Fig. 7B).

Simulated δ¹⁸O values within the shallow aquifer system varied between about −13‰ to −18‰ (Fig. 9), which is consistent with the groundwater isotopic compositions measured by Noyes et al. (2021). The δ¹⁸O values in the Paradox Formation were fixed at 5‰. Within the lower aquifer system and underlying crystalline basement, mixing between ¹⁸O-enriched fluids from the Paradox Formation and relatively δ¹⁸O-depleted meteoric recharge resulted in a net range of isotopic composition of between ~−8‰ to +5‰ (Fig. 7C). Within the lower aquifer, simulated values of δ¹⁸O for the high-permeability and fault scenarios came closest to matching the average δ¹⁸O value reported by Kim et al. (2022a). For the low-permeability crystalline basement scenario (10⁻¹⁸ m²), the δ¹⁸O values increased with depth to ~9‰ at the base of the model domain (Fig. 7C). Figure S1 compares temporal trends in δ¹⁸O at the upper surface of the model near the top of the La Sal Mountains and Lisbon Valley within the Burro Canyon Formation, the lower aquifer system, and the crystalline basement (see white dots in Fig. 2A). For the high-permeability scenario, the effects of transient thermohaline convection cells on simulated δ¹⁸O can be seen in the temporal trends in the isotopic composition of fluids within the crystalline basement; the transient thermohaline convection cells have a much shorter period than the climate forcing cells (Fig. S1C). In the Burro Canyon Formation (Fig. S1A), where flow variations are controlled by fluctuations in the water table, there are longer period and lower amplitude δ¹⁸O variations than in the deeper aquifer (Fig. S1B). Long period, thermohaline convection developed in the intermediate permeability scenario within the lower aquifer system (Fig. S1B). Temporal variations in δ¹⁸O within the crystalline basement (Fig. S1C) are observed for the high permeability scenario (10⁻¹⁴ m²).

Figure 10 presents computed temperatures for scenarios 1–5. Simulated temperatures are influenced by both convective heat transfer effects and the thermal conductivity contrast between the Paradox Formation and other units (Table 1). The absence of the Paradox Formation by the La Sal Laccolith created a complexity of simulated temperature patterns (Fig. 10A). The Paradox Formation is cut by the La Sal Laccolith between x = 0–8 km created complexity in simulated conductive temperature patterns (Fig. 10A). The bulk thermal conductivity Paradox Formation (5.0 W-m/°C) is about twice that of the La Sal Laccolith unit. Increases in laccolith and crystalline basement permeability resulted in convective cooling beneath the La Sal Mountains (note the change in position of the 50 °C isotherm between Figs. 10A–10C). Within Lisbon Valley, conductive heat transfer erased convective effects along the fault zones (Figs. 10D and 10E; Person et al., 2007). Figure 7D compares simulated temperature profiles at 0 km, 20 km, and 31 km along the model cross section to oil well temperature measurements collected during oil well shut-in tests conducted by Allis (2013). The broad range of temperatures below 2 km depth (43–104 °C) and changes in temperature gradients are largely due to the effects of thermal conductivity contrasts between the Paradox Formation and other units rather than convective effects (Fig. 7D). We are unaware of hot springs reported within the Paradox Basin. The change in slope of temperatures with depth occurs within the thermally conductive Paradox Formation.

This study demonstrates that laccoliths cutting sedimentary confining units created important pathways (hydrologic windows) for groundwater recharge into the lower aquifer system of the Paradox Basin and underlying crystalline basement rocks. Our hypothesis of deep groundwater circulation beneath the Colorado Plateau is not new. Crossey et al. (2009) used ³He/⁴He and ⁸⁷Sr/⁸⁶Sr data to argue that crystalline basement-hosted springs within the great unconformity in the Grand Canyon region are associated with deep-flow system scavenging, mantle-derived ³He, and radiogenic Sr. Our findings that the upper aquifer system in the Paradox Basin is being recharged via the La Sal Mountain block are supported, in part, by the findings of Gardner et al. (2020). Interpreting the δ¹⁸O data, these authors concluded that recharge to the Glenn Canyon Group, which includes the Navajo Sandstone, is being supplied via the fractured La Sal Laccolith rocks on the western side of the La Sal Mountains (see their fig. 10; Gardner et al., 2020). The observed salinity, groundwater residence times, and isotopic composition of pore fluids measured within the lower aquifer system by Kim et al. (2022a, 2022b) are consistent with modeled scenario results that assigned crystalline basement permeability in the intermediate–high range (10⁻¹⁶ m² to 10⁻¹⁴ m²).

To test our hypothesis that the crystalline basement underlying the Paradox Basin is relatively permeable, we developed a simple one-dimensional analytical model of triggered seismicity. Ake et al. (2005) reported triggered seismicity after ~110 days of continuous brine injection into the Redwall Limestone at a rate of 1290 L/min within Paradox Valley. They indicated that the average formation pressure at the wellhead rose from ~42 MPa (hydrostatic) to 80 MPa (ΔP = 38 MPa). The seismicity that was triggered occurred to ~5 km lateral distance and to a depth of up to ~1.2 km beneath the Redwall Limestone (Fig. 11A). To approximate the downward propagation of a pressure front beneath the Redwall Limestone, we used the following analytical model:

where h(d,t) is the computed, time-dependent anomalous head; d is depth below the injection horizon; erfc is the complementary error function; h_o is the value of elevated head in response to fluid injection (at d = 0, t > 0); t is time; S_s is specific storage (m⁻¹); and D_h is the hydraulic diffusivity (K/S_s; m²s⁻¹), where K is hydraulic conductivity, which is a function of permeability (K = kρ_fg/μ_f; k is permeability, ρ_f is fluid density, μ_f is fluid viscosity, and g is gravity; ms⁻²). At time zero, heads/anomalous pressures are hydrostatic, i.e., 0 MPa, h(d,t = 0). For t > 0, the head at d = 0 was instantaneously increased to 38 MPa (3800 m). We report our results in equivalent anomalous pressures (MPa) rather than heads. Figure 11B presents earthquake foci as well as computed anomalous pressures after 110 days using basement permeabilities of 10⁻¹⁸ m², 10⁻¹⁶ m², and 10⁻¹⁴ m². We assumed that S_s was 10⁻⁶ m⁻¹ within the crystalline basement. Pressure anomalies greater than 1 MPa occurring at a depth of 1 km are more than sufficient to trigger seismicity (Ge et al., 2009). The analytical solution results are most consistent with a crystalline basement permeability of 10⁻¹⁴ m².

Flow within the crystalline basement to depths of 7 km is certainly possible according to Ingebritsen and Manning (2010) and Manning and Ingebritsen (1999), who suggest that permeability at 10 km depth can be as high as 10⁻¹⁶ m² in geothermal and metamorphic environments. Precambrian basement rocks have relatively high permeabilities in the upper 1 km (10⁻¹⁷ to 10⁻¹⁴ m²) and lower permeabilities at greater depths (<10⁻¹⁸ m²), based on noble gas residence time tracers (Ferguson et al., 2023). In mountainous regions, topographically driven flow can drive meteoric fluids to depths of up to 5 km, based on stable water isotopes (McIntosh and Ferguson, 2021). The studies above primarily used geophysical and geochemical/isotopic datasets to arrive at their conclusions. Hydraulic tests conducted within deep boreholes also indicate relatively high crystalline basement permeability (Stober and Bucher, 2015a).

We computed the solute mass flux that exits Lisbon Valley (x = 31 km) for all sensitivity study simulations (Fig. 4C). The solute mass flux within various aquifers in the upper aquifer system varied between 1.2 × 10⁴ kg/yr to 4.6 × 10⁴ kg/yr. Solute mass flux migrating out of Lisbon Valley beneath the Paradox Formation was ~1.5 × 10² and 2 × 10⁴ kg/yr for the low- and intermediate-permeability scenarios (10⁻¹⁸ m² to 10⁻¹⁶ m²), respectively. For the high basement permeability scenario (10⁻¹⁴ m²), the solute mass flux exiting Lisbon Valley was 8.5 × 10⁴ kg/yr, similar to that exiting the upper aquifer system. Reitman et al. (2014), using a three-dimensional mathematical model (SUTRA), estimated that the mass of salt discharging into the Colorado River from the Honaker Trail and Cutler formations within Gypsum Valley was 2.3 × 10⁵ kg/yr. Where is this salt going? Within the upper aquifer system, the salt load is likely migrating toward the Dolores River to the north, which has an annual Cl flux of 1.3 × 10⁸ kg/yr (Hite and Lohman, 1973). Within the crystalline basement, a component of the salt load may be migrating toward lowlands along the Colorado River. Crossey et al. (2009) noted that along the Colorado and Little Colorado rivers in the Grand Canyon, there is discharge of Na-Cl–rich fluids with total dissolved solids of up to 50 ppt. We hypothesize that some of these saline fluids may be derived from the Paradox Basin. The lower aquifer system potentiometric surface near Mexican Hat is lower than the elevation of the San Juan River (Fig. 1B), which suggests that groundwater flow is migrating southwestward toward the Grand Canyon.

In some settings, faults act as seals (Stober and Bucher, 2015b). This study did not find that faults were the locus of significant groundwater transfer between the upper and lower aquifer systems. This is because the underpressured Paradox Formation was able to capture fluids migrating down faults beneath Lisbon Valley. Copper mineralization associated with fault zones within Lisbon Valley clearly indicates that faults focused vertical fluid flow in the geologic past (Jacobs and Kerr, 1965; Huntoon, 1986; Chan et al., 2000; Bailey et al., 2022). We argue that the topography-driven groundwater flow system is a relatively recent phenomena, perhaps only established in the past 2–5 m.y. (Murray et al., 2016; Kim et al., 2022b). During the Eocene, when 2 km of the Mancos Shale capped the Paradox Basin, thermohaline convection was likely the dominant mechanism driving fluid flow vertically along fault zones and may have been an important mechanism for ore mineralization.

Inspection of the time series of computed δ¹⁸O values in Figure S1 suggests that the hydrologic system approached dynamic equilibrium conditions within the crystalline basement after 1 m.y. Within the Burro Canyon Formation (Fig. S1A), simulated temporal trends in δ¹⁸O are out of phase and have a lower amplitude relative to the Pleistocene recharge signal (blue line in Fig. S1C; Loosli et al., 1998). Simulated temporal variations within the lower aquifer system and basement are controlled by the interplay between forced and thermohaline convection. It is worth noting that simulated thermohaline convection cells within the crystalline basement are approximations of actual conditions. Haline convection cells are sensitive to grid discretization (Post and Kooi, 2003) as well as heterogeneity of spatial permeability (Gerdes et al., 1995) not represented in our simulations. The long simulation times required to reconstruct the paleo-hydrogeology of Lisbon Valley over the Pleistocene prevented us from considering additional grid refinement.

We used subsurface heat and mass transport models constrained by geochemical/isotopic data from Kim et al. (2022a, 2022b) to understand the hydrologic interactions between the La Sal Mountains and Paradox Basin near Lisbon Valley. An important component of La Sal Mountain recharge enters the upper aquifer system laterally through the mountain block. A fresh–saline water mixing zone develops within the Honaker Trail and Cutler formations within the upper aquifer system. The anisotropic faults (k_z > k_x) in our model acted mainly as a barrier to lateral flow. Fluid-impelling mechanisms within the lower aquifer system and underlying crystalline basement include both topography- and density-driven flow (haline convection). Underpressures form within the low-permeability Paradox Formation due to erosion and sediment decompaction.

Importantly, we found that the La Sal Mountains act as a hydrologic window into the lower aquifer system and underlying crystalline basement. For scenarios where the crystalline basement was relatively permeable (10⁻¹⁶ m² to 10⁻¹⁴ m²), meteoric fluids mixed with brines of the Paradox Formation. Models that included a permeable crystalline basement were largely in agreement with isotopic tracers and salinity data reported by Kim et al. (2022a, 2022b). The presence of faults did not significantly modify fluid exchange between the upper and lower aquifer systems. This was due to underpressuring within the Paradox Formation (Fig. 1D). We hypothesize that the downward hydraulic gradient observed beneath Lisbon Valley is the result of a long-distance hydrologic connection to crystalline basement rocks that outcrop along the Colorado River at lower elevations perhaps as far away as the Grand Canyon. Meteoric recharge through hydrologic windows may have reintroduced microbial communities into previously sterilized sediments at the bottom of the Paradox Basin (McIntosh et al., 2023).

This study highlights the importance of groundwater circulation through the relatively permeable crystalline basement and its interactions with overlying/adjacent sedimentary basin units. Sedimentary basins should no longer be thought of as closed hydrologic systems. They have porous lower boundaries through which solutes, heat, and microorganisms (Crossey et al., 2016; McIntosh et al., 2023) are transported. Our findings may also have implications for the involvement of crystalline basement in the formation of Mississippi Valley-type ore deposits (Wilkinson, 2010).

We thank Steve Ingebritsen for his assistance in revising this manuscript. Support from the W.M. Keck Foundation (grant no. 989941) and the National Science Foundation (National Science Foundation Frontier Research in Earth Sciences grant no. 1925974 and Division of Earth Sciences grant no. 1830172) is gratefully acknowledged. Support to Jennifer McIntosh under National Science Foundation Frontier Research in Earth Sciences Subsurface Microbe-Rock-Fluid Systems grant no. 2120733 is gratefully acknowledged.