Groundwater Circulation Within the Mountain Block: Combining Flow and Transport Models With Magnetotelluric Observations to Untangle Its Nested Nature

Mountains are vital water sources for humans and ecosystems, continuously replenishing lowland aquifers through surface runoff and mountain recharge. Quantifying these fluxes and their relative importance is essential for sustainable water resource management. However, our mechanistic understanding of the flow and transport processes determining the connection between the mountain block and the basin aquifer remains limited. Traditional conceptualizations assume groundwater circulation within the mountain block is predominantly shallow. This view neglects the role of deep groundwater flowpaths significantly contributing to the water, solute, and energy budgets. Overcoming these limitations requires a holistic characterization of the multiscale nature of groundwater flow along the mountain‐to‐valley continuum. As a proof‐of‐concept, we use a coupled groundwater flow and transport model to design a series of numerical experiments that explore the role of geology, topography, and weathering rates in groundwater circulation and their resulting resistivity patterns. Our results show that accumulating solutes near stagnation zones create contrasting electrical resistivity patterns that separate local, intermediate, and regional flow cells, presenting a target for magnetotelluric observations. To demonstrate the sensitivity of magnetotelluric data to features in our resistivity models, we use the MARE2DEM electromagnetic modeling code to perform forward and inverse simulations. This study highlights the potential of magnetotelluric surveys to image the resistivity structure resulting from multiscale groundwater circulation through relatively impervious crystalline basement rocks in mountainous terrains. This capability could change our understanding of the critical zone, offering a holistic perspective that includes deep groundwater circulation and its role in conveying solutes and energy.


