The southern Mesilla Basin, located within western Texas and southern New Mexico, is one of two major sources of groundwater for the city of El Paso, Texas, and provides ∼30% of the region’s domestic groundwater needs. Groundwater is also used for agriculture in the non-urbanized regions of the basin. The basin is one of the southernmost basins of the Rio Grande rift where extension has overprinted older features, including extensive Eocene trachyandesitic volcanism and Laramide deformation. An increase in groundwater salinity is observed from north to south within the basin. Some previous researchers have suggested that this salinity change is due to runoff and recharge from agricultural activity. We use a combination of gravity and groundwater geochemistry in an integrative study to determine the possible influences of faults and other subsurface structures on groundwater salinity and quality. Gravity studies suggest the presence of other fault systems within the basin that serve as conduits for deeper, warmer, more Si-rich waters in the northern part of our study area, and as recharge zones for Ca-rich surface runoff from the carbonates of the Franklin Mountains in the eastern portion of our study area. The high Cl/Br ratios found in 90% of wells suggest that the salinity increase is primarily due to dissolution of evaporites within the basin rather than to deep-basin brines. Agricultural activity and water interaction with igneous and carbonate bedrock also are minor influences. This study highlights the unique resolving power of combining geophysical and geochemical techniques in understanding groundwater chemistry and its relationships with local geological structures. Such an approach can be readily applied to other systems with similar geologic and hydrological settings for groundwater exploration and resource management.

The Mesilla Basin is one of the southernmost basins of Rio Grande rift system. The basin extends south from ∼32.5°N in New Mexico to the U.S.–Mexico border, with a surface area of ∼2900 km2 (Fig. 1). The basin serves as an important source of groundwater for the entire region—which encompasses the cities of Las Cruces, New Mexico, and El Paso, Texas, with a combined population of >1,050,000 (U.S. Census Bureau, 2015)—including for agricultural irrigation for ∼720 km2 of agricultural land, especially in times of drought when availability of Rio Grande water is limited. Our study area focuses on the southeastern portion of the basin, which forms the present valley of the Rio Grande and extends ∼27 km from near the Texas–New Mexico border at 32°N to the U.S.–Mexico border. It is bounded to the east by the Franklin Mountains and to the west by La Mesa surface (Fig. 1).

Although the Rio Grande valley is currently the topographically lowest part of the Mesilla Basin, the thickest valley fill (∼1100 m) is found ∼10 km west of the Rio Grande (Hawley and Kennedy, 2004). The basin contains a mixture of alluvial, fluvial, and igneous material. The southern boundary of the basin is controlled in part by the Eocene-age Cerro de Cristo Rey intrusion (Fig. 1). A number of smaller intrusions that are geochemically similar to Cerro de Cristo Rey (Hoffer, 1970; Barnes et al., 1991) are found throughout the southern basin (Fig. 1). The northern edge of the Jurassic Chihuahua trough (rift) that affected subsequent Laramide deformation also appears to have been located in the southern Mesilla Basin (Lawton, 2004; Seager, 2004).

Groundwater within the basin primarily flows to the southeast and exits through the Paso del Norte (Hawley and Kennedy, 2004) near the point where the U.S. states of New Mexico and Texas and the Mexican state of Chihuahua meet. Groundwater salinity within the basin gradually increases toward the south (Gelhar and McLin, 1979; Hibbs and Merino, 2007). The Canutillo well field supplies groundwater for the city of El Paso (Fig. 1). The field supplies between 31,000 ML/yr (megaliters per year) during normal years and up to 43,000 ML/yr during drought conditions when no water is available from the Rio Grande (El Paso Water Utilities, 2007). This is ∼30% of the city’s total groundwater needs. Numerous private wells for agricultural and residential water use are also located in the Mesilla Basin.

Several major and numerous minor faults have been mapped within our study area (e.g., Witcher et al., 2004; Hawley and Kennedy, 2004; Khatun et al., 2007) (Fig. 1), but the relationship between these faults and groundwater quality is not well understood. The age of the last movement on the faults within the Rio Grande valley is unknown, although faults on the La Mesa surface (Fig. 1) appear to have offset the surface 3–6 m in the past 750,000 yr (Machette et al., 2000).

The goal of this study was to determine how subsurface geologic features influence water salinity and quality within the southern Mesilla Basin. First, we collected geophysical data (primarily gravity) to better determine the locations of faults, extent of igneous intrusions, and depth to bedrock within the study area. Second, we collected new geochemical data for 65 water wells within the basin, as well as compiled extensive water-quality information from previous groundwater studies. These data allowed us to comprehensively analyze changes in water chemistry and quality and their relation to surface and subsurface geological features.

Geological Setting

Normal faulting in the Mesilla Basin is related to extension within the Rio Grande rift that began in the Quaternary (Seager and Morgan, 1979) and has remained active but less intense in the past 2–3 m.y. (Seager et al., 1984). The rift has been subjected to two stages of deformation. The first stage of deformation began in the region by 24 Ma (Henry and Price, 1986) with ENE-WSW–oriented extension creating broad, northwest-trending grabens. A second pulse at 12–15 Ma produced north-south–trending basins (Keller and Cather, 1994; Langford et al., 1999), including the southern Mesilla Basin. Maximum fault displacement within the basin likely occurred from 4 to 10 Ma (Hawley and Kennedy, 2004). Infilling of the basin occurred from 4 to 0.7 Ma (Hawley and Kennedy, 2004). The northwestern portion of our study area contains ∼250 m of Cenozoic fill (Hawley and Kennedy, 2004), shallowing to 20–80 m of fill at the southern edge of the study area.

