Geological and geophysical studies of the structure and stratigraphy of the northwestern Hueco Bolson Aquifer, El Paso, Texas

The Hueco Bolson Aquifer contributes ∼30% to 40% of the yearly drinking water supply to El Paso’s population when sufficient water can be taken from the Rio Grande (El Paso Water Utilities, http://www.epwu.org/water/water_resources.html [accessed 30 June 2016]). In drought years, when river water is minimal, the Hueco Bolson Aquifer supplies an even greater percentage of the local water supply (El Paso Water Utilities, http://www.epwu.org/water/water_resources.html [accessed 30 June 2016]). The Hueco Bolson Aquifer contains both freshwater and brackish water, and El Paso has built one of the largest desalination plants in North America (El Paso Water Utilities, www.epwu.org/water/desal_info.html [accessed 20 December 2017]) in order to efficiently process brackish water from the Hueco Bolson Aquifer. Thus, it is important to understand the structural and stratigraphic controls on the quantity and quality of water within the Hueco Bolson Aquifer, where the groundwater level has declined by 60 m in some areas under the Hueco Bolson Aquifer over a period of <90 yr (Heywood and Yager, 2003).

The Franklin Mountains and the Hueco Mountains define the eastern and western boundaries of the northern Hueco Bolson (basin), which extends southward into Mexico (Fig. 1). The northern boundary is a groundwater divide separating the Hueco Bolson from the Tularosa Basin, located just north of the Texas–New Mexico State line. However, our study is limited to the El Paso region as outlined in Figure 1.

Work by Marrufo (2011) suggested that faults serve as the barrier between freshwater and saline water in the central Hueco Bolson Aquifer; however, we do not know if the faults are the controlling mechanism for the freshwater–brackish water contact in the northern part of the Hueco Bolson Aquifer. Other faults also likely serve as conduits for water flow within the Hueco Bolson Aquifer. In addition, basin stratigraphy plays a role in the location of fresher water lenses within the northwestern Hueco Bolson Aquifer. Increasing our understanding of the structure and stratigraphy of the Hueco Bolson Aquifer will allow us to better estimate the volume of freshwater and brackish water for sustainable resource management. Finally, many faults within the Hueco Bolson show evidence for Quaternary displacement (e.g., Collins and Raney, 2000). Better knowledge of their locations and offsets can be used to update seismic hazard analyses of the region. The most recent felt earthquake in El Paso, on 6 March 2012 (U.S. Geological Survey, 2012), is believed to have occurred along one of these faults.

The Hueco Bolson is one of a series of north-trending, interconnected asymmetrical grabens that form the Rio Grande rift system (e.g., Chapin, 1971; Seager et al., 1984; Collins and Raney, 1994; Keller and Baldridge, 1999). The Hueco Bolson consists of two subbasins, the northern subbasin, the focus of this study, and a southern subbasin located ∼60 km southeast of our study area. The northern subbasin contains predominately north-striking faults (Fig. 1), while the southern subbasin contains northwest-striking faults. The Rio Grande rift has experienced two phases of extension, producing two different styles of basins (Seager et al., 1984; Morgan and Golombek, 1984). Early extension (mid-Oligocene to early Miocene) resulted in formation of broad and relatively shallow basins with low-angle normal faults, while the later extension (mid-Miocene to Quaternary) produced high-angle, fault-bounded, relatively deep and narrow basins (Morgan and Golombek, 1984; Keller and Baldridge, 1999). Extension of the Rio Grande rift is still occurring, with evidence for a least four earthquakes along the East Franklin Mountains fault, the fault bounding the western edge of the Hueco Bolson Aquifer, within the past 64 k.y. (McCalpin, 2006).

Basin fill in the Hueco Bolson is primarily composed of unconsolidated to poorly consolidated fluvial, lacustrine, and eolian sediments of Miocene to Holocene age (Fig. 2; Hadi, 1991; Anderholm and Heywood, 2003; Heywood and Yager, 2003; Buck et al., 1998), but little is known about older Cenozoic sediments due to lack of outcrop and core data (Collins and Raney, 1991; Hawley and Kennedy, 2004). The basin fill deposits can be broadly divided into lower basin fill, upper basin fill, and surficial deposits (Fig. 2; Collins and Raney, 1991). The Pliocene–Pleistocene Fort Hancock and Camp Rice Formations (Strain, 1966) represent the upper part of the upper basin fill. The southern Hueco Bolson subbasin contains the thickest lower basin fill, while the northern Hueco Bolson subbasin contains a thicker sequence of upper basin fill (Collins and Raney, 1994). This difference in subbasin depositional history suggests that the northwest-striking faults were most active in the early formation of the subbasins, while the north-striking faults were more active during deposition of the upper basin fill (Collins and Raney, 1994).

The Fort Hancock Formation mostly consists of lacustrine and proximal to distal alluvial-fan deposits (Gustavson, 1990). Fluvial facies, likely deposited by a meandering stream system, have also been reported (Willingham, 1980; Riley, 1984). A regional unconformity separates the overlying Camp Rice Formation from the underlying Fort Hancock Formation (Strain, 1966). The Camp Rice Formation mostly consists of sand and gravel of braided stream systems and alluvial fans. It also contains lacustrine and flood deposits in minor amounts (Collins and Raney, 1991). Series of Pleistocene gravel units, the Miser, Madden, Gills, Ramey, and Balluco gravels (Fig. 2), were deposited on piedmont slopes and arroyo terrace environments that overlie the Fort Hancock and Camp Rice Formations (Albritton and Smith, 1965). These are covered by windblown fluvial and alluvial deposits of the present Rio Grande system.

Heywood and Yager (2003) divided the sediments of the Hueco Bolson into fluvial, alluvial-fan, lacustrine-playa, and recent alluvial facies. Based on grain-size analysis and well logs, Marrufo (2011) identified an additional playa-margin facies. The fluvial facies consists of relatively coarse-grained sand and gravel interbedded with silt and clay that are primarily related to the ancient Rio Grande River system (e.g., Heywood and Yager, 2003; Hutchison, 2004; Marrufo, 2011). The predominant geologic formation of the fluvial facies is the Camp Rice Formation. The alluvial-fan facies of the Camp Rice Formation consists of poorly sorted gravel and sand mainly derived from the present-day Franklin Mountains. The Fort Hancock Formation contains lacustrine-playa and playa-margin facies related to Lake Cabeza de Vaca, the terminus of the ancestral northern Rio Grande until its integration with the lower Rio Grande ca. 2.25 Ma (Gustavson, 1990). These lacustrine-playa and playa-margin facies occur in the central, eastern, and southeastern parts of the basin beneath the fluvial and alluvial-fan facies.

Marrufo (2011) demonstrated the success of using a combination of water well logs, well cuttings, and microgravity data to locate faults that appear to separate zones of freshwater and saline water in the Hueco Bolson Aquifer within central El Paso. The logs and cuttings from Marrufo’s study also helped to define a sedimentological and stratigraphic framework that could be applied to other parts of the Hueco Bolson Aquifer. In our current study, we collected additional information from water wells and microgravity data from the northwestern Hueco Bolson Aquifer (Fig. 3) in an effort to trace the faults and stratigraphic facies across the basin.

This study is an extension of previous gravity surveys carried out in El Paso and surrounding areas (e.g., Figuers, 1987; Burrows, 1984; Hadi, 1991; Burgos, 1993; Marrufo, 2011; Avila et al., 2016). In our study, we collected 502 new gravity data points (Fig. 3, green symbols) and combined them with existing gravity data available from the University of Texas at El Paso (UTEP) gravity database (https://research.utep.edu/default.aspx?tabid=37229 [accessed 20 December 2017]); Avila et al., 2016; Avila, 2016).

We collected gravity data with a 75–300 m spacing using a Lacoste and Romberg model G-1115 gravimeter and Topcon GB 1000 global positioning system (GPS) with ground plane. This gravimeter has an accuracy of 0.01 mGal. The GPS base points were submitted in receiver independent exchange (RINEX) format to the National Geodetic Survey via the Web-based Online Positioning User Service (OPUS). OPUS processes the GPS points with respect to three reference sites, and the processed points are sent back to the user via e-mail (Mader et al., 2003). We postprocessed our data using OPUS solutions, which allowed us to achieve a high-accuracy estimate of elevation (<5 cm) and associated latitude and longitude.

A secondary base station was established near Painted Dunes Golf Course (Fig. 3, blue diamond) by tying it to the absolute base station situated at the Kidd Memorial Seismological Observatory on the UTEP campus. Gravity points were collected in loop style, and survey loops were closed at the golf course base station every 2–3 h to reduce the effects of tides and instrument drift (Lowrie, 1997). The repeatability of the gravity value at base stations was ∼0.01–0.04 mGal. The golf course base station was tied at the beginning and end of each day to the absolute base station at UTEP. We converted the gravity readings from dial units to gravity units and then applied corrections for tides/drift, latitude, elevation (free-air), terrain, and density (Bouguer).

Figure 4 shows a complete Bouguer anomaly map that was generated using Geosoft’s Oasis Montaj software with a 200 m grid size and minimum curvature interpolation technique. An average crustal density of 2670 kg/m³ was chosen for the Bouguer correction. The terrain correction used digital elevation models from the U.S. Geological Survey.

The complete Bouguer anomaly map shows a strong correlation with the structural highs and lows known from previous studies (e.g., Burrows, 1984; Hadi, 1991; Burgos, 1993; Avila et al., 2016). High gravity anomalies are related to the high-density rocks of the Franklin Mountains, and low anomalies are associated with low-density basin fill deposits. Gravity values decrease from west to east, indicating increasing basin fill to the east. The eastern portion of the study area shows the lowest gravity anomaly (∼–169 mGal). The East Franklin Mountains fault on the east side of the Franklin Mountains and other smaller faults are reflected in the gravity anomaly map (Fig. 4). The gravity signatures of some of these small-scale structures are not necessarily identifiable in the complete Bouguer anomaly map alone.

To interpret these smaller anomalies, we created several short profiles comparing simple Bouguer anomaly values to topography in the flat-lying portions of our study area as shown in Figure 5. Profiles were created so that all gravity values were located within 50 m of the profile line. Figure 5, with data taken along McCombs Street (Figs. 3 and 4), generally shows that the simple Bouguer anomaly values change inversely with changes in topography. For example, as topography increases between 0 and 500 m and 2000 and 3000 m along the profile, the simple Bouguer anomaly values decrease. However, between distances of 1300 and 1800 m, the anomaly values increase by 0.3 mGal, although the elevation changes very little (<0.5 m). One possible explanation for the observed anomaly increase could be if we assumed that points T19 and T20 were on the footwall block of a concealed normal fault that lies somewhere between points T20 and T21. We used this assumption (i.e., changes in simple Bouguer anomaly values that are not correlated with changes in topography could be due to concealed normal faults) to locate possible positions of normal faults along our profiles. At least four of these possible faults (faults 8, 9, 10, and 22) align very well with the strike and sense of motion along faults mapped at the surface, and we extended these faults as indicated by the purple lines in Figure 4.

A residual Bouguer anomaly map (Fig. 6) was obtained by fitting a third-order polynomial surface to the complete Bouguer anomaly map for the entire northern Hueco Bolson and subtracting the polynomial to eliminate deeper-seated regional features in the area (Avila, 2016). The residual map shows the lowest values (<–10 mGal) in the south-central portion of the study area extending to ∼31.90°N. There is a broad region of moderately low residuals (–1 to –4 mGal) throughout much of the central study area that decrease at the extreme northeastern edge of the study area, indicating the southern edge of the Tularosa Basin.

The horizontal gradient magnitude method (HGM) was applied to the residual Bouguer anomaly values to help further delineate faults. We used the method of Grauch and Johnston (2002) in this analysis. HGM values should be highest (i.e., show the greatest change in gradient) over the edge of a fault if the fault is near vertical or is not adjacent to another near-vertical fault or structure. The HGM map is shown in Figure 7.

The highest HGM values are associated with the East Franklin Mountains fault (Fig. 7), but other mapped surface faults do not appear to be associated with high HGM values, perhaps because they are closer together and have less vertical offset. The only HGM highs that appear to suggest additional faults are in the region of the East Franklin Mountains fault step-over (fault 29) and in the region just west of faults 16 and 17, which we have labeled fault 16N.

We selected two profiles in the study area along which to model gravity data (Fig. 4) using Geosoft’s GM-SYS software, which is based on the forward modeling technique of Talwani et al. (1959). The profiles were constructed to pass through the regions of densest gravity observations. Regional geology, subsurface lithostratigraphy, density, and thickness of formations are reasonably well understood from previous studies (e.g., Burrows, 1984; Collins and Raney, 1991; Hadi, 1991; Burgos, 1993; Hawley and Kennedy, 2004; Granillo, 2004; Marrufo, 2011; Avila et al., 2016; see Table 1 for density information). Geologic maps and cross sections by Collins and Raney (1994, 2000), Hawley et al. (2009), Marrufo (2011), and Avila et al. (2016) were used for geologic contacts and fault constraints. The model results are shown in Figure 8.

Profile A-A′ is an east-west–trending profile located ∼2 km south of the Texas–New Mexico State line at the northern edge of the study area (Fig. 4), and it overlaps with the hydrogeologic cross section B-B′ of Hawley et al. (2009). The profile passes across several faults, including the East Franklin Mountains fault (Fig. 8). Faults I, II, III, IV, VII, XI, and XV are faults previously mapped by Collins and Raney (2000; see also Fig. 4). Faults I and III were also noted by Hawley et al. (2009) on their cross sections. We interpreted three other gravity anomalies of ∼0.5 mGal as possible northeast-dipping normal faults (faults 3, 4, and 18; Fig. 4). Many of the faults cut across the profile at an oblique angle. Note that the deepest part of the basin is more than 13 km east of the East Franklins Mountain fault, reaching a depth of ∼1450 m near fault 4 (Fig. 8). Based on oil and water well information, Hawley et al. (2009) suggested that the valley fill is 1300 m thick near fault I and over 1500 m thick near fault 4, values that are consistent with our density model. Avila et al. (2016) modeled gravity data along a profile located ∼6 km north of profile A-A′ and suggested that there was 900 m of basin fill adjacent to the East Franklin Mountains fault and ∼1800 m of basin fill near fault 3.

Faults that appear to have offset the contact between the upper and lower basin fill layers in the density model include faults I, II, III, IV, 4, and XI and an unnamed fault lying between faults 4 and III. Fault 18 appears to offset the lower basin fill, but not the upper basin fill. Fault XV may have a minor offset of the upper basin fill, as suggested by the shape of the upper-lower basin fill contact (Fig. 8). Faults VII and 3 do not appear to offset this contact.

Profile B-B′ (Fig. 8) is also a west-east–trending profile (Fig. 4). This profile is located just north of a point where the East Franklin Mountains fault steps over ∼2 km to the west and where intrabasin faults begin to change strike from north-south to northwest-southeast. The profile lies within 2–3 km of the hydrogeological profile C-C′ of Hawley et al. (2009). Avila et al. (2016) modeled gravity data along a northwest-southeast–trending profile located ∼5 km south of profile B-B′. Faults VII, VIII, IX, X, and XVI are faults previously mapped by Collins and Raney (2000). The profile cuts across faults 9, 10, 14, and 15, which are based on our interpretation of gravity anomalies (Fig. 4). As previously noted, we infer that faults 9 and 10 are extensions of faults IX and X, respectively, based on their proximity and similarity of strike and dip (Fig. 4). The cross section of Hawley et al. (2009) also shows faults that correspond to the positions of faults 9, VII, and VIII. Very little basin fill is modeled adjacent to the East Franklin Mountains fault (Fig. 8). The basin fill thickens significantly ∼4 km east of the East Franklin Mountains fault, reaching a maximum of 1400 m fill located ∼10 km east of the fault. Avila et al. (2016) estimated a maximum of 1800 m of basin fill ∼8 km east of the East Franklin Mountains fault along their profile located to the south of B-B′.

Our model suggests the presence of a fault (fault 29) located ∼800 m east of the East Franklins Mountains fault that offsets lower basin fill deposits. The presence of this fault is also suggested by the HGM analysis (Fig. 7). The next fault to the east offsets both the lower and upper basin fill and may be an extension of faults 16 and 17 (fault 16N), which is supported by the HGM analysis (Fig. 7). Faults 14 and 15 appear to be very minor features, if they exist at all. The density profile suggests the presence of a fault between faults 14 and XVI (Fig. 8) that is not apparent on the HGM map; we label it fault 30. Faults 30, 16N, XVI, 10/X, 9/IX, VIII, and VII in our density model offset both upper and lower valley fill (Fig. 8), while fault 29 may only offset the lower valley fill.

Digital well logs obtained from 11 El Paso Water Utilities (EPWU) water wells (THNH1, THNH2, TH43A, E-3A, MNST3, MNST5, MNST6, FBT04, 29B, 40A, and 601) were used in this study (Figs. 3 and 4). Marrufo (2011) and Doser and Langford (2008) analyzed grain-size variation from cuttings in wells 601, 605, 610, 615, and 509A (Fig. 3), located to the south of our study area, and interpreted depositional environments based on grain-size distribution and log responses. The EPWU provided simple cuttings descriptions based on the unified soil classification system for most logs we analyzed in this study. These descriptions do not help much in inferring depositional environments. The proximity of our study area to that of Marrufo (2011) allowed us to use cuttings from well 601 as a reference, and we compared these to log responses (Fig. 9).

We used the “electrosequence analysis” concept (Rider, 2002) to determine geological information from the geophysical logs. We used gamma logs (Fig. 9, left logs) because only these records were available for all wells; however, multiple logs are desirable for this technique. We identified depositional facies through log responses and then grouped them together to form facies associations. After establishing characteristic facies, we identified surfaces such as flooding surfaces, erosional surfaces, or unconformities and used these for well-log correlation. For example, we interpreted sandstone resting on the top of shale and shale sitting on the top of sandstone as erosional and flooding surfaces, respectively (see Fig. 9). Figure 9 also shows abrupt breaks in the log that indicate changes in lithology or a structural break (unconformity, fault, etc.).

We inferred five facies (channel fill, proximal alluvial fan, distal alluvial fan, floodplain-deltaic, and sandy-silty) from logs of the Hueco Bolson, and then we grouped them into four facies associations that corresponded to different depositional environments (Table 2). Log characteristics, stacking patterns, and surfaces used to characterize depositional facies are shown in Figure 9 and Table 2. Cuttings from channel fill facies show rounded to well-rounded gravels that grade upward into sands. On gamma-ray logs, these are distinguished by bell-shaped patterns with a sharp base, low to moderate gamma-ray readings, and high resistivities (Fig. 9, ∼180 ft [55 m]). Proximal alluvial-fan facies consist of a variety of poorly graded gravels, silty gravels, and clayey gravels, sometimes interbedded with silt and clay. These facies show serrated to slightly coarsening-upward gamma-log signatures (e.g., ∼250 ft [76 m]; Fig. 9). Distal alluvial-fan facies consist of well-graded gravel, clayey gravel, clayey sands, and clay. Gamma-log responses for the facies show serrated, bell, funnel, and dumbbell patterns (e.g., ∼220 ft [67 m]; Fig. 9). The funnel-shaped pattern is especially distinctive of this facies. Floodplain-deltaic facies are characterized by very thick, multiple funnel-shaped gamma-log patterns stacked on top of each other and separated by high gamma peaks (e.g., 450–530 ft [137–162 m]; Fig. 9). The sandy-silty facies shows either a serrated or blocky pattern (with internally serrated structure) in gamma logs (e.g., 770–810 ft [235–247 m]; Fig. 9). More details on the facies and their associated patterns may be found in Budhathoki (2013).

Due to the limitation of available logs, the distances between wells, and the depositional complexity of the area, it was not easy to correlate all facies across the region. In order to overcome this problem, we grouped the preceding five facies into four facies associations according to environment of deposition. Figure 10 is a schematic representation of the depositional environments of the study area.

Facies association 1 is interpreted to have been deposited in a fluvial environment. The channel fill facies is the dominant facies of this association. The association consists of thin- to thick-bedded, single- to multistory or isolated channel fill deposits with episodes of flooding. Due to its proximity to the Franklin Mountains, larger grain size, stacking patterns, and the lack of observed point bar deposits, the channels are interpreted to be predominantly associated with braided systems.

Facies association 2 is mainly dominated by gravel deposits in the proximal part of alluvial fans and gravels to coarse-grained sand deposits in the more distal part of the fans. The gravels in the distal alluvial fans are variably rounded, with limestone and granite matrix, and contain finer-grained sands than fluvial deposits (Doser and Langford, 2008). Doser and Langford (2008) also observed coarse-grained sands with sparse gravel clasts and felt this facies was comparable to exhumed intervals of the Santa Fe Group found in the Hueco and Mesilla Bolsons.

Floodplain and deltaic facies are grouped together in facies association 3 and interpreted to represent playa-margin and deltaic deposits. This sequence is mainly composed of clayey silts and silty clay with sands and gravels. Doser and Langford (2008) and Marrufo (2011) reported carbonate nodules and white silty sands and sandy silts in this unit.

Facies association 4 contains sandy-silty facies and shows repeated coarsening-upward cycles in the gamma-log response records. From the cuttings analysis, 3–30-m-thick, variegated silty clays and clayey silts have been reported for this association. When sediment enters into a water body, the suspended load settles out, leading to accumulation of fine-grained sediments. The thick clay indicates a stable period of time, whereas alternating sands and clays or silty-clayey successions are the result of multiple sheetflood events or repetitive stream invasion (Mader, 1985).

Four structural cross sections (for locations, see Fig. 4) are shown in Figures 11 through 14, with their tops at ground elevation. Well logs were correlated using Schlumberger’s Petrel software suite. Seven key erosional surfaces (ES1–ES7) were identified based on gamma and resistivity signatures that could be correlated between wells.

Figure 11 shows a cross section of wells extending from the foothills of the Franklin Mountains to east of Railroad Drive and incorporates wells THNH1, TH43A, 29B, 40A, and FBT04. Structure along this cross section should be similar to that observed on gravity/density profile B-B′ (Fig. 8). Faults 14, 15, and 16N (all east dipping) are consistent with the down-to-the-east offsets of all erosional surfaces of 60–260 ft (18.3–79.3 m) between wells THNH1 and TH43A, but profile B-B′ (Fig. 8) suggests that fault 16N is likely responsible for the observed offsets between wells. There is also a down-to-the-east offset in erosional surfaces 1–6 of 40–145 ft (12.2–44.2 m) between wells TH43A and 29B, consistent with profile B-B′ (Fig. 8), where we inferred fault 30 to be located. East-dipping fault XVI (Figs. 4 and 11) appears to offset erosional surfaces 1–5 between wells 29B and 40A from 20 to 55 ft (6.1–16.8 m). Profile B-B′ (Fig 8) indicates that the major faults located between wells 40A and FBT04 are faults 9/IX (west dipping) and 10/X (east dipping). Figure 11 indicates that three erosional surfaces appear downdropped to the east, three are downdropped to the west, and one appears to have little offset between wells 40A and FBT04. This suggests there has been a complicated history of movement along these faults.

The cross section shown in Figure 12 extends from north to south and passes through wells E-3A, THNH2, and TH43A (Fig. 4). Fault 18, observed on profile A-A′ (Fig. 8), likely crosses the profile between wells E-3A and THNH2. The downward offset of all erosional surfaces of 120–335 ft (36.6–102.1 m) to the south between these wells is consistent with movement on fault 18, which dips to the southeast. Downward offsets of erosional surfaces to the north between wells THNH2 and TH43A could be related to fault 30 (Fig. 8, profile B-B′), the strike of which is poorly determined.

Figure 13 shows a northwest-southeast–oriented cross section that includes wells E-3A, MNST3, MNST5, and FBT04 (Fig. 4). All erosional surfaces between wells E-3A and MNST3 suggest down-to-the-southeast offset of 155–245 ft (47.2–74.7 m), consistent with movement along fault XI as indicated on profile A-A′ (Fig. 8). Between wells MNST3 and MNST5, the well-log data suggest small offsets of 0–65 ft (0–20 m) of the erosional surfaces, consistent with Figure 4, which indicates the cross section is located parallel to major faults in the region. The five deepest erosional surfaces between wells MNST5 and FBT04 have been downdropped 30–65 ft (9–20 m) to the northwest, but there is little offset (<20 ft [6 m]) of erosional surfaces 1 and 2. The northern ends of faults X and IX (Figs. 4 and 8) should cut between these wells, and similar to what was observed on the cross section of Figure 11, there appears to be a complicated history of movement along these two faults.

The final cross section (Fig. 14) trends north-south and extends from well MNST6 through wells MNST3 and 40A, ending at well 601 in the south (Fig. 4). Figure 4 indicates that this cross section should parallel major faults within the region, except at well 40A, where it may cross the southern end of fault XVI. The cross section indicates ∼25–90 ft (7.6–27.4 m) of consistent southward offset of all erosional surfaces between wells MNST6 and MNST3. Between wells MNST3 and 40A, erosional surfaces 5–7 are displaced downward to the north, and erosional surfaces 1–3 are displaced downward to the south, and between wells 40A and 601, there is 10–110 ft (3–33.5 m) of southward offset of the surfaces. The offset between wells 40A and 601 may be related to fault XVI, which cuts obliquely across the profile.

Based on these structural cross sections, faults 16N, 18, XI, and XVI, and possibly fault 30, appear to have consistently displaced all seven erosional surfaces in the same direction and have offset at least three of the seven surfaces by ≥30 m. These faults are thus expected to have the greatest influence on groundwater movement within the basin. In addition, faults X/10 and IX/9 appear to have a complex history of differential movement that may serve as a barrier to groundwater flow.

The salinity of the Hueco Bolson Aquifer varies from well to well, and the groundwater quality generally degrades with depth (Hutchison, 2004). Our study area has few wells, and they are located far apart. This makes it difficult to judge the true distribution of the total dissolved solids (TDS) across the study area. It has also been observed that the chloride limit exceeds drinking water standards before TDS limits are reached in some of the wells (Hutchison, 2004). For this reason, chloride concentration data can be used alternatively to map the distributions in salinity (Hutchison, 2004). A chloride concentration between 250 and 500 mg/L is considered as transitional between freshwater and brackish water.

We compared the chloride concentrations from Hutchison (2004) to our fault interpretations, as shown in Figure 15, to determine how the presence of faults may affect the occurrence of freshwater and saline water in the study area. Hutchinson’s interpretations were based on water chemistry information from ∼50 wells in our study area and used a three-dimensional kriging technique with a grid size of ∼500 m to develop a series of concentration versus depth maps. Figure 15A shows the 250, 500, and 1000 mg/L chloride concentration contours at ∼1084 m elevation (∼135 m below the surface), and Figure 15B shows the contours at ∼731 m elevation (∼488 m below the surface). Faults that have been imaged by gravity (HGM analysis and/or density modeling from this study; Avila et al., 2016; Marrufo, 2011) or surface mapping (Collins and Raney, 2000) and that are consistent with well logs (from this study; Hawley et al., 2009; Marrufo, 2011) are indicated by bold lines. Note that fault 30 is consistent with the density profile along profile B-B′’ (Figure 8) and well logs, but it was not recognized in the HGM or simple Bouguer analysis. Its exact location and strike are thus not well known. In Figure 15, we suggest that fault XVI is linked to fault XI through fault 22, and that fault XVI may also be linked to fault M1 by fault 12 (dashed line).

At shallow depths (Figure 15A), brackish to saline water (chloride concentrations 250 to over 1000 mg/L) is found in the region between the East Franklin Mountains fault and the approximate location of fault 30. Most freshwater (chloride concentrations <250 mg/L) south of 31.93°N is located between the East Franklin Mountains fault and faults X and IX. North of 31.93°N, freshwater is found between faults 30 and VII. Small pockets of more brackish water (chloride concentrations from 250 to 750 mg/L) are located between faults XVI and X, in a region southwest of faults XVI and 12, and near a westward bend in fault M1. A small region between faults M1 and X has a chloride concentration between 250 and 500 mg/L (transitional water).

Marrufo (2011) noticed that clay volume in the Hueco Bolson increases toward the east from well 601 to 610 (Fig. 3), creating hydraulic anisotropy and impeding the flow of brackish water into the freshwater zone across faults IX and X. We also observed a similar trend between wells MNST5 and FBT04 (across fault IX; Fig. 13) and between MNST3 and MNST2 (across fault X; Budhathoki, 2013), where the volume of clay increases gradually to the east. In addition, Figures 11 and 13 suggest differing amounts and senses of offset of the erosional surfaces across faults IX and X that would impede the upwelling of brackish water from the east.

The high chloride concentrations to the west of fault 30 (Figure 15A) suggest that faults in this part of the study area, including fault 16N, may be acting as conduits for saline water, with fault 30 serving as a possible barrier to the flow of saline water to the east. Figure 11 indicates that fault 16N offsets all seven erosional surfaces by 18 m or more, suggesting penetration of faulting into the deeper aquifer where saltier waters are found. This region of higher salinity also is located where the East Franklin Mountains fault steps over to the west. Intensive fracturing often occurs within fault step-overs, and this may have enhanced hanging-wall permeability, allowing upwelling of saline water from deeper sources. Alternatively, Druhan et al. (2008) suggested that the East Franklin Mountains fault may have significantly uplifted the Fort Hancock Formation in this region, and that the increase in salinity could be due to halite dissolution from lacustrine sediments within the Fort Hancock Formation. The residual Bouguer gravity anomaly map (Fig. 6) shows a high that corresponds to the region of higher salinity, consistent with possible uplift of the Fort Hancock Formation. However, this does not explain the presence of pockets of more saline water to the east that are not associated with this anomaly high.

Bends in the XI-22-XVI fault system could be causing leakage of brackish water to the east from the deeper aquifer. The pocket of brackish water associated with fault M1 appears to be related to a single well and is not observed at a depth of 488 m (Fig. 15B). The transitional water found between faults M1 and X could be related to southward leakage of brackish water along two minor (<2.5 km long) faults mapped by Marrufo (2011).

Figure 15B shows that freshwater at 488 m depth is restricted to the region of the study area lying generally between fault M1 and the East Franklin Mountain fault, where the residual Bouguer anomaly (Fig. 6) is <–4.7 mGal, indicating the deepest part of the basin. The region of higher salinity observed south of fault XVI at 135 m depth (Fig. 15A) has now increased in size, and its center is offset ∼2 km to the west (Fig. 15B), suggesting that the shallower pocket of brackish water may be coming from a different source. All water north of 31.95°N exceeds a chloride concentration of 1000 mg/L. This northern boundary to the freshwater may be controlled by a basement structure that influences both the step-over in the East Franklin Mountain fault and the change in strike of intrabasin faults to the east.

The location of freshwater and brackish water within the northern Hueco Bolson is controlled both by stratigraphic and structural changes. Faulting near the Franklin Mountains in the northwestern Hueco Bolson appears to have fractured and/or uplifted deeper basin sediments, leading to upwelling of saline-rich water, despite its location close to sources of recharge from the mountains. A deeper structural feature affects the localization of faulting across the entire northern subbasin, causing a westward step-over in the East Franklin Mountains fault and a change in the strike of intrabasin faults (from north to northwest) to the east. This change in fault strike likely causes more fracturing within the deeper subsurface, allowing the migration of saline water from both sides of the basin into the northern portion of the aquifer.

Not surprisingly, the greatest thickness of freshwater is found in the deepest part of the basin between the East Franklin Mountains fault and fault M1, as indicated by the residual Bouguer gravity anomaly (Fig. 6). However, stratigraphy also plays a role in the localization of freshwater within the upper part of the aquifer, where an increase in clay volume across faults X and IX serves to create a barrier to flow of saline water from the east.

Based on our comparisons of gravity data, water well logs, and faults that have been mapped at the surface in the less urbanized parts of our study area, we were able to identify intrabasin faults with moderate amounts of displacement (100–200 m) of upper or lower basin fill. A comparison of simple Bouguer gravity anomaly changes to changes in elevation between gravity stations seems to be a quick first-order step in recognizing smaller intrabasin faults. However, this interpretation technique can be misleading if the gravity stations are too far apart, the associated data collection errors are large, or the data from several surveys with different accuracies are combined. Several of these intrabasin faults were also imaged using the more computationally intensive HGM technique on residual Bouguer gravity anomaly data and by constructing structural cross sections based on well logs.

This research was partly supported by a Terra Nova scholarship from the West Texas Geological Society to P. Budhathoki. A. Ruiz, from El Paso Water Utilities, provided well logs for this study. We thank G. Kaip, N. Mankhemthong, F. Ziwu, V. Avila, M. Sandoval, M. Lucero, and A. Gebregiorgis for assistance with the field work. S. Marrufo also provided gravity data used in this study. We thank M. Baker for discussions related to well-log analysis, and C. Montana for handling computer problems. We appreciate comments from an associate editor and reviewer that helped improve this analysis.