Introduction
Mountains are an abundant source of surface and subsurface water vital to humans and ecosystems (Viviroli et al., 2007).Natural and anthropogenic stresses, a product of unprecedented climate change, population growth, and land use, critically impact these water sources, increasing the risk of water scarcity in lowland areas, especially in arid and semi-arid environments where water supply is insufficient to meet the increasing demand (Messerli et al., 2004;Viviroli et al., 2011Viviroli et al., , 2020)).The sustainable use of "unconventional" water resources, like brackish aquifers, has emerged as an alternative to offset the ever-increasing water deficit in these arid environments (Robinson et al., 2019;Stanton & Dennehy, 2017;Stanton et al., 2017).However, these sources are intimately linked and cannot be assessed independently, as the quantity and quality of brackish resources depend on the mountain source that replenishes them.In other words, sustainable management of water resources in arid regions requires a detailed understanding of the flow and transport along the mountain-to-valley transition, particularly the unobserved and typically underappreciated flow field within the mountain block.
The mountain block is an area with significant topographic relief primarily formed by a bedrock matrix with fractured and unfractured regions (Wilson & Guan, 2004).How these mountain systems recharge lowland aquifers has been an active area of research for groundwater hydrologists over the past decades, with significant advances made in the mechanistic understanding of this process (Markovich et al., 2019).One characteristic approach divides the water that enters the lowland aquifer originating from the mountain block, known as mountain-front recharge (MFR), into two main components: (a) water-table recharge at the mountain-front zone (direct MFR) and (b) mountain-block recharge (MBR) (Wilson & Guan, 2004).MBR represents only the portion of MFR that comes from the regional groundwater sources moving within the mountain block (Markovich et al., 2019;Wilson & Guan, 2004).Both of these processes provide significant amounts of water to lowland areas and can be the primary source of recharge in arid regions (Scanlon et al., 2006(Scanlon et al., , 2012)).Over the past two decades, the hydrologic community has made fundamental advances in MBR science (Markovich et al., 2019).However, there is still considerable uncertainty when it comes to identifying and measuring the depth of groundwater circulation through mountain blocks (Condon et al., 2020;Markovich et al., 2019).This uncertainty is exacerbated by the scarcity of data for hydrogeologic parameters, such as hydraulic conductivity, porosity, and flow rates for deeper groundwater flow systems in various geologic settings (Carroll et al., 2020).
Traditional conceptualizations of groundwater circulation through the mountain block have evolved over the years (Markovich et al., 2019;Somers & McKenzie, 2020).Early ones assumed the groundwater flow only occurred in the permeable soil above the bedrock (Fan et al., 2019;Mosley, 1979;Tani, 1997); however, this assumption has been challenged in the light of new data.For example, irrigation experiments performed in the Panola Mountain Watershed in Georgia (Tromp-van Meerveld et al., 2007) and the western Cascade Range in Oregon (Graham et al., 2010) found considerable infiltration of water into the bedrock (91% in the Panola Mountain and 27% in the Cascade Range).This infiltration can circulate at different depths within the mountain block.A recent study for the Copper Creek watershed in Colorado used a hydrologic model and baseflow age estimates from dissolved gas tracers to constrain these depths, suggesting that circulation can exceed 100 m (Carroll et al., 2020).Another study estimated circulation depths from 0.2 km to over 1.2 km using quartz-silica geothermometry for two watersheds in New Mexico and Colorado, thus challenging existing hydrogeological models in these areas (Frisbee et al., 2017).Similar depth ranges of groundwater circulation were presented in another study that used the stable isotopes of water (δ 2 H and δ 18 O) in wells across North America to determine the extent of meteoric water circulation, indicating that water can flow deeper than 2 km in areas with high topographic relief (McIntosh & Ferguson, 2021).Similarly, analog studies using residence time data from environmental tracers have shown the important contribution of old groundwater to streamflow (Meyers et al., 2021;White et al., 2019White et al., , 2021)).Person et al. (2024) argued that in mountainous terrains of the Paradox Plateau, USA enhanced mountain from recharge occurs where the crystalline basement cropped out.These findings are consistent with a permeable (10 14 -10 16 m 2 ) crystalline basement and supported groundwater age dates ( 14 C, 81 Kr) and salinity inversions within lowland basal sedimentary aquifer systems (Kim et al., 2022).Hydrothermal investigations within the desert southwest of the United States suggest MBR fluid circulation depths of 3-6 km (Barroll & Reiter, 1990;Mailloux et al., 1999;Pepin et al., 2014).Lastly, simulations performed in the Swiss Alps suggest meteoric water recharge through faults of depths of ∼10 km to explain thermal anomalies (Alt-Epping et al., 2021).
From a modeling perspective, significant contributions have been made in determining the features controlling deep groundwater circulation.We know that water flowpaths through the mountain block are multi-scale and nested in nature.This behavior results from the interaction between local, intermediate, and regional flowpaths driven by the topographic features of the mountain (Tóth, 1963).While local and intermediate flowpaths recharge streams in sub-catchments within the mountain block, regional flowpaths move large amounts of water, energy, and solutes through deep systems contributing to MBR (Gomez & Wilson, 2013;Markovich et al., 2019).The latter has substantial implications on streamflow generation and the water budgets at catchment scales (Fan, 2019;Flerchinger & Cooley, 2000;Frisbee et al., 2017;Fujimoto et al., 2016;Voeckler et al., 2014).Frisbee et al. (2011) showed that the streamflow is not an accumulation of near-surface hillslope responses but a combination of the latter and flow through bedrock.Studies have used 2D and 3D numerical simulations of groundwater flow through mountainous terrains to investigate the most prominent features affecting this nested nature.These studies have determined that the vertical extent of local and regional flowpaths and their relative recharge are highly dependent on the topographic relief, the depth of the water table, and the hydraulic parameters of the mountain block (Gleeson & Manning, 2008;Welch & Allen, 2012, 2014;Welch et al., 2012).The interaction of these flowpaths also produces a mixing between young and old water within the mountain block that is translated to the lowland aquifers, producing longer tails in the residence time distributions (Frisbee, Phillips, et al., 2013;Frisbee, Wilson, et al., 2013;Gomez & Wilson, 2013).This behavior has also been shown to depend on topographic relief and variations in hydraulic parameters (Cardenas & Jiang, 2010;Jiang et al., 2010).Lastly, the interplay between flowpaths produces stagnant zones in the meeting branches of local and regional flowpaths, where solutes and heat accumulate within the mountain block (An et al., 2014;Jiang, 2012;Jiang et al., 2011Jiang et al., , 2014)).As expected, topographic relief and hydraulic conductivity also play an essential role in the location of the stagnation zones within the mountain block (Jiang et al., 2011).Samples from hydrologic systems (i.e., streams, lakes, and wells) encapsulate the effects of a myriad of flowpaths, representing a wide range of spatial and temporal time scales (Frisbee, Wilson, et al., 2013).To interpret these samples, we typically divide the catchment into active and inactive zones (e.g., Mayo et al., 2003), conceptualizing the effective circulation occurring within the shallow active region and therefore ignoring the role of deeper flowpaths (Frisbee et al., 2017), which might not be a reasonable assumption.The samples from these systems represent the interaction between the local, intermediate, and regional flow systems.Conceptualizing a shallow circulation through the mountain block will likely result in an underestimation of the transport processes occurring within the system, leading to deceptive simulation results that do not conform with the samples taken.
The use of geophysical measurements is promising, as they constitute an extensive, cost-efficient, and noninvasive vertical exploration tool of shallow and deep subsurface characteristics.Novel inversion techniques (e.g., Binley et al., 2010Binley et al., , 2015;;Ferré et al., 2009;Rubin & Hubbard, 2005) have shown that geophysical surveys can go beyond the mapping of petrophysical properties (e.g., electrical conductivity, seismic velocities, and dielectric constant) to the inference of spatiotemporal variations of hydraulic parameters (i.e., porosity and hydraulic conductivity) and state variables (i.e., solute concentration and water content, among others), providing fundamental information about hydrogeological parameters and processes (Rubin & Hubbard, 2005).The ability to detect groundwater through geophysical methods depends on the degree of saturation and how this saturation causes contrasts of bulk petrophysical properties for different lithological settings.For instance, in fully saturated systems, the p-wave velocities tend to be larger, resistivity values tend to be smaller, and permittivity values positively correlate with the water content (Kirsch, 2009).The range of probable values for each property relies on the lithological setting, the chemistry of water, and the method used (Kirsch, 2009).The use of electrical resistivity-based geophysical methods has allowed the identification of water infiltration through the vadose zone (Daily & Ramirez, 1995) and fractured bedrock (Cassiani et al., 2009), the quantification of spatial patterns of groundwater recharge (Behroozmand et al., 2019;Cook et al., 1992), the estimation of subsurface hydraulic parameters and state variables (Herckenrath et al., 2013;Hinnell et al., 2010;Pollock & Cirpka, 2012), mapping submarine and subglacial groundwater saline aquifers (Gustafson et al., 2019(Gustafson et al., , 2022)), mapping the saline lithology in the Dead Sea Basin (Meqbel et al., 2013), studying fluid transport in the lithology of the crust (Egbert et al., 2022), and the exploration of geothermal systems (Blake et al., 2016;Marwan et al., 2021;Peacock et al., 2016Peacock et al., , 2020)), among other applications.The vertical investigation extent of a survey strongly depends on the geophysical method used.For example, while electrical resistivity tomography (ERT) data allow for a shallow exploration (generally less than 500 m), magnetotelluric (MT) data can be sensitive to depths of order 100 m to more than 10 km (Figure 1a), depending on the frequency bandwidth of the data (Chave & Jones, 2012).Furthermore, the current borehole information across the United States, extracted from the Geochemical Database for the Brackish Groundwater Assessment (Qi & Harris, 2017;Stanton et al., 2017), shows abundant information above 150 m deep, leaving deeper systems unmeasured all across the U.S. (Figure 1b).The latter shows the potential of using MT to provide additional information in the deeper parts of the hydrological systems.
Information gathered from MT surveys in mountainous terrains can provide crucial information about the subsurface fluid flow and geology of the system.With fluid flow in mind, a crucial question is whether electromagnetic data can detect regional to local groundwater flow systems in mountain blocks.A study on the Dosit River in China (Jiang et al., 2014) seems to point in that direction.MT surveys performed in a cross-section of the Dosit River show contrasting resistivities that decrease with depth in areas where brackish groundwater is present (see Figures 2 and 3 in Jiang et al. (2014) and Figure S2 in Supporting Information S1).The authors argue that the differences in resistivities are not related only to the geology but to the saline groundwater circulating at those depths.This study provides important insight into how the nested nature of the flow system in mountainous terrains produces distinct resistivity patterns on MT surveys.
This work explores how MT surveys can provide information about the nested nature of flow in mountain-to-valley systems recharging lowland aquifers through regional circulation paths.To achieve this objective, we propose a handful of idealized mountain-to-valley conceptual models similar to the geologic setting of the Tularosa Basin in south-central New Mexico (see Figures S2 and S3 in Supporting Information S1).Water chemistry measurements performed in this area points to recharge from the Sacramento Mountains (Mamer et al., 2014), where water in the basin becomes older and more saline further away from the mountain block (Figure S2 in Supporting Information S1).We explore various numerical experiments that account for different types of hydrological systems, exploring the role of topography, geology, and weathering rates.The groundwater flow and transport models we used are fully coupled by equations of state that connect variations in density and viscosity with changes in temperature and pressure.We also investigated the environmental age distribution at discharge sections of the watershed and compared the predicted distribution to previous observations made in the Tularosa Basin by Eastoe and Rodney (2014) and Mamer et al. (2014) (Figure S2 in Supporting Information S1).Lastly, we used a modified version of Archie's Law proposed by Glover et al. (2000) to explore how flow and transport affect the electrical resistivity of the system.2D  (Bernard, 2003), frequency domain electromagnetics (FDEM) with 3.5 m coil instrument, time domain electromagnetics (TDEM) with 50 and 500 m loop soundings at multiple stations along a transect (Bear et al., 2013), and magnetotellurics (MT) with 1 km spacing between stations.MT can survey deeper than 2 km, and its investigation depth depends on the conductivity, duration and sampling rate of the deployment, and the upper and lower frequency limits (Chave & Jones, 2012;Spies, 1989).(b) Histogram for groundwater borehole logging depth in the U.S. (Qi & Harris, 2017;Stanton et al., 2017).forward and inverse models of resistivity computed from flow and transport were investigated using MARE2DEM, a freely available adaptive finite element code for the 2D EM modeling (Key, 2016).