Basin fill in the study area consists of the Oligocene to lower Pleistocene Santa Fe Group and Pleistocene to Holocene Rio Grande alluvium (Hawley and Kennedy, 2004). The Rio Grande alluvium is estimated to have an average thickness of ∼46 m (Leggat et al., 1963), and the Santa Fe Group a thickness of 460 to 760 m (Hawley and Kennedy, 2004; Witcher et al., 2004). The lower Santa Fe unit formed prior to uplift of the present-day mountains and consists primarily of playa deposits and eolian sands. Coarse-grained material is found only at the edges of the basin (Hawley and Kennedy, 2004). The middle Santa Fe unit was deposited when uplift of the surrounding mountain ranges began and consists of intertonguing of alluvial and basin-floor materials (Hawley and Kennedy, 2004). Gypsum-selenite and mirabilite-thenardite evaporites and calcic cements are found in both the lower and middle Santa Fe (Hawley and Kennedy, 2004). The ancestral Rio Grande was present in the region by 2–3 Ma (Hawley et al., 2009) and brought fluvial material from beyond the local basin, as observed in sediments of the upper Santa Fe unit (Sellepack, 2003), which lack lacustrine deposits (Hawley and Kennedy, 2004). The Rio Grande shifted from the east side of the Franklin Mountains to the west side of the mountains ca. 2 Ma, but did not begin to incise its present course until ca. 0.7 Ma (Mack et al., 2006). Thus, basin fill alternates both vertically and horizontally between playa–alluvial fan deposits, eolian deposits, and fluvial deposits.

The Mesilla Valley fault zone is the primary fault system within our study area. The location of the Mesilla Valley fault zone shown in Figure 1 is based on water-well information from Hawley and Kennedy (2004) and Witcher et al. (2004), but is also consistent with the results of Sellepack (2003) who used well logs, cuttings, and measured sections of outcrop of the Santa Fe units. There is as much as 250 m of displacement of the lower Santa Fe unit across the fault zone (Hawley and Kennedy, 2004). Sellepack (2003) suggested that at least two other north-south faults are located east of the Mesilla Valley fault zone, while gravity studies by Khatun et al. (2007) suggested that three faults (the River, I-10, and Three Sisters faults) are located east of the Mesilla Valley fault zone.

It is important to emphasize that recent rift features are superimposed on a number of older structures related to pre-Cenozoic deformation. Lawton (2004) placed our study area at the northern edge of the Chihuahua trough rift system in the Late Jurassic. Consequently, in the Laramide, the study area appears to have formed a transitional zone between the more “thin-skinned” thrust faulting found in northern Chihuahua (e.g., in the Sierra de Juárez located just south of our study area) and the higher-angle reverse faulting of southwestern New Mexico (e.g., Seager, 2004; Carciumaru, 2004; Scharman, 2006).

Eocene intrusions of similar age and composition (Barnes et al., 1991; Hoover et al., 1988; Hoffer, 1970; Garcia, 1970; Lovejoy, 1976) are found in a number of localities within the southern Mesilla Basin (Fig. 1), including the most prominent intrusion, Cerro de Cristo Rey, which appears to be a plug-like body formed from the fractional crystallization of subduction-related basaltic magmas at subcrustal depths (Barnes et al., 1991). Preliminary geophysical studies of some of the smaller intrusions (e.g., Baker et al., 2012; Montana et al., 2012) indicate that they are considerably more extensive at depth. These bodies have likely played a significant role in the formation of many recent features in the region, including stream and arroyo morphology, hydrology, and soils.

Other lithologic units that may influence surface and groundwater chemistry include Cretaceous shales and limestones found in outcrops surrounding the Eocene intrusions (Fig. 1) and the Paleozoic and Precambrian rocks of the western Franklin Mountains (Fig. 1). Table 1 provides a brief lithologic description of these units.

Hydrological Framework

Shallow wells in the Mesilla Basin produce water from the Rio Grande alluvium aquifer for domestic or agricultural purposes, but the water is highly variable in salinity (e.g., Leggat et al., 1963). Gelhar and McLin (1979) showed that total dissolved solids (TDS) in the alluvium aquifer increases from ∼1000 mg/L in the northern Mesilla Basin to 8000 mg/L in the southern basin. The alluvium aquifer is recharged by seepage from the Rio Grande and irrigation water (Sheng, 2013). Individual wells with high salinity or iron values appear to tap abandoned river channels and swamp or bog deposits (Arunshankar, 1993). The Rio Grande alluvium is hydrologically connected to the underlying aquifers within the Santa Fe Group, leading some researchers (e.g., Sheng, 2013) to propose that use of Rio Grande water and saline groundwater from the Rio Grande aquifer for flood irrigation has led to the increased salinity observed within the Santa Fe Group aquifers in the southern Mesilla Basin.