Conceptual Model
Our groundwater system is a simplified conceptualization of a mountain-to-valley transition within an extensional tectonic setting.We used the Sacramento Mountains and the Tularosa Basin in south-central New Mexico as inspiration (Figures S2 and S3 in Supporting Information S1).This modeling approach is intended as a sensitivity study within an idealized rift basin, and it is in the spirit of previous efforts that draw fundamental insight from simplified domain conceptualizations (e.g., Forster & Smith, 1988;Freeze & Witherspoon, 1967;Garven & Freeze, 1984;Person & Garven, 1994).Bedrock in the Sacramento Mountains is comprised of Paleozoic carbonates, sandstone, shale, and gypsum about 1.2 km thick overlying Proterozoic metasedimentary and crystalline basement (Kelley et al., 2014;Newton & Land, 2016).The mountain block represented within our model is idealized and does not consider evaporite beds such as gypsum.No layers of Paleozoic strata are considered, and we assume that basement permeability is controlled by fracture networks.Our alluvium is divided into an alluvial fan (domain III in Figure 2a), a distal alluvial fan (domain I in Figure 2a), and an old alluvium that is cemented and hence we did not include a high permeability zone adjacent to the fault (domain II in Figure 2a).The horizontal extent of the conceptual model (Figure 2a) is 65 km divided between alluvium (L v = 52 km) and a mountain block (L m = 13 km).A fault zone characterizes the transition from mountain to valley at the piedmont of the mountain range.The fault zone has a dip angle of 55°and a thickness of 100 m.The alluvium surface has a constant slope of S v = 0.0001 from the proximal piedmont deposits near the mountain front to basin floor deposits toward the left boundary (∂Ω l ).The mountain block is simplified to follow a Tóthian-like system (Tóth, 1963(Tóth, , 2009) ) described with a sinusoidal function.The complete surface topography (i.e., boundary ∂Ω s ) in the model is represented with the following stepwise expression, where λ is the valley wavelength [L], A is the valley amplitude [L], and S v and S m are the alluvium and mountain regional slope, respectively [ ].In this case, we assume that the surface topography, Z s (x), coincides with the water table, consistent Tóth's model (Tóth, 1963(Tóth, , 2009) ) where the water table is a subdued replica of the topography.From the surface of the lowest valley, the domain extends a distance B vertically [L].We assume B = 5 km with five valleys.The domain is bounded laterally by ∂Ω l and ∂Ω r , and at the bottom by ∂Ω b .
We also assume a decreasing permeability (κ z and κ x ; [L 2 ]) and porosity (ϕ; [ ]) with depth for both the alluvium and the mountain block.The only exception is the proximal piedmont alluvium (III), which is assumed to have a constant permeability higher than the rest of the systems, as suggested by Paola et al. (1992).For permeability in the remaining subdomains, a number of models have been used to describe this behavior (Ingebritsen & Manning, 2010;Jiang et al., 2009;Manning & Ingebritsen, 1999;Saar & Manga, 2004).Here, we use the piecewise permeability decay model presented in Saar and Manga (2004), which uses an exponential function for the shallow system and a potential function for the deep system.This model allows for a smaller decrease in permeability at lower depths, preventing singularities deep in the bedrock.
where κ s is the permeability at the surface, A s [m 1 ] and B d [ ] are the decaying exponents that represent the decreased rate of permeability with depth, d l is a defined depth (1 km for all the simulations) that separates the shallow from the deep systems, D d = 1 km, and A d can be calculated as follows, Water Resources Research Equation 2 calculates the permeability in the z direction, and we explore the effects of anisotropy in the shallow system with κ x = ηκ z , where η = 10 is the anisotropy value.We assume anisotropic behavior only occurs above d l in the upper alluvium (I ) and in the mountain block (VI) (see Figure 2b).The value of A s = d 1 l log(κ s ) log(κ d l ), where κ d l is the permeability at the distance d l that is extracted from Saar and Manga (2004) and Ingebritsen and Manning (2010).
For porosity, we assume κ z /κ s = ϕ/ϕ s ) m , where ϕ s is the porosity at the surface and m is a medium and processdependent coefficient.This expression has been used in earlier studies (see Bernabé et al., 2003;Cardenas & Jiang, 2010;Jiang et al., 2010).Replacing the permeability values with porosity in Equation 2 yields the expression below.The behavior of the porosity within the system is illustrated in Figure 2c.

Model for Fluid Flow and Transport of Heat and Solutes
The simulations presented in this work solve the following governing equations for groundwater flow and transport of heat and solutes (Bear, 1972;Fetter, 2001;Nield & Bejan, 2013), Beginning with the groundwater flow Equation 5, q is the Darcy's flux , n is an outward vector normal to the boundary [ ], ρ f and ρ f,0 are the water density and the water reference density [ML 3 ], respectively, μ f is the dynamic viscosity of water [ML 1 T 1 ], g is the gravity acceleration constant [LT 2 ], and p is pressure [ML 1 T 2 ].As an initial condition in the system, we use the hydrostatic equation.As boundary conditions, we use a constant head boundary on the left (∂Ω l ), a no-flow condition on the right (∂Ω r ), and bottom (∂Ω b ) boundaries, and a pressure boundary is set at the surface (∂Ω s ).
In the heat (T; [Θ]) transport Equation 6, is the surface temperature profile (at ∂Ω s ), Γ is the atmospheric lapse rate, C pf and C ps are the water and solid heat capacity, ρ s is the solid density and q Hb is the basal heat flux imposed along the bottom boundary (∂Ω b ).The system has a T profile as the initial condition, and a thermal insulation boundary condition was set in the lateral boundaries (∂Ω l , ∂Ω r ).
In the solute transport Equation 7, c is the volumetric concentration of salinity (NaCl) [ML 3 ], D = D ij is the dispersion-diffusion tensor [L 2 T 1 ] defined as (Bear, 1972) with α T and α L the transverse and longitudinal dispersivities [L], D m the effective molecular self-diffusion coefficient, ξ m = ϕ 1/3 is the fluid tortuosity (defined here with the Millington and Quirk (1961) model), and δ ij is the Kronecker delta function.In Equation 7, R = k mt (c max c) is a first-order reaction solute mass rate of dissolution of precipitation (solute source term), used in Lemieux et al. (2008a) and Provost et al. (2012) to mimic the buildup of salinity due to fluid-rock interactions.In this expression, c max is the maximum allowable fluid concentration, and k mt is the mass transfer reaction rate [T 1 ] that is temperature dependent and can be described through the Arrhenius expression (Langmuir, 1997), where A 0 is a temperature-independent coefficient [T 1 ], E 0 is the activation energy [L 2 MT 2 ] and R 0 is the ideal gas constant [ML 2 T 2 Θ 1 ].To obtain a reasonable value of A 0 we explore a range of values of k mt → k mt 0 used by Lemieux et al. (2008b) and Provost et al. (2012) and solve for A 0 in Equation 9 assuming a mean temperature of 160°C for those values (temperature at a depth of 5 km).It is worth noting that although this representation might not be appropriate for all environments, it allows for the consideration of realistic and applicable salinity scenarios.As an initial condition, we set the concentration to be equal to c max , because the system will flush out the excess of solute as the simulation progresses.As boundary conditions, we set a no flux condition on the right and bottom boundaries (∂Ω r and ∂Ω b ), a symmetry boundary condition on the left boundary (∂Ω l ), and an open boundary condition is set at the surface (∂Ω s ), with a reference concentration c 0 .
Lastly, we use empirical approximations for density and dynamic viscosity derived by Kestin et al. (1981), based on Rowe and Chou (1970), as they provide a reasonable approximation to experimental data (Adams & Bachu, 2002).The empirical expressions vary the density and viscosity according to variations in pressure (p), temperature (T ), and volumetric concentration of solutions (mainly NaCl) (c).These equations are presented in more detail in Adams and Bachu (2002) Equations 1 and 15.

Model for Residence Time Distribution
We model the residence time distribution (RTD) of water molecules Ψ(x, τ) [T 1 ] using the approach described in Gomez and Wilson (2013).From a numerical perspective, it is more stable to simulate the evolution of the where τ (τ ≥ 0) is age [T] and H(t) is the Heaviside step function.

Model for Bulk Electrical Resistivity
We estimate the electrical conductivity (the inverse of resistivity) for each simulation result using a modified version of Archie's Law, published by Glover et al. (2000).This modification is a mixing model that estimates the bulk electrical conductivity (σ r, bulk , S/m; [M 1 L 3 T 3 A 2 ]) of the medium with its porosity (ϕ) and both, fluid and solid, electrical conductivities (σ r,f , σ r,s ), as where m r is the cementation exponent (common values range from 1.5 to 2.5 for sedimentary rocks containing saline aqueous fluids) and p r = log(1 ϕ) m r / log(1 ϕ) .Observations show that the fluid's electrical conductivity primarily depends on salinity and temperature (Arps, 1953;Pepin, 2019;Sen & Goode, 1992;Ucok et al., 1980).Empirical models have been used to calculate fluid resistivity (ρ r,f (T, c) = 1/σ r,f (T, c)); some of the models use solute concentration and temperature (Pepin, 2019;Sen & Goode, 1992;Ucok et al., 1980), and others use only temperature with an assumed seawater salinity (Becker et al., 1982;Constable et al., 2009).We tested each model, comparing them with information on fluid resistivity collected by Pepin (2019) from reported measurements (Ho et al., 2000;Lide & Haynes, 2009;National Research Council, 1930;Noyes, 1907;Quist & Marshall, 1968;Ucok et al., 1980;Zimmerman et al., 1995).We decided to use the model presented in Sen and Goode (1992) for the comparison, as it provides a reasonable approximation to the data while having low model complexity (see the comparison in Figure S4 of the Supporting Information S1).The model form is the following, where M is the molality, and the temperature (T ) is in Celsius.
Glover's model (Equation 11) shows that as porosity decreases, solid resistivity becomes the dominant value in the system.However, the change in bulk resistivity depends on the water salinity.The transition between fluid and solid resistivity preponderance on the total bulk resistivity happens below a porosity of 5%, suggesting that we will be able to see saline water in systems with relatively low porosity (see Figures S5 and S6 in Supporting Information S1).
Lastly, we create a forward model using MARE2DEM (Key, 2016), a code that uses an Occam's inversion approach (Constable et al., 1987), to perform forward and inverse MT modeling.In that sense, we created forward models using the bulk resistivity (ρ r,bulk ) from the simulation results and populating the surface with three different arrangements of magnetotelluric (MT) receivers that explore 49 frequencies between 10 2 and 10 3 Hz (skin depth, δ, between 5 and 112 km with the lowest frequency and assuming resistivities of 1 and 500 Ω m, respectively, using the equation δ = 503 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ρ r,bulk /f √ from Chave & Jones (2012)).We then created synthetic noisy observations by adding a 5% relative Gaussian error.That is, a random error with mean zero and standard deviation of 5% of the simulated value for the variable of interest (bulk resistivity and phase) at each frequency.This Water Resources Research 10.1029/2023WR035906 means the error floor of our data is within the range of 5% for all of our synthetic observations, mimicking goodquality field data.We inverted both TE and TM modes as bulk resistivity and phase data to test whether an MT survey can capture the dynamics of the systems we are exploring.For the inversions, we used a telescopic meshing strategy that combines structured and unstructured meshes (see Figure S7 in Supporting Information S1).Within the mountain-to-valley system (x ∈ [ 51, 13] km and z ∈ [ 3, 1.6] km), we used a structured mesh with quadrilateral elements of size 0.1 × 0.05 km (width × depth).The idea is to use mesh with enough resolution to capture the complex spatial patterns expected from a nested regional groundwater flow system.For reference, the horizontal separation between the local topographic ridges is approximately 2.8 km.Then, enclosing the domain of interest, we use an unstructured mesh large enough to minimize artificial boundary effects (x ∈ [ 100, 100] km and z ∈ [ 100, 1.6] km).The final mesh has 68,273 elements (24,376 structured and 43,897 unstructured).Finally, we set a root-mean-square (RMS) misfit target of 1.0 for all our synthetic inversions.

Scenarios Explored
Although the conceptual model presented above is idealized, we strive to explore realistic scenarios that resemble shallow, intermediate, and deep circulation systems.With this in mind, we investigate two local and regional topographic systems, changing the regional slope (S m ) and the amplitude (A) of Equation 1.For permeability in the shallow region (depths larger or equal than d l ), we set the surface value (κ s ) and the value at d l (κ d l ), from this information, we back-calculate A s in Equation 2 using two values for κ d l that account for fast and slow decaying permeability cases and then use the anisotropy (η) to calculate κ x .In the deep region (depths larger or equal than d l ), we use two values for the deep decaying exponent (B d ) (the resulting permeability for both systems can be seen in Figure S8 of the Supporting Information S1).For solute concentration, we look at two reaction rates (k mt 0 ) that encompass some of the variability presented in Lemieux et al. (2008b).To estimate the resistivity fields, we use the porosity of the system, a constant value for the cementation factor (m r ), and explore constant values for solid resistivity (ρ r,s ), assuming distinct medium resistivity values for the system regions described in Figure 2. Lastly, we ran the simulations for t = 10 6 years to achieve steady-state conditions (less than 0.1% change), running a transient solution with the coupled model helps in a smooth convergence of the simulation.
The values for all parameters used in our analysis are presented in Table 1.We check the following physical requirements in all our simulations.(a) Recharge rates between 10 2 and 10 0 m/yr, as presented in other studies (Gleeson et al., 2011;Schaller & Fan, 2009;Wilson & Guan, 2004), (b) temperature profiles that follow realistic scenarios coming from previous simulations (Smith & Chapman, 1983), and (c) age distributions similar to the ones found in analogous systems (Figure 4) (Eastoe & Rodney, 2014;Frisbee, Phillips, et al., 2013;Frisbee, Wilson, et al., 2013;Mamer et al., 2014).