Due to its low water salinity (e.g., <1000 mg/L), the aquifer within the middle Santa Fe unit is the most heavily pumped aquifer for drinking-water, industrial, and irrigation use (Wilson and White, 1984). The Canutillo well field produces from the middle Santa Fe unit and from an eolian unit of the lower Santa Fe (Hawley and Kennedy, 2004). The water with the lowest TDS values in the Santa Fe units is generally observed from the Canutillo well field northward.

Permeability varies greatly in the Mesilla Basin, with horizontal permeability as much as ten times greater than vertical permeability (Witcher et al., 2004). Clay lenses restrict vertical flow (Nickerson, 1989), leading to an overall southeastern flow of groundwater from north to south that exits the basin at the Paso del Norte (Hawley and Kennedy, 2004). In the Santa Fe units, hydraulic conductivity ranges from <0.03 to 30.5 m/day (Witcher et al., 2004).

Witcher et al. (2004) suggested several sources for groundwater salinity within the basin, including dissolution of carbonate and gypsiferous sediments found in the middle and lower Santa Fe units. They also suggested that Paleozoic and Cretaceous rocks found at the edges of the basin could act as conduits for deeply circulating, saline groundwater. There is evidence for upward flow of water in deeper wells (>180 m) at the edge of the basin in the northeastern portion of our study area (Nickerson and Myers, 1993), and several artesian wells emanated from the deep aquifer prior to its extensive development (Leggat et al., 1963). Geochemical tracers (major ions, stable O and H isotopes, sulfur and strontium isotopes) have been used to indicate that the possible inter-basin groundwater flow from the Jornada del Muerto Basin to the Mesilla Basin through connected aquifers could be responsible for the salinity increase in the northeastern Mesilla Basin (Langman and Ellis, 2013). Finally, upward circulation of geothermal water along fault zones and fractures may also lead to increases in salinity, as suggested by the presence of geothermal (>26 °C) waters at several wells within the study area (Witcher et al., 2004; Nickerson, 2006).

Geophysical Surveys

Shallow structure (<1000 m depth) within the basin is reasonably well constrained by water-well logs (Hawley and Kennedy, 2004). The extensive urbanization, agricultural activities, and cultural noise (e.g., power lines, pipelines, fences) in our study area limit our ability to use geophysical techniques such as electrical soundings to image the deeper subsurface (see Arunshankar, 1993). We have found the gravity method to be one of the most effective techniques to resolve deeper structure in the study area, especially in the heavily urbanized regions south of 31°53′N.

Khatun et al. (2007) conducted a microgravity study of the structure of the Mesilla Basin using a LaCoste-Romberg Model G gravimeter with elevation control from differential GPS or leveling surveys from existing benchmarks. They collected ∼1200 data points with spacings of 60–200 m focusing on the central portion of our study area (Fig. 2). This study identified several potential north-south–trending faults within the Mesilla Basin located east of the Mesilla Valley fault, including one paralleling the present Rio Grande (River fault), one paralleling Interstate Highway 10 (I-10 fault), and one along the western edge of the Franklin Mountains (Three Sisters fault) (Fig. 1). The northern portion of the River fault and southern portion of the Three Sisters fault are consistent with the locations of faults suggested by Sellepack (2003) based on well log and outcrop studies.