Role of Topography and Geology in Groundwater Dynamics and Residence Times
This section explores the topographic features and role of geology in groundwater system dynamics and residence times.In all our simulation results, the recharge rates in the mountain block systems vary between 0.35 and 0.9 m/ yr in the peaks of the mountains and the discharge rates from 0.1 to 0.6 m/yr in the valleys, with the highest recharge rates occurring in the high regional slope systems.These rates are comparable to values explored in other systems (Gleeson et al., 2011;Schaller & Fan, 2009).
First, we focus on the flow system with high local and regional relief systems with low permeability (panel g in Figures 3-5).In this mountain-to-valley transition system, we can see a distinct separation of local and regional flowpaths, where stagnation points are made visible by the lower Darcy velocities in the meeting branches of the local and regional flowpaths (panel g in Figure 3).The regional flow travels 0.5-2 km deep within the mountain block until it reaches the fault.In the fault, the water moves upwards in the first 400 m and then switches downward, recharging the alluvial piedmont and the lower alluvium.The recharge of the alluvial piedmont comes preferentially from the first local relief systems of the mountain block.This flow behavior strongly impacts the age distribution of this mountain-to-valley transition system (panel g in Figure 4), as water flowing through the local flowpaths ranges from less than 1 year to 100 years old, with increasing ages as depth increases.In the alluvium, the old regional groundwater flow mixes with younger water coming from the local flowpaths, increasing the water ages with increased distance from the piedmont into the alluvium, reaching ages older than Water Resources Research 10.1029/2023WR035906 1,000 years.Finally, the production of solutes is strongly related to the water age, as greater ages involve a longer contact time of water with minerals that could potentially dissolve, increasing the total dissolved solute (TDS) content of the groundwater.The solute concentrations presented for these systems (Figure 5) are divided into fresh (<1,000 mg/L), brackish (1,000-10,000 mg/L), and saline (>10,000 mg/L) water to make an easier comparison (with values from Stanton et al. (2017)).The solute concentration increases with depth, following the multi-scale nature of the flow field in the mountain block with brackish and saline water present in the stagnation points (panel g in Figure 5).In the alluvium, the solute distribution is a consequence of the interaction between the local and regional flow systems from the mountain block, with freshwater present at the alluvial piedmont and brackish water appearing nine km away from the mountain block, where young water mixes with older water.The appearance of brackish water on the alluvium surface is a consequence of the solute reaction rate, the interaction of the mountain flowpaths, and the change in permeability between the alluvial piedmont and the upper alluvium sections.The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The black lines show uniformly spaced streamlines.The streamlines in the mountain block show the distinct separation between the local and regional flowpaths recharging the lowland aquifers.The flow through the fault is preferentially downward in all the simulations.

10.1029/2023WR035906
Reducing the local topographic relief of the mountain block (panel e in Figures 3-5) produces shallower regional flowpaths with depths of 200 m, where the stagnation points disappear.In these systems, the regional groundwater flow travels at a higher rate through the mountain block, recharging the lowland aquifers through the alluvial piedmont and the lower alluvium.The water age increases linearly with depth, following the regional groundwater streamlines, and freshwater is present in the first kilometer of the mountain block due to the fastmoving regional flow.On the other hand, reducing the regional topography (panel c in Figures 3-5) leads to systems with enhanced local flowpaths that can reach 2 km deep.In these cases, the upper alluvium is predominantly recharged by the first peak of the mountain system, leaving the regional flow system to recharge the deeper parts of the lower alluvium.This behavior influences the water ages as younger water penetrates deeper into the mountain block, moving the stagnation points deeper where the concentration of solutes increases.Lastly, increasing the permeability of the mountain block (panels b, d, f, and h in Figures 3-5) produces faster and shallower regional groundwater fluxes, reducing the residence times of the overall mountain block and increasing the amount of freshwater present in the system.Furthermore, the recharge of water in the alluvial piedmont and The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The colors separate water ages.Each color denotes the water that has 60% of water in that specific range.In the mountain block, the water ages increase with increased depth.In the alluvium, the water ages combine young water from local flowpaths and old water from regional sources.

Water Resources Research
10.1029/2023WR035906 upper alluvium shifts to a predominantly regional source, which increases the amount of brackish groundwater present in the upper alluvium.In addition to the presented results, we also estimated the temperature distribution (Figure S9 in Supporting Information S1).The resulting temperature fields follow the multi-scale nature of flow in the mountain block and increase with depth, agreeing with results presented in previous studies (Gleeson & Manning, 2008;Jiang et al., 2009;Smith & Chapman, 1983;Tóth, 1963).
The latter simulation results explore a range of probable geologic and topographic features.The resulting behavior of the high local and regional relief systems agrees with age dating and geochemical analysis measurements performed in the Tularosa Basin (Figure S2 in Supporting Information S1).For instance, Carbon-14 age estimates show water ages between 1,340 and 27,500 years in the alluvium, increasing with increased distance from the mountain block (Eastoe & Rodney, 2014;Mamer et al., 2014).Similarly, TDS measurements in the area show freshwater present in the alluvial piedmont and brackish and saline water in sections of the basin-floor alluvium further away from the mountain system, attributed to the dissolution of gypsum (Newton & Land, 2016; ).The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The contours separate freshwater, brackish water, and saline water.The solute concentration increases with depth, delineating the multi-scale nature of flow in the mountain block.In the alluvium, the interaction of local and regional flowpaths produces the presence of brackish water in the upper alluvium.NMBGMR et al., 2016;Orr & Myers, 1986).However, not all mountain-to-valley transition systems have a high accumulation of solutes close to the surface.By decreasing the reaction rate of one of the systems, freshwater penetrates deeper into the mountain block, causing freshwater to extend into the upper part of the alluvium (Figure 6).In these systems, the downward movement of water through the fault is apparent, recharging fresh and brackish water to deeper parts of the alluvial valley.To illustrate the importance of multiscale circulation within the mountain block, we estimated water, solute, and energy fluxes passing through the fault from the mountain block to the alluvium (Figure S10 in Supporting Information S1).Water flux decreases with depth due to decreasing permeability (Figures S10a and S10b in Supporting Information S1).For systems with high reaction rates (red and blue lines), the solute fluxes are characterized by two peaks.The shallowest one is associated with local and intermediate circulation cells and occurs at approximately 50 m.This peak is more subtle for the cases with low permeability (blue lines).The deepest peak is linked to regional groundwater flowpaths that move brackish and brine water from the mountain block to the alluvium and occur at approximately 500 and 750 m for the low-and high-permeability cases, respectively.Similarly, the peak flux is lower for the low-permeability cases.After the second peak, we see a decrease of about an order of magnitude as we move to the deeper parts of the mountain block.However, these deeper contributions of solute mass still enter the fault and are eventually redistributed to the alluvium.For low reaction rate systems (green lines), we see a sharp decrease in solute flux with a shallow minimum at around 450 m.Because these systems do not have enough time to accumulate significant amounts of solutes in the upper regions of the mountain block (Figure 6b), we see freshwater entering the lower sections of the mountain block that are then moved through local and regional flowpaths to the alluvium.Lastly, the heat flux (Figure S10d in Supporting Information S1) only shows strong variations in the first 200 m (associated with local flowpaths), with heat fluxes varying between 0.05 and 0.25 W/m 2 .The lowest values are associated with low regional slope and low relief systems, whereas the highest values correspond to high regional slope and low relief systems.As we move deeper into the domain, these fluxes decrease to 10 2 W/m 2 .
In summary, the multi-scaled nested flow in mountain systems strongly depends on the topography and geology of the mountain block, where regional flowpaths drive water and solutes further into the alluvium in systems with high local and regional relief and high permeability.The interaction of local and regional flowpaths produces stagnation points in the mountain block.These points are delineated by brackish and saline water in systems with high dissolution rates, analogous to mountains with geologic material that can be easily dissolved, like evaporites.The following section explores how resistivity patterns give insight into the multi-scale nested nature of flow in the mountain block and how magnetotelluric (MT) measurements can be used to study this behavior.

Can We Image the Nested Nature of Hydrologic Systems Using Electromagnetic Methods?
As suggested by Glover et al. (2000), the resistivity patterns are a mixture that depends on fluid conductivity and lithology linked through the porous media (Equation 11).We also see that fluid conductivity is linked to heat and solute concentration, where systems with high temperature and solute concentration tend to be more conductive and, thus, less resistive.Our results show a clear distinction in the resistivity values of fresh and saline water in parts of the system with high porosity values (Figure 7).These patterns vary with permeability, where systems with low permeability and porosity are associated with increasing resistivity with depth in the mountain block.As a result, we see pockets of low resistivity close to the stagnation zones, where regional circulation meets the local flowpaths accumulating solutes.These zones appear below the valleys of the mountain system.On the other hand, ).The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The contours separate freshwater, brackish water, and saline.The resistivity fields show the location of stagnation points present in the mountain block, separating the local and regional flowpaths.
Water Resources Research 10.1029/2023WR035906 high permeability systems are linked to a slower decaying porosity.The results show a distinct decrease in resistivity as water becomes saline, starting at the interface between the local and regional flowpaths in the mountain block.Lastly, the bulk resistivity in the alluvium is primarily conductive, with resistivity values lower than 100 Ω m.This behavior is linked to a higher porosity, the presence of saline and brackish water, and a lower solid resistivity in the upper part of the alluvium.However, we also see contrasting resistivity values in the alluvial piedmont, where freshwater is being recharged from the mountain block.
As we mentioned earlier, the accumulation of solutes is linked to dissolution rates.If we reduce the dissolution rates of the system, we reduce the amount of solute buildup close to the stagnation zones (Figure 8).Thus, we also lose our ability to see notable changes in the bulk resistivity of the mountain block and the alluvium, especially in systems with low permeability, where most resistivity values are associated with the solid phase.These results also show the importance of geology, as systems with clay, shales, or other conductive materials could be mistaken for brackish or saline water when the resistivity is larger than 10 Ω m.
Apparent field electrical resistivity measurements in the Dosit River in China (Jiang et al., 2014) and the Tularosa Basin in New Mexico, USA (Newton & Allen, 2014;Newton & Land, 2016;Orr & Myers, 1986;Zohdy et al., 1969) have shown similar results.In these studies, the authors detected brackish and saline water in areas with low resistivity values and freshwater in areas with high resistivity values.It is worth remembering that the electrical resistivity strongly depends on the solid phase resistivity, porosity, and chemical composition of water (see Figures S5 and S6 in Supporting Information S1 illustrate the tight connection between the bulk resistivity and the porous media and the fluid resistivities).Our results suggest that geophysical methods, especially electromagnetic surveys, can serve as a tool for imaging mountain flow separation, where local and regional flow paths can be distinguished given the right solute concentration and geologic conditions.
Although these results are promising, there have not been, to our knowledge, electromagnetic measurements done in mountain systems that address the behavior we see in our simulations, partly due to the difficulty in performing such observations in these settings.Another concern is whether MT field data, which involve electromagnetic diffusion, can resolve the complex resistivity patterns shown in predictions in Figures 7 and 8. To that end, we designed a series of experiments on one of the groundwater flow and transport synthetic models, where we performed forward and synthetic data inverse MT simulations for variable receiver spacing densities using MARE2DEM.The main goal of this synthetic inversion approach is to test if MT data can detect and image the resistivity patterns found in the system, especially the Tóthian patterns found in the mountain block.
The forward and inverse models were tested on the model with low permeability and high topographic features (Figure 7g).We first performed an inversion using only the solid resistivity (porous media without fluids such as water or air) to see if we get back the general behavior of the system.Next, we changed the number of receivers in the system, first spacing the stations every 200 m, then we used 26 receivers spaced ∼1.5 km, and lastly, we used the same 26 receivers placed on top of structures we are interested in mapping in the mountain block.In our inversions, we used a unique initial condition of 1 Ω m.This strategy was reasonable given the good performance of the inversion results, as illustrated by all models reaching the target RMS misfit of one (see Figure S11 in Supporting Information S1).In addition, as described in the methods section, we used a telescopic modeling mesh that extends to a depth of 100 km.The idea is to avoid artificial boundary effects by refining the mesh in the area of interest (Figure S7 in Supporting Information S1).Furthermore, this depth selection is consistent with the skin depth range of 5-112 km estimated for our analysis.The oversampled inversion can capture the general behavior of the system; however, as expected with this type of inversion problem, the recovered model presents some smearing due to limited resolution (Aster et al., 2016); see, for example, the resistivity patterns in the alluvium (Figure 9).When fluid is considered in the resistivity calculation, the resistivity values in the alluvial valley change considerably due to a combination of higher porosity and saline water.The oversampled inversion, in this case, can capture the overall changes in resistivity in the alluvium and the location of the small pockets of low resistivity in the mountain block close to the stagnation zones.The inversion also recovers part of the fault's low resistivity values.The inversions show a higher resistivity than the true models at the bottom of the system.These overestimated values can be attributed to a plurality of reasons that include the extent of the MARE2DEM model, the initial conditions for the inversions, the non-uniqueness in the inversion arising from limitations in the inherent resolution of MT data, the Gaussian noise added to the data, along with the impact of the smoothing operator used to stabilize the inversion process.Decreasing the number of receivers (and hence increasing the receiver spacing) hinders the ability of the inversion to delineate some of the low resistivity pockets in the mountain block.However, placing sensors directly over these pockets allows us to characterize some of the structures of interest within the mountain block (Figure 9h).This modeling approach, coupled with optimization methods for sensor placement design (Hu et al., 2017;Robertson & Thiel, 2019) and expert advice, can give the best initial location of receivers in future studies.
These results show that an MT survey will be able to capture the overall behavior of the resistivity patterns in the investigated systems, especially the Tóthian patterns in the mountain block.Furthermore, the use of this method will allow us to see the pockets of saline water produced by the interaction of local, intermediate, and regional flowpaths, thus giving us a sense of how deep local and regional circulation occurs in mountain systems.