Data from Khatun et al. (2007), Baker et al. (2012), and Montana et al. (2012), as well as from several unpublished studies, form part of the gravity database for western Texas and northern New Mexico that has been carefully compiled, processed, quality checked, and maintained by the University of Texas at El Paso (UTEP, (triangles, Fig. 2). We used these data in our analysis, and collected ∼400 new data points to fill gaps in the preexisting data set.

Gravity Data

We collected new gravity data using a Lacoste-Romberg Model G gravity meter tied to the absolute gravity base station on the UTEP campus. One survey, with station spacing of ∼100 m (western half of line Q-Q′, Fig. 2), was designed to extend from the La Mesa topographic surface across the Mesilla Valley fault zone to tie with previous surveys within the river valley. Other surveys were conducted to fill gaps in gravity data for the northwestern and southeastern portions of our study area. Each gravity data collection loop was completed in <4 h. Elevations were determined through differential GPS.

We corrected the gravity data for tidal effects, dial constant, and drift. We then corrected for latitude, free air, and terrain. Terrain corrections used digital elevation models from the U.S. Geological Survey (USGS). We used the North American Datum of 1983 (NAD83) for these corrections. Finally, we applied the complete Bouguer correction with a reduction density of 2670 kg/m3. Data from the UTEP gravity database were reprocessed using the NAD83 datum to allow merging with the new data. More details of the data reduction process have been described by Avila (2016) and Hiebing (2016).

Groundwater Sampling and Analysis

We collected 65 groundwater samples from irrigation, domestic, and municipal drinking-water wells throughout the study area (Fig. 3A). For large wells that are run continuously or for long periods of time, grab samples were obtained from the piping nearest to the wellhead. For wells that were not continuously pumped, they were run for 5–10 min until all of the monitored parameters were stabilized prior to sampling. Sample location, date, depth, pH, nitrate content, temperature, and TDS content were recorded with a YSI Professional Plus multimeter in the field. A table summarizing this information may be found in Table S11.

Samples were prepared for major element analysis by inductively coupled plasma–optical emission spectrometry (ICP-OES) and ion chromatography (IC) at UTEP. In the lab, water samples were filtered using a 0.45 µm cellulose acetate filter to remove sediment and particulates and placed in two 250 mL acid-washed high-density polyethylene Nalgene bottles. The contents of one bottle were acidified with three drops of concentrated nitric acid for major cation analysis, and the contents of the other bottle were archived without acidification for immediate anion analysis. For major cation concentrations (Ca, Mg, Na, K, and Si), ∼15 mL of acidified sample was used for analysis on a Perkin Elmer 5300DV ICP-OES system. The USGS M-210 and NIST 1640a standards were analyzed three to five times during each run to assess measurement precision. Percent error of the standards was no greater than 10%. For major and trace anion concentrations (Cl, SO4, Br, etc.), the non-acidified filtered sample was diluted with deionized water approximately ten times. The accurate dilution factor for each sample was calculated with sample weights. These samples were analyzed using a Dionex ICS-2100 IC system. An in-house water standard was measured at least twice during each run to ensure accuracy. In general, standard errors were no greater than 12%. Selected water samples were also analyzed for trace element concentrations (Fe, As, V, etc.) at the Pennsylvania State University. Complete details of the groundwater sampling and analysis procedures have been described by Hiebing (2016) and Garcia (2017).

River water samples collected by Nyachoti (2016) and Hiebing (2016) (Fig. 3A; Table S1 [footnote 1]) and analyzed using the same procedures at outlined above were used for comparison purposes with the groundwater information. In addition to groundwater data collected by Hiebing (2016), we compiled extensive geochemical information on groundwater from Witcher et al. (2004), Leggat et al. (1963), Nickerson (2006), the Texas Water Development Board’s (2016) database, and the U.S. Geological Survey (2016), which were used in our analyses (see Fig. 3B; Table S2 [footnote 1]).

Gravity Anomaly Interpretations

The complete Bouguer anomaly map for the study area was interpolated using minimum curvature with a grid size of 250 m (Fig. 2) and compared to locations of known existing faults in this region (orange dashed lines in Fig. 2 are faults inferred by Hawley and Kennedy [2004] and Witcher et al. [2004] from water-well log information; white dashed lines are faults inferred by Khatun et al. [2007] based on densely spaced gravity data). Note that there is a general north-south trend to anomalies within the region. Anomaly highs are associated with the Franklin Mountains, and lows with the northwestern portion of the study area where sediment is thicker. The highest anomaly values are observed in the southeastern portion of the study area where Eocene intrusions are found. From here, anomaly values decrease toward the northeastern portion of the study area (∼31°57′ N), although upper Paleozoic rocks crop out in this region (Fig. 1).

We obtained a Bouguer residual anomaly map (Fig. 4) for the study area by fitting a third-order polynomial to the complete Bouguer anomaly map for a larger region of southern New Mexico and western Texas and then subtracting the polynomial to eliminate longer-wavelength, deeper-seated regional variations in gravity. This residual map reveals a complicated series of anomaly highs and lows within the study area.

High residual anomaly values are associated with most outcrops of Eocene intrusions in the southeastern study area, with suggestions that some of the smaller outcrops (e.g., the River, Westerner, Three Sisters intrusions) may be linked to Cerro de Cristo Rey by feeder dikes. North of Transmountain Drive (∼31°57′N), a distinct anomaly low is observed. This region is covered with a thin veneer of alluvial material (see Fig. 1), but Hawley and Kennedy (2004) (their section J-J′) suggested that this region is underlain by Paleozoic rocks. South of this low is an exposed northwest-southeast–striking anticline of upper Paleozoic rocks, and north of the low is a north-south–striking syncline of the same exposed Paleozoic units (Hawley and Kennedy, 2004) (Fig. 1).

A residual high located near the town of Westway, Texas (∼31°58′N, 106°35′W), does not correspond to any known subsurface features. It is possible that this feature is a continuation of the anticlinal structure observed in outcrop to the southeast, or it could be a remnant of lower Paleozoic rocks that have been thrust over the upper Paleozoic rocks by an eastward-dipping fault. Cross sections by Carciumaru (2005) and Scharman (2006) located 5–6 km to the southeast of the anomaly high both show an eastward-dipping thrust fault along the mountain front.

Anomaly lows in the middle of the study area tend to follow the course of the Rio Grande and step to the southeast at the southern end of the basin (∼31°51′N). Several studies (e.g., Khatun et al., 2007; Sellepack, 2003) have suggested that an east-west–striking accommodation structure may exist in this region that transfers stress southeastward to the adjacent Hueco Basin fault system located east and south of the Franklin Mountains. The anomaly high in the southwestern portion of the study area reflects the shallowing southeastern edge of the Mesilla Basin, consistent with water-well data (Hawley and Kennedy, 2004).

Gravity data in the northwestern part of the study area (Fig. 4) are too sparse to pinpoint the location of the Mesilla Valley fault zone or Western fault both mapped by Witcher et al. (2004). To the southwest, the Mesilla Valley fault zone bends around the eastern edge of a gravity high. The River fault, as mapped by Khatun et al. (2007), appears to be delineated by a transition from low to moderate residual values in the central study area. The Three Sisters and I-10 faults, as mapped by Khatun et al. (2007), do not correlate well with trends observed in the residual values.

A horizontal gradient magnitude (HGM) map (Fig. 5) was constructed using methods outlined by Grauch and Johnston (2002). Abrupt changes across near-vertical features (such as faults or the edges of intrusions) produce the highest HGM values. The map shows a complicated pattern of small variations in gradient that likely reflects the fact that many changes in basement structure may not be cut by near-vertical boundaries. Distinct gradient highs are associated with the Cerro de Cristo Rey, Three Sisters, River, and Westerner intrusions, but not with several intrusions located east of the Three Sisters. There appear to be east-west–trending highs in the southern study area where we expect that an accommodation zone is located. A gradient high is also associated with the southern Canutillo well field. Lack of data for the Mesilla Valley fault zone does not allow us to draw conclusions about its structure. The HGM map suggests that the River fault may be composed of two separate branches in the northern part of the study area. Brown solid lines in Figure 5 indicate our new interpretations of fault locations based on the HGM and residual gravity maps. In most cases, we have shifted faults by <1 km from their locations in previous interpretations.

Forward Modeling of Gravity Profiles

We modeled the gravity anomaly data along four profiles (Fig. 2) using the GM-SYS gravity modeling software package (sold by Geosoft, based on the forward modeling techniques of Talwani et al. (1959). Three of the profiles (J-J′, K-K′, and NW-SE) correspond to portions of hydrologic cross sections published by Hawley and Kennedy (2004) (based on water-well logs and cuttings) that we have used to help constrain the upper portions (∼1000 m) of our density models. The fourth profile (Q-Q′) was constructed to obtain information on structure within the middle part of our study area and takes advantage of closely spaced data that we collected extending from the present valley floor to the La Mesa surface. The density models that we obtained for these profiles are shown in Figure 6. The geologic units used in the modeling and their corresponding densities are given in Tables 1 and 2. The densities are based on the gravity studies of Avila et al. (2016) and Imana (2002).

The southernmost profile, K-K′ (Fig. 6), was constructed based on Hawley and Kennedy’s hydrogeologic profile K-K′. Note that we have overlain a portion of Hawley and Kennedy’s profile (extending from 425 to 1200 m above sea level) on the upper part of our density model. Three faults are required to fit the observed gravity data. The Western fault appears to be a continuation of a fault that Witcher et al. (2004) showed as terminating ∼2.5 km north of K-K′. The second fault is the Mesilla Valley fault that Hawley and Kennedy (2004) showed as crossing their hydrogeologic profile. The third fault, labeled RF, is located ∼2.5 km east of the Mesilla Valley fault and could be the southern end of the River fault inferred by Khatun et al. (2007). We do not see evidence for the I-10 fault or Three Sisters fault in this profile; the probable positions of these faults are marked by bars. An intrusion is required to match the increase in gravity observed at ∼12 km along the profile. This may be the base of the intrusion associated with the Westerner outcrop (W in Fig. 2). The density model suggests that the lower Santa Fe unit is ∼500 m thick between 0 and 8 km distance along the strike of K-K′ (Fig. 6), ∼100 m thicker than previously estimated by Hawley and Kennedy (2004). The density model suggests that the Western fault has offset the basement by ∼100 m but has produced very little offset of the Santa Fe units. The main Mesilla Valley fault appears to offset the basement and lower Santa Fe unit by 300 m, while the River fault offsets the basement by ∼100 m.

Profile Q-Q′ crosses both the Western fault and the Mesilla Valley fault zone, which are required to match the observed gravity data (Fig. 6). Note that this profile also has changes in strike. Two faults (labeled RFW and RFE) are also required to the east of the Mesilla Valley fault zone but do not appear to significantly offset the upper Santa Fe unit. These are interpreted to represent the two strands of the River fault (Fig. 5). Although profile Q-Q′ does not correspond to a hydrogeologic cross section constructed by Hawley and Kennedy (2004), the sediment thicknesses are consistent with those shown to the south along their profile K-K′ and to the north along J-J′. The Western fault and Mesilla Valley fault zone appear to offset both the basement and the lower Santa Fe unit by ∼140 m. The two strands of the River fault appear to also offset basement by ∼140 m, but not the lower Santa Fe unit.

Profile J-J′ (Fig. 6) extends along a portion of hydrogeologic profile J-J′ of Hawley and Kennedy (2004) (their profile extends in depth from 320 to 1200 m above sea level). Water wells constrain the thickness of the Santa Fe units at several places along the profile. Note that the deepest part of the lower Santa Fe unit on the western side of our profile (between 0 and 3 km distance along J-J′, Fig. 6) is ∼800 m thick, ∼400 m thicker than shown by Hawley and Kennedy (2004), but they show no water well that penetrated to bedrock in this region (see Fig. 6). The Cretaceous unit along this profile is modeled as being thinner than to the south (consistent with the interpretation of Hawley and Kennedy [2004]), suggesting thinning northward across the limit of thin-skinned Laramide deformation. The Western fault, not shown on Hawley and Kennedy’s (2004) profile, is required to offset bedrock by ∼160 m to match the observed gravity data but does not appear to offset the Santa Fe units. The Mesilla Valley fault zone offsets bedrock by ∼400 m and the Santa Fe units by ∼200 m. Strands of the River fault offset bedrock by 80–160 m but do not appear to offset the Santa Fe units. The I-10 and Three Sisters faults appear to align with the edges of an inferred igneous intrusion. Profile J-J′ also crosses a gravity anomaly low (Figs. 2 and 4) observed near the Franklin Mountains at 31°57’N. The gravity data are consistent with a shallow (∼200 m) basin containing upper and/or middle Santa Fe units underlain by upper Paleozoic units.

The final two-dimensional gravity profile that we modeled extends from northwest to southeast (Fig. 2) along a major portion of hydrogeologic cross section NW-SE constructed by Hawley and Kennedy (2004) (their profile extends in depth from 300 to 1200 m above sea level). Note that the profile has several changes in strike (Fig. 2). The profile shows valley fill thinning from ∼750 m at the northwestern end to <90 m at the southeastern end, in agreement with Hawley and Kennedy’s cross section. The most significant decrease in fill thickness occurs across the Mesilla Valley fault zone. The density model includes a slight offset (∼100 m) in units along the Western fault and a total offset of 700–800 m along the Mesilla Valley fault zone. Neither the River fault or I-10 fault are required by the model to match the observed gravity, however this profile is located near the southern end of both faults where offsets would be expected to be low. The gravity data support the bedrock high and shallow “Sunland Paleo Valley” at the southern end of the study area (between 18.5 and 21.5 km distance along profile NW-SE, Fig. 6) as mapped by Hawley and Kennedy (2004). Unlike the Westerner outcrop (profile K-K′), the River outcrop, a small Eocene intrusion located on the west bank of the Rio Grande, does not appear to lie directly above a major subsurface intrusion.

Geochemical Information

In this section, we focus on how groundwater chemistry may be influenced by faulting and bedrock composition. In general, the TDS values of the 65 groundwater samples (with well depths from 30 m to ∼400 m) collected in this study range from ∼500 mg/L to ∼17,000 mg/L, averaging ∼1400 mg/L. A piper diagram showing major element concentrations of these samples (Fig. 7) indicates the presence of dominant Na-Cl and Ca-Na-SO4-Cl types of water in the aquifers of this study, consistent with previous studies in the Mesilla Basin (Witcher et al., 2004). Mineral saturation index calculations show that most of these groundwater samples are oversaturated with respect to carbonate minerals such as dolomite and calcite but undersaturated with evaporative minerals such as halite and gypsum.

River chemistry studies (Gaillardet et al., 1999) have suggested that surface-water chemistry is dominantly controlled by the lithology within the rivers’ catchments. We applied a similar classification (carbonates, silicates, and evaporites; Fig. 8) scheme to groundwater (green symbols in Fig. 8) and river water (blue symbols) collected in the study area by Hiebing (2016) and Nyachoti (2016). The ratios of major elements in these samples (Fig. 8) indicate that groundwater chemistry in this region is dominated by the dissolution of evaporite minerals (e.g., gypsum and halite) and interactions with aluminosilicate minerals that make up the majority of the Santa Fe units.

Note that several groundwater and most of the river samples are characterized by high Ca/Na and Mg/Na ratios (Fig. 8), which might indicate the dissolution of carbonate rocks. For example, Ca/Na ratios in groundwater samples (Fig. 9) are generally higher in the northern (north of ∼31°57′N) part of the study area. This is likely due to groundwater interaction with the Paleozoic carbonates near the basin margins (Fig. 1), as well as infiltration of surface water flowing over the Paleozoic carbonates of the northern Franklin Mountains. The few high Ca/Na ratios in the southern part of the study area may be related to interaction of water with Cretaceous carbonate bedrock that is <100 m from surface (see profile NW-SE in Fig. 6).

Silicon (Si) concentrations across the study area range from 6.75 mg/L to 31.75 mg/L (Fig. 10). The primary unit penetrated by wells within the study area are the sands and silty clays of the middle Santa Fe unit. Silicon concentrations tend to increase with the age of the groundwater due to intensive water-rock interaction along the flow path in the aquifer. The relatively low or normal Si concentrations in most of these samples suggest that the majority of water in the basin is generally young, with recent recharges that have not remained within the basin long enough to dissolve a significant amount of Si from silicate minerals. However, the highest concentrations occur in the west central part of the study area (Fig. 10) within 2 km of the Mesilla Valley fault zone. Well water that has high Si concentration may indicate deeper, older water upwelling along the Mesilla Valley fault zone.

Ancient seawater, which may be contained in deep sedimentary basins, exhibits low Cl/Br ratios resulting from evaporation past the point of precipitation of halite, which preferentially excludes Br from its lattice (Davis et al., 1998). Typically, deep residual brine waters also exhibit low sulfate (SO4) concentrations due to sulfate-reducing reactions that remove dissolved SO4. Hence, deep basin brines are generally associated with low Cl/Br ratios and low SO4 concentrations. However, such geochemical indicators can be modified by dissolution of evaporite minerals. For example, sulfate levels are also affected by the dissolution of gypsum (CaSO4•H2O). During gypsum dissolution, both sulfate and calcium are released into the groundwater and their concentrations increase. The dissolution of halite increases Cl/Br ratios. If upwelling brines from the bedrock are passing through gypsum- and halite-dominated, evaporite-rich formations such as the middle Santa Fe unit, higher sulfate and calcium concentrations, in combination with high Cl/Br ratios, are expected. Moreover, irrigation return waters from agricultural areas have also been found to contain high levels of sulfate and high Cl/Br ratios (Szynkiewicz et al., 2011). These waters seep into the underlying shallow aquifer units where they mix with the groundwater and elevate sulfate levels and Cl/Br ratios.

Most brackish wells (TDS >1000 mg/L) within the study area show high Cl/Br ratios (>200) (Fig. 11). This indicates that saline groundwater in this area is not a product of deep basinal brine upwelling. Instead, the saline source is evaporite minerals, mainly gypsum and/or mirabilite-thenardite, as well as halite, and is being dissolved in the groundwaters in the middle to lower Santa Fe units. There are a few wells within the study area, however, that exhibit a combination of lower SO4 (≤470 mg/L) concentrations and lower (<200) Cl/Br ratios (Fig. 11) that could indicate upwelling of deeper basin brines. A majority of these wells are located within the Mesilla Valley fault zone (Fig. 11); a few lie in the southernmost part of the study area where bedrock is <90 m below the surface (Fig. 11). It should be noted that higher concentrations of Br can also be associated with human-impacted surface-water runoff that contains higher amounts of Br from soaps, cleaners, and swimming pool chemicals. This might be another source of anthropogenic Br for the shallow wells in the southernmost portion of the study area, which is the ultimate groundwater discharge zone for the entire valley (Witcher et al., 2004) before it is discharged through the Paso del Norte into the Hueco Basin. Additional studies, including the use of anthropogenic isotope tracers (such as boron isotopes; Garcia, 2017), could reveal more information to distinguish the natural versus anthropogenic sources of Br in groundwater in the southernmost part of the study area.

Groundwater chemistry in this study shows a systematic spatial variation with respect to the distance of wells from known and inferred faults in the Mesilla Basin (faults shown in Fig. 5). For example, groundwater well samples with high TDS (with Cl >500 ppm) are generally located within ∼500 m of known faults (Fig. 12). In addition, samples from these wells also show unusually high concentrations of trace elements such as Fe (up to ∼180 ppb) and As (up to 100 ppb). The associations of high TDS, Fe, and As concentrations with known fault locations suggest that the fault zones in this area act as conduits for groundwater flow that are most likely connected to the deeper part of the lower Santa Fe unit. The dissolution of halite within the lower Santa Fe unit is responsible for the high Cl concentrations, and the deeper and reducing conditions likely lead to dissolution of sulfide minerals that contributes to the high dissolved Fe and As concentrations in this group of samples.

Witcher et al. (2004) classified geothermal water within the Mesilla Basin as water having a temperature >26 °C. Our field results (Table S1 [footnote 1]) show that 102 wells within study area have geothermal water (>26 °C), and only 20% of these wells are >300 m deep. Only five of these wells are not located in the region between the Western and I-10 faults. In addition to the higher temperatures, 53% of these wells have TDS values in >500 mg/L, which is within the typical range (500–5000 mg/L; Witcher et al., 2004) for Mesilla Valley geothermal water. Higher TDS concentrations in geothermal waters are attributed to a combination of lateral or vertical movement of water within a basin and the resulting dissolution, reprecipitation, and evapotranspiration occurring closer to the surface (Witcher et al., 2004). Although not all of these geothermal wells show the typical chemistry of deeply derived basinal brine water, they show that warmer water is upwelling from deeper levels within the lower Santa Fe unit, with the Mesilla Valley fault zone and the River fault acting as likely conduits for groundwater flow.

Combined studies of gravity and groundwater chemistry data suggest that the Mesilla Valley fault zone and the River fault act as conduits for the upwelling of deeper water within the Mesilla Basin. Wells near these faults tend to have higher temperatures, elevated Si, Cl, Fe, and As concentrations, and lower SO4 and Cl/Br ratios than those further from the faults, resulting from water-rock interactions in the lower Santa Fe unit. Consistent with this interpretation, several studies of the Rio Grande near the Mesilla Valley fault zone and the River fault have suggested that the river chemistry is directly affected by the upwelling of deeper groundwater in this region (e.g., Phillips et al., 2003; Hogan et al., 2007). Recent isotope studies of the Rio Grande (Szynkiewicz et al., 2015a, 2015b; Nyachoti, 2016; Garcia, 2017) have shown consistently elevated (234U/238U) ratios (where the parentheses indicate an activity ratio) in the river during low-flow seasons (October to April), indicative of upwelling of deep groundwater that is characterized by high (234U/238U) ratios due to intensive alpha recoil effects in aquifers. Indeed, the Mesilla Valley fault zone and the River fault intercept the course of the Rio Grande at several locations in this area (Fig. 1), and the impacts of groundwater upwelling on Rio Grande salinity at several key locations in the Mesilla Basin have been studied by a number of previous researchers (e.g., Phillips et al., 2003).

The higher Ca/Na ratios (Fig. 9) observed in the northeastern part of the study area suggest the influence of groundwater interaction with carbonates of the northern Franklin Mountains. The I-10 fault may be serving as a recharge zone for surface runoff of carbonate-rich water into the aquifer system. Some of the lowest Si values are also observed near the I-10 fault (Fig. 10), indicating more modern recharge in this region.

Groundwater chemistry at the southern edge of the basin indicates significant water interaction with the shallow carbonate and igneous bedrock in this region, and possible interaction with human-impacted surface water enriched in Br. This is not unexpected due to the shallowing of the basin to <90 m depth and resulting rock-water interactions with the Eocene igneous intrusions and Cretaceous carbonate bedrock. A more detailed study of the local water chemistry and bedrock geology is required to unravel some of this observed complexity.

Analysis of gravity data shows a complexity of shallow subsurface structures with abrupt variations in lithology type and thickness (Fig. 6) compared to the deeper Hueco Basin (e.g., Avila et al., 2016). Gravity analysis has aided in determining the subsurface extent of Eocene intrusions (see Fig. 4), but more extensive three-dimensional (3-D) modeling at a regional scale is required to better image features that may be related to the edge of the Chihuahua trough, to the transition between thin-skinned and thick-skinned Laramide deformation and to more recent rift extension. We are in the process of collecting more gravity data to adequately image the structure of the Mesilla Valley fault zone (Fig. 2), one of the major controls on basin structure and hydrology.

Our density models (Fig. 6) are in good agreement with most features of the hydrogeologic cross sections of Hawley and Kennedy (2004) that extend to depths of 500–1000 m from the surface. On some cross sections, the gravity data suggest thicker sediment packages, more extensive Eocene intrusions, or additional faults relative to Hawley and Kennedy’s interpretations, but none of our interpretations violate observations obtained directly from water-well information.

A combined study of gravity and groundwater geochemistry information suggests that the Mesilla Valley fault zone (as mapped by Hawley and Kennedy [2004]) and the River fault (inferred by Khatun et al. [2007] and Sellepack [2003]) serve as conduits for deeper (e.g., from the lower Santa Fe unit), warmer, more Si- and Cl-rich waters with lower Ca/Na and Cl/Br ratios relative to shallow groundwater and low SO4 in the northern part of our study area. The I-10 fault (inferred by Khatun et al. [2007]) may serve as a recharge zone for high-Ca/Na surface runoff from the dissolution of carbonates of the northern Franklin Mountains, while the Three Sisters fault (inferred by Khatun et al. [2007] and Sellepack [2003]) may be an older feature that served as a zone of weakness for the localization of Eocene intrusions. The high Cl/Br ratios found in ∼90% of wells in the study area (Fig. 11), however, indicate that water salinity is not a product of deep residual basin brines, but is a result of dissolution of halite, gypsum-selenite, and/or mirabilite-thenardite found in the middle and lower Santa Fe units. Groundwater chemistry in the southern Mesilla Valley indicates a complex interaction between shallow subsurface Eocene igneous intrusions and Cretaceous carbonates, coupled with upflow of deeper waters that likely flow along east-west–trending accommodation structures that eventually discharge the groundwater through the Paso del Norte into the Hueco Basin.

Gravity modeling suggests that the greatest offset of the Santa Fe units has occurred along the Mesilla Valley fault zone and Western fault. Although other faults (e.g., the River, I-10, and Three Sisters faults) in the southern basin offset bedrock and appear to play an important role in water flow, they likely played a greater role in Laramide or early Cenozoic deformation. Lack of gravity data for the region surrounding the Mesilla Valley fault zone made it difficult to image its structure. We are currently collecting gravity readings in this area to remedy this problem. Analysis of gravity data on a regional scale with the use of 3-D modeling techniques should assist in revealing subsurface structures associated with Jurassic and Cenozoic rifting, Laramide compression, and Eocene to recent volcanism.

This study highlights the unique resolving power of combining geophysical and geochemical techniques in understanding groundwater chemistry and its relationships with local geological structures. Such an approach can be readily applied to other systems with similar geologic and hydrological settings for groundwater exploration and resource management.

We thank M. Baker, M. Engle, G. Kaip, and C. Montana for field and technical support. Students who assisted in the field work and data processing and analysis included: Y. Abushalah, R. Alfaro, L. Baker, R. Guzman, J. Lucero, M. Lucero, M. Moncada, S. Nyachoti, J. Ochoa, A. Ramirez, C. Reyes, R. Rivera, M. Sandoval, and F. Ziwu. We appreciate comments by the associate editor, reviewer A. Al-Taani, and an anonymous reviewer that helped us in revising this manuscript. This research was supported in part by National Science Foundation grant EAR-1349091 to L. Ma. M. Hiebing received scholarship support from the Hunt and Rowling scholarship funds, the West Texas Geological Society, and the Southwest Section of the American Association of Petroleum Geologists. V. Avila was supported by a scholarship from the Consejo Nacional y Tecnología (CONACYT).

1Supplemental Files. Table S1: Geochemical data for water-well and river-water sampling locations shown in Figure 3A. Table S2: Groundwater geochemistry data from other sources. Please visit or access the full-text article on to view the Supplemental Files.
Science Editor: Raymond M. Russo
Associate Editor: Roberto Stanley Molina-Garza
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.