Dosit River Resistivity
As an additional analysis, we extracted the N-S cross-sectional area of the Dosit River, China, explored by Jiang et al. (2014), and used the coupled 2D groundwater flow and transport model applied in our idealized simulations to model the bulk resistivity of the system.This area was selected because Jiang et al. (2014) previously performed an MT survey of the region covering 38 km of the cross-section that can be used as a qualitative comparison for the bulk resistivity results.
The model was built using a smoothed cross-section of a Digital Elevation Model (DEM) downloaded from the SRTM data set of the region (Earth Resources Observation And Science (EROS) Center, 2017).We also used a decaying permeability and porosity with depth, using a surface permeability, κ s = 5 × 10 15 m 2 , a permeability at d l , κ d l = 1 × 10 16 m 2 , a decaying permeability exponent, B d = 3.2, a surface porosity, ϕ s = 0.3, and a first-order reaction rate, k mt 0 = 1 × 10 13 .The cross-section results show an accumulation of solutes near the river discharge, where the groundwater flow coming from both mountain systems converges (see Figure 10).This system also separates fresh, brackish, and saline water and gives a contrasting feature on the electrical resistivity patterns, similar to the patterns observed in the idealized case scenarios.

Water Resources Research
10.1029/2023WR035906 Measurements of total dissolved solutes at the surface of the river give 2,037 mg/L, other measurements near the area show values of concentration higher than a 1,000 mg/L below 150 m deep and values up to 8,499 mg/L at 1,400 m deep near the river basement (Jiang et al., 2014), similar to the results of concentration calculated in our simulations.Although the model used does not include a detailed lithology of the site, and the resistivity values obtained are higher than those presented in the article, the resistivity patterns are similar to those obtained by the authors (Figure S1 in Supporting Information S1).The system shows low resistivity near the river discharge and in the basement, where brackish and saline water is accumulating, and higher resistivity values where freshwater is located.These results reinforce the idea that the low resistivity values present at the system's bottom are likely attributed to variations in the groundwater solute concentration, as the authors also mention.

Discussion
Our simulations follow the nested nature of hydrologic systems expected in topographic-driven flow through mountainous systems (as presented in Tóth ( 2009)), with water age and brackish groundwater patterns within the alluvium comparable to the ones found in the Tularosa Basin (Eastoe & Rodney, 2014;Mamer et al., 2014;Newton & Land, 2016).This range of patterns in the alluvium is a consequence of the mixing between local and regional flowpaths in the mountain block that are sensible to the topographic and geologic characteristics of the mountain system.In all the simulations explored, we can see the importance of the regional groundwater flow in providing water and solutes to the alluvial aquifers and delineating the depth of local flowpaths in the mountain block.These regional flowpaths reach the alluvium with a high variation in depth, from hundreds of meters to kilometers deep, as other studies have suggested (e.g., Barroll & Reiter, 1990;Carroll et al., 2020;Frisbee et al., 2017;Mailloux et al., 1999;McIntosh & Ferguson, 2021;Pepin et al., 2014).
Our results also suggest that the water chemistry and lithological and geological conditions present in the mountain block provided the best conditions for the use of magnetotelluric surveys to map the nested groundwater circulation patterns.Mapping these patterns is crucial for characterizing the depth and the chemical and weathering processes occurring in the critical zone (Ackerer et al., 2021;Anderson et al., 2007;Condon et al., 2020;Holbrook et al., 2014;Riebe et al., 2017), and can also provide estimates of the mountain block recharge (Markovich et al., 2019;Wilson & Guan, 2004).Furthermore, inverse modeling of MT surveys, as presented here, could help to constrain the permeability structure of the crystalline basement to depths of up to 10 km, assuming that the survey acquired reliable information at depth and proper hydrogeophysical inversions methods are implemented (e.g., Herckenrath et al., 2012;Hinnell et al., 2010;Rubin & Hubbard, 2005).This is an important opportunity as there is a dearth of deep (>1 km) subsurface permeability information (Carroll et al., 2020;Ingebritsen & Manning, 2010;Manning & Ingebritsen, 1999).The available deep borehole data shows a wide standard deviation for metamorphic and igneous rocks forming the crystalline basement (Stober & Bucher, 2007).We expect a lateral variation in crystalline basement formation resistivity between the recharge and discharge areas as salinity builds up along the flow path at depths due to water-rock interactions (specifically fluid inclusion dissolution; Nordstrom, Ball, et al., 1989;Nordstrom, Lindblom, et al., 1989).
As presented in this study, the contrasting resistivity pattern comes from the interaction between fluid chemistry and the subsurface matrix propertied (porosity and resistivity), where brackish and saline water is less resistive.Mapping these "unconventional" water resources is critical in arid and semi-arid environments where the water supply is inadequate to meet the ever-increasing demand (Kang & Jackson, 2016;Qadir et al., 2007;Stanton et al., 2017;Viviroli et al., 2007).Stanton et al. (2017) note that there is about 300 × 10 3 km 3 of brackish water within the United States.As shallow, young groundwater resources become increasingly depleted, brackish water appears to represent an essential resource for municipalities.In Israel and Spain, desalinated brackish water is currently being used in the production of high-value crops (Aparicio et al., 2017;Ghermandi & Minich, 2017).Similarly, the Tularosa Basin has been proposed as a potential source of water for Las Cruces and Alamogordo, NM, and hosts the U.S. Bureau of Reclamation's National Inland Desalination Research Center, a critical testbed to assess the potential of brackish aquifers as an "unconventional" resource (Bureau of Reclamation, 2022;Newton & Land, 2016;Stephens, 2015).Magnetotelluric field campaigns may be of benefit in identifying deep subsurface brackish water resources within unconventional aquifer systems near the margins of basins.
Our study focuses on an idealized two-dimensional system without exogenous factors affecting the MT surveys.Field MT applications, however, must consider potential complicating factors such as 3D electrical structure and potential distortions arising from any strong near-surface gradients in conductivity (e.g., Ledo, 2006;Ledo et al., 1998).We refer the reader to Chapter 6 in Chave & Jones (2012) for a detailed look at these complicating factors in interpreting MT field data and strategies for deployment.Furthermore, additional considerations should account for the cation exchange capacity (CEC) of clays that can play a crucial role in the electrical conduction of shales and shaly sands (Schön, 2015), using other formulations of Archies' law, like the one presented by Waxman and Smits (1968), as found in other studies (Armadillo et al., 2020;Rizzello et al., 2022).

Conclusions
This work explores the influence of topographic and geologic variations on active groundwater circulation of mountain-to-valley transition systems and discusses optimum electromagnetic geophysical survey designs needed to detect the resultant flow paths.Our simulations show that groundwater ages and solute concentration patterns strongly depend on topography, geological and lithological structures, and dissolution rates.In the mountain block, solutes accumulate near the stagnation zones located in the meeting branches of local, intermediate, and regional flowpaths close to the valleys of the mountain.These solutes are then transported to the alluvial valley, where brackish water is diffused close to the surface in systems with high dissolution rates.This behavior is also observed in the distribution of groundwater ages in the alluvium, where a combination of young water from local flowpaths and old water from regional mountain circulation converged in the simulations.Although our systems are idealized, similar behavior can be seen in the Tularosa Basin in New Mexico, where young freshwater from the Sacramento Mountains system follows local and regional circulation patterns, increasing the water ages and dissolving gypsum that increases the total dissolved solids content in the alluvium (see Figure S2 in Supporting Information S1 and Mamer et al. (2014) and Newton and Land (2016)).These results show the importance of characterizing deep circulation in mountain systems, as water movement through relatively deep systems (>500 m) transports heat and solutes close to the alluvial valley surface.
The distribution of solutes in the mountain block and the alluvium strongly affects the bulk resistivity patterns.Contrasting resistivity values occur in parts of the system where there is a change from freshwater to saline and where porosity becomes low.These changes in resistivity occur near the stagnation zones in the mountain block, allowing us to map the depth where regional circulation meets the local flowpaths.We hypothesize that rigorous inversion techniques that couple electromagnetic geophysical surveys and groundwater flow and transport models could allow us to isolate the effects of nested flow processes in mountain systems.These techniques can also be used to estimate hydraulic parameters and state variables of these systems where information is scarce (Carroll et al., 2020).Note that this contrast of resistivity values can only be detected in systems with considerably high dissolution rates.
Lastly, we perform an idealized magnetotelluric inversion in one of the resulting scenarios to determine if electromagnetic surveys can detect resistivity patterns.Our results show that realistic low-noise MT data and 2D inversion can capture the depth and overall extent of the low resistivity structure in the mountain block.Our results also show that even sparsely spaced receivers (1.5 km spacing) can recover the low resistivity zones when carefully positioned with respect to local topographic high and low points.Our results give us confidence that carefully planned and executed MT measurements would be able to capture these structures in real systems.It is worth mentioning that the complexity of natural mountain systems can produce 3D circulation patterns that are not captured by our idealized cross-sectional scenarios.Future modeling and field studies should explore the potential implications three-dimensional topographic structure on the characterization of subsurface resistivity patterns with MT techniques.This research was supported by the National Science Foundation (NSF) (awards EAR-1830172, OIA-2020814, and OIA-2312326)  This manuscript has been co-authored by staff from UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE).The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes.DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http:// energy.gov/downloads/doe-public-accessplan).

Figure 1 .
Figure 1.(a) AfterBinley et al. (2015), comparison of horizontal and vertical survey scales for electrical resistivity tomography (ERT) with 1, 5, and 10 m electrode array spacing(Bernard, 2003), frequency domain electromagnetics (FDEM) with 3.5 m coil instrument, time domain electromagnetics (TDEM) with 50 and 500 m loop soundings at multiple stations along a transect(Bear et al., 2013), and magnetotellurics (MT) with 1 km spacing between stations.MT can survey deeper than 2 km, and its investigation depth depends on the conductivity, duration and sampling rate of the deployment, and the upper and lower frequency limits(Chave & Jones, 2012;Spies, 1989).(b) Histogram for groundwater borehole logging depth in the U.S.(Qi & Harris, 2017; Stanton et al., 2017).

Figure 2 .
Figure 2. (a) Conceptual model of the mountain-to-valley transition.The mountain system is a subdued replica of the water table described by Z s (x), and a fault characterizes the transition from mountain to valley at the base of the mountain range.The Roman numbers separate the system regions: (I) upper basin-floor alluvium, (II) lower alluvium, (III) upper proximalpiedmont alluvium, and (IV) mountain block.(b) Permeability variations with depth, the red and blue lines denote the permeability variation in the z and x direction, respectively.(c) Porosity variations with depth.The solid and dashed lines in (b) and (c) represent the two permeability scenarios explored for the mountain block (region IV).The dotted lines represent the permeability and porosity used for the alluvium (region II).

Figure 3 .
Figure 3. Spatial variability of the Darcy flux magnitude.The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The black lines show uniformly spaced streamlines.The streamlines in the mountain block show the distinct separation between the local and regional flowpaths recharging the lowland aquifers.The flow through the fault is preferentially downward in all the simulations.

Figure 4 .
Figure 4. Spatial variability of groundwater age.The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The colors separate water ages.Each color denotes the water that has 60% of water in that specific range.In the mountain block, the water ages increase with increased depth.In the alluvium, the water ages combine young water from local flowpaths and old water from regional sources.

Figure 5 .
Figure 5. Solute concentration for a high reaction rate k mt 0 = 1 × 10 13).The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The contours separate freshwater, brackish water, and saline water.The solute concentration increases with depth, delineating the multi-scale nature of flow in the mountain block.In the alluvium, the interaction of local and regional flowpaths produces the presence of brackish water in the upper alluvium.

Figure 6 .
Figure 6.Solute concentration for low permeability (a and c) and high permeability (b and d), and high dissolution rate (a and b) and low dissolution rate (c and d) scenarios.The contours separate water between fresh, brackish, and saline.Lower dissolution rates cause freshwater to penetrate deeper into the mountain block and alluvium.

Figure 7 .
Figure 7. Bulk resistivity predicted from the model results using Equation11for a high reaction rate k mt 0 = 1 × 10 13 ).The first column (panels a, c, e, and g) denotes scenarios with low permeability, and the second column (panels b, d, f, and h) denotes scenarios with high permeability.The rows present different degrees of topographic relief, starting with low relief and low regional slope (panels a and b), high relief and low regional slope (panels c and d), low relief and high regional slope (panels e and f), and high relief and high regional slope (panels g and h).The contours separate freshwater, brackish water, and saline.The resistivity fields show the location of stagnation points present in the mountain block, separating the local and regional flowpaths.

Figure 8 .
Figure 8. Bulk resistivity predicted from the model results using Equation 11 for low permeability (a and c) and high permeability (b and d), and high dissolution rate (a and b) and low dissolution rate (c and d) scenarios.The contours separate freshwater, brackish water, and saline.Lower dissolution rates cause the loss of contrasting resistivity values within the mountain block.

Figure 9 .
Figure 9. Synthetic inversion tests assessing the ability of MT data to image the resistivity of the groundwater system features.The columns denote the forward (a, c, e, and g) and inverse (b, d, f, and h) models.In (a) and (b) the inversion of the idealized system is presented with receivers spacing every 200 m.The following three rows denote the same system with receivers spacing every 200 m (c and d), ∼1.5 km (e and f), and 26 receivers placed on top of structures of interest (g and h).The inverted white triangles show the location of the receivers in each inversion.The low resistivity pockets in the mountain block are captured with good sensor placement.

Figure 10 .
Figure 10.Simulations of the Dosit River transect.The panels show the results of solute concentration (a), residence times (b), and bulk resistivity (c).The contours denote the separation of fresh, brackish, and saline groundwater.The vertical dashed lines show the area where Jiang et al. (2014) performed their MT survey (Figure S1 in Supporting Information S1).
and the AGU Horton Research Grant awarded to D. Gonzalez-Duque.Additional support was provided by the U.S. Department of Energy, Office of Science, Biological and Environmental Research.This work is a product of two programs: (a) Environmental System Science Program, as part of the Watershed Dynamics and Evolution (WaDE) Science Focus Area at Oak Ridge National Laboratory and the IDEAS-Watersheds project, and (b) Data Management Program, as part of the ExaSheds project.