{'fraction_non_alphanumeric': 0.0258339210615678, 'fraction_numerical': 0.012068019747668678, 'mean_word_length': 5.5667238421955405, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaporite karst geohazards in the Delaware Basin, Texas: review of traditional karst studies coupled with geophysical and remote sensing characterization

*staffordk@sfasu.edu Citation:


INTRODUCTION
Permian evaporite karst of the Gypsum Plain in southeastern New Mexico and west Texas is extensive, including manifestations that range from simple, shallow hypergene caves to complex, hypogene caves that may reach significant depths. Although these regions are sparsely populated because of the arid environment and limited surface water resources, recent advances in hydrocarbon extraction techniques, primarily hydraulic fracturing (fracking) technologies, has resulted in increased infrastructure development and utilization within the region. With increased oil and gas development, karst geohazards are routinely encountered that range from subsurface instability during drilling operations and void interception during pipeline construction to road subsidence. Road failures are of critical concern because of potential economic impact associated with traffic delays and property damage; therefore, identification and delineation of karst geohazards has become an increasing priority within areas of intense petroleum International Journal of Speleology, 46 (2), 169-180. Tampa, FL (USA) May 2017 exploration and development in these karstified evaporite terrains. Common geohazards proximal to infrastructure within the Gypsum Plain include road subsidence associated with soil piping into solutional conduits / fractures, collapse failure of shallow karst features, suffosion induced by buried utilities, and suffosion enhanced by road construction (Fig. 1).
Traditional, surface-based analyses of karst development provide the most accurate assessment of speleological processes because visual observation of void development and morphology provide context for regional karst evolution. However, traditional surveys often provide limited information on extent and density of subsurface karst phenomena, because surficial karst manifestations often represent only a small fraction of total karst development. Therefore, coupling of traditional surveys with non-invasive, remotely-sensed characterization of karst phenomena can be used to better delineate potential regions of increased geohazard risk both locally and across broad regions. In this study, traditional karst surveys (i.e., geomorphic surface and cave mapping) are coupled with geophysical techniques (i.e., resistivity analyses) and aerial remote-sensing techniques (i.e., digital elevation models and image analyses) to better characterize gypsum karst development, spatial distribution, and the relationship to potential geohazards associated with thoroughfares heavily impacted by increased traffic associated with oil and gas operations. Karst geohazards associated with infrastructure throughout the Gypsum Plain: A) road failure as a result of suffosion; B) sinkhole collapse associated with increased suffosion into solution conduit at margin of road; C) vertical suffosion structures associated with buried utility cable (fiber optic cable was observed at 2.5 m depth within feature); and D) increased suffosion adjacent to concrete reinforcement at pavement shoulder (i.e., toe wall). White and black scale bars in figures are ~50 cm long.

GEOMORPHIC AND GEOLOGIC ENVIRONMENT OF THE GYPSUM PLAIN
The evaporite karst region of southeastern New Mexico and west Texas, commonly referred to as the Gypsum Plain (Hill, 1996), encompasses ~2800 km 2 that is dominated by outcrops of the Ochoan Castile Formation and to a lesser extent residual portions of the Salado and Rustler formations in Eddy County, New Mexico and Culberson County, Texas, within the Delaware Basin (Fig. 2). The Gypsum Plain is located along the northern edge of the Chihuahua Desert where annual precipitation averages 267 mm with most precipitation occurring as late summer monsoonal storms (Sares, 1984); however, anomalously high precipitation events can occur at decadal intervals as single storm events that are less than 48 hours in duration yet exceed 100 mm. Annual temperature averages 17.3°C, with an average annual low and high of 9.2°C and 25.2°C, respectively (Sares, 1984).
The Delaware Basin was formed as an intracratonic basin with limited connectivity to the open marine environment as a result of assimilation of Pangea during the Late-Paleozoic (Dickenson, 1981). During Guadalupian time, the Delaware Basin was dominated by peripheral development of the Capitan Reef (Fig. 2) while mixed evaporite / carbonate / clastic sequences were deposited on the surrounding platform and shelf environments; the interior of the Delaware Basin during the Guadalupian was dominated by episodic, siliclastic deposition associated with sealevel lowstands (Scholle et al., 2004). During the Ochoan, the Delaware Basin became increasingly restricted and the Capitan Reef carbonate factory shut down as the basin interior transitioned into a density-stratified, hypersaline sea that was infilled by Castile strata (Scholle et al., 2004). Subsequent to basin-infilling, mid to late Ochoan deposition of the evaporite-dominated Salado Formation and mixed evaporite / carbonate strata of the Rustler Formation capped the region, including the Delaware Basin and surrounding shelf and platform environments (Scholle et al., 2004). These evaporite units provided regional seals that promoted isolation, development, and accumulation of petroleum resources across the greater Permian Basin that are exploited today.
Throughout the Mesozoic and Cenozoic, the region was largely subaerially-exposed with siliciclastic deposition and erosion dominating, although a brief period of carbonate deposition occurred during transgression of the Cretaceous Western Interior Seaway (Hill, 1996). Pangea rifting exhibited minimal effect on the Delaware Basin region, but Late Mesozoic compression associated with the Laramide Orogeny tilted strata 3-5° to the east / northeast and uplifted the region significantly above sea-level, where it remains today (Horak, 1985). Igneous activity of the Trans-Pecos magmatic province elevated the regional geothermal gradient during the late Paleogene including emplacement of igneous dikes throughout the northern Delaware Basin. However, most brittle deformation associated with Basin and Range extension in the area likely did not occur until early Neogene (Hentz & Henry, 1989); conjugate joint / fault sets oriented at ~N75°E and ~N15°W induced by extension are common throughout the Gypsum Plain (Nance, 1993). Significant Quaternary climate change sculpted the modern geomorphic surface of the Gypsum Plain as the region transitioned from an early Pleistocene cool, wet climate to the dry, arid conditions that dominate the region today (Hill, 1996).

TRADITIONAL KARST ANALYSES OF THE GYPSUM PLAIN
Karst development within the greater Delaware Basin is common and widespread, most notably, the world famous carbonate caves of the Guadalupe Mountains (e.g., Carlsbad Cavern, Lechuguilla Cave), but also cave development in mixed carbonate / evaporite strata of the Northwestern Shelf (e.g., Coffee Cave) (Stafford et al., 2008a) and evaporite karst development within the Delaware Basin (e.g., Gypsum Plain) (Stafford et al., 2008b). Although karst phenomena have been documented in the Salado and Rustler formations (Stafford et al., 2009), the most extensive evaporite karst development documented in the Gypsum Plain is within the Castile Formation outcrop region and thus is the emphasis of this study. Traditional karst surveys have been conducted by the authors throughout the Gypsum Plain that include geomorphic surface mapping and delineation of karst phenomena based on systematic traverse surveys. Numerous caves large enough for human entry have been surveyed and mapped to delineate their spatial extent and speleogenesis; analyses indicate that 30-40% of enterable caves are of hypogene origins, including the majority of caves greater than 50 m in length. Stafford and others (2008c) showed that surficial manifestations of Castile karst development within the outcrop region are highly clustered with more concentrated development in the western, and specifically the northwestern, portions of the outcrop region. Karst phenomena within the Castile outcrop region can be broadly classified into four categories: 1) Intrastratal Dissolution; 2) Hypogene Caves; 3) Hypergene Caves; and 4) Suffosion Caves. The Castile Formation is primarily a varved sulfate / carbonate unit up to 480 m thick in the subsurface of the eastern Delaware Basin (Hill, 1996), where interbedded halite also occurs, but thins to a solution edge on the updip western margin proximal to the Delaware Mountains (Fig. 2). Halite interbeds have been removed by intrastratal dissolution in shallow subcrop and outcrop regions within the Gypsum Plain; sulfate / calcite laminae reflect annual variation in salinity at time of deposition related to dry / wet seasons, respectively (Scholle et al., 2004).
Castile hypogene karst can be subdivided into widespread intrastratal dissolution and focused hypogene caves. Surficially, intrastratal dissolution manifests as both negative and positive relief structures, forming either solution grabens and subsidence troughs or resistant brecciated hills (Stafford et al., 2008b,d). Widespread intrastratal dissolution has removed halite interbeds from the Castile outcrop and subcrop region as solutionallyaggressive, ascending fluids delivered from underlying siliciclastic strata (Guadalupian strata of the Bell and Cherry Canyon formations) migrated updip towards the west in the past (Stafford et al., 2008b, d). This International Journal of Speleology, 46 (2), 169-180. Tampa, FL (USA) May 2017 resulted in differential collapse and settling across the region including development of laterally-extensive, blanket breccias that provide high permeability flow paths for migration of fluids (Stafford et al., 2008b). In more intensely affected regions, solution grabens up to 30 m deep have developed while subsidence valleys infilled with Cenozoic clastics are common proximal to persistent fluvial bodies (Maley & Huffington, 1953).
Throughout the Gypsum Plain, and greater Delaware Basin, breccia pipe development commonly results in surficial manifestations of resistant topographic highs. Similar to blanket breccias, solutionally-aggressive fluids sourced from underlying siliciclastics have created vertical stoping structures that can extend hundreds of meters laterally and vertically that are commonly affiliated with fractured regions where preferential vertical pathways enabled fluid migration, both ascending, low-density, undersaturated fluids and descending, high-density, saturated fluids (Anderson et al., 1978;Anderson & Kirkland, 1980). Often these breccia pipes become re-cemented as fluid flow paths change with time and thus create masses that are solutionally more resistant during surface denudation (Fig. 3A). However, it is not uncommon for these features to retain vertical, high permeability zones, especially at the distal breccia zone margins (Fig. 3B).
Hypogene caves occur as a continuum between more classically defined, multi-storey maze caves hundreds of meters in survey length in strata that are complexly folded and fractured to single, riser features that can extend to depths of 100 m in poorly fractured rock (Stafford et al., 2008b). All hypogene caves exhibit the morphometric suite of rising flow features that have been defined as characteristic of hypogene speleogenesis (Klimchouk, 2007) and show little lithologic control (Fig. 3D). Hypogene caves cut across lithologic boundaries with regularity, showing little to no preferential development with transitions between varved gypsum, recrystallized (selenite) gypsum or other secondary sulfate fabrics. Most hypogene caves occur proximal to the underlying contact with siliclastic strata, where ascending, solutionallyaggressive fluids are sourced (Stafford et al., 2008b); therefore, documented hypogene caves occur in greater frequency proximal to the western margin of the Gypsum Plain where surface denudation, and subsequent cave breaching, has been greatest and total thickness of the Castile Formation has been greatly reduced. It is probable, that significant hypogene karst development continues near the lower contact of the Castile Formation further to the east, but surface denudation has not proceeded far enough to induce system breaching to enable human exploration. Anecdotal reports by drilling operations in the region commonly report bit drops near the base of the Castile Formation, which further supports this theory.
Hypergene karst is disseminated broadly across the Gypsum Plain but commonly occurs in higher concentrations where rock is more intensely fractured (Stafford et al., 2008b). Hypergene caves in the western portion of the outcrop area exhibit strong structural control along fractures, are laterally-limited and have rapid average aperture decreases away from cave entrances, suggesting that overland flow to these features is near saturation when entering the subsurface and has little remaining solutional potential. While most of these caves appear to be dominantly-formed by vadose processes (Fig. 3E) in the west, phreatic tubes (Fig. 3C) become increasingly common in the eastern portions of the Gypsum Plain. Often these phreatic caves show greater lengths, occasionally more than 200 m, that transition from partially-filled to completely-filled with water. They are often heavily armored by clay-rich mud suggesting that undersaturated fluids are able to travel greater distances before reaching saturation, thus enabling significantly greater hypergene conduit development. The presence of hypergene conduits at relatively shallow depths (less than 20 m) within ridges that are in close proximity to solution troughs more than 20 m deep indicate that Gypsum Plain hydrogeology is highly partitioned.
Bedrock exposures are common across the Gypsum Plain, but many regions are mantled by mixed alluvial gravels and gypsic soils that have resulted from the transition from the cool, wet climate of the Pleistocene to the hot, dry climate of today (Hill, 1996). Extensive surveys on the extent of solutional valley and sinkhole fills have not been conducted across the Gypsum Plain, but excavations by the authors have revealed that gypsic soil thickness can vary from decimeters thick to more than five meters thick over lateral distances of less than ten meters. Within these soils, preferential soil piping and solution of gypsum fractions produces lateral and vertical pipes. Often, soil piping is connected with solutional conduits at depths and more commonly with solutionally-widened fracture zones (Fig. 3F). Therefore, soil caves form by a combination of suffosion and solution of soluble fractions, often with little correlation to surrounding surface geomorphology. Most of these caves are less than five meters in length with heights under one meter, but when associated with larger solutional caves that provide outlets for suffosion products, soil caves may extend for tens of meters and even create "soil chambers" over five meters tall (Fig. 3G).
Further complicating karst development within the Gypsum Plain are varying degrees of diagenetic overprinting, most notably evaporite calcitization (Kirkland & Evans, 1976). Ascending light hydrocarbons from deep sources have migrated upwards through the Castile Formation providing the source material for sulfate reduction (Lee & Williams, 2000); however, the exact mechanism of sulfate reduction is still debated in the region as being either bacterial sulfate reduction (BSR) or thermal sulfate reduction (TSR) (Stafford et al., 2008d). As a byproduct of sulfate reduction, evaporite calcitization is widespread but most frequently associated with breccias. In many areas, calcitized zones also host native sulfur deposits, some of which have been economically mined in the past (Phillips, 1917;Wallace & Crawford, 1992), that attest to a diagenetic environment that was restricted and prevented migration of hydrogen sulfide byproducts out of the system. Hill (1996) suggested that secondary selenite masses within the Castile Formation resulted from oxidation of native sulfur ore bodies as flow regimes transitioned within the Delaware Basin and enabled oxygen-rich waters to migrate into sulfur bodies; this process was likely the result of hypergene processes altering mineralized zones created by hypogene processes within the Gypsum Plain.

AIRBORNE-/ SATELLITE-BASED, REMOTELY-SENSED CHARACTERIZATION
Recent advances in airborne-and satellite-based remote sensing are now enabling rapid, large-scale characterization of geologic variability and karst phenomena as remote-sensing technologies rapidly improve. Traditionally, airborne photogrammetry has provided imagery for remote characterization of geomorphic surfaces and vegetation patterns; however, computer-based, geospatial analyses have enabled this data to be more easily accessed, analyzed and accurately georeferenced to surface locations. Color InfraRed (CIR) imagery has enabled variations in moisture content and vegetation to be more easily assessed (Yang et al., 2013), while Light Detection and Ranging (LiDAR) methodologies continue to expand the abilities of airborne platforms, where pulses of ultraviolet light are used to accurately map elevation changes (Liu, 2008). New satellite platforms offer exciting new possibilities for karst terrain assessment; however, the resolution of commercially available hyperspectral and multispectral data is currently insufficient beyond broad scale characterization.

Methodology for Airborne-/ Satellite-Based, Remotely-Sensed Data
Both CIR and LiDAR data used in this study were collected by TxDOT (Texas Department of Transportation) as part of geohazard assessment of infrastructure within northern Culberson County, Texas. CIR and LiDAR data were collected with a minimum horizontal resolution of 50 cm and 10 cm vertical resolution; data coverage included more than 55 km of 300-meter-wide swaths centered along roadways (i.e., 150 m normal to roads) in the study area (~18 km 2 in total area). Data was collected with a horizontal datum of NAD83 (2011) and vertical datum of NAVD88 and geoid model GEOID12A. A surface adjustment factor of 1.00025 was applied to all data to adjust for flight path in relationship to imaged position.
LiDAR data was processed in ArcGIS with 64.76% of data representing bare earth or ground sampling; average horizontal resolution of processed, bare earth LiDAR data was 0.3 m. For karst delineation, a digital elevation model (DEM) was created and topographic lows calculated through raster subtraction of original DEMs from filled DEMs to delineate sinkholes and  (Stafford et al., 2008c). The data was filtered to remove any features <10 cm deep because that data could not be differentiated from potential noise within vertical data resolution; resulting features delineated were smoothed and merged where they exhibited overlap and likely represented individual closed depressions. This produced a total of 8,349 closed depressions that were filtered based on geology and anthropogenic structures identified through CIR image analyses to remove features that exhibited low probability of being related to karst development in the region; 5121 karst sinkholes remained after filtering of DEM, producing an average density of 300 sinkholes / km 2 with more than 90% of features being less than 5 m 2 is size. Additionally, bare earth LiDAR return intensity was used as a proxy for assessment of geology through intensity classification (Bryant, 2012).
CIR imagery was used to characterize natural/versus anthropogenic features throughout the study area. CIR imagery was also used to identify trends within the karst landscape related to vegetation patterns; lineaments representing planar brittle deformation features (joints and fractures) were digitized based on expression of vegetation patterns. Similarly, vegetation patterns were used to assess relative moisture content within gypsic soils to identify regions of enhanced shallow fluid migration. Normalized Difference Vegetation Index (NDVI) was calculated from CIR imagery to identify the relative "health" of vegetation (Yang et al., 2013), which is directly related to available moisture in arid desert environments like the Gypsum Plain; regions with a high NDVI index are interpreted as regions of greater moisture flux and probable sites of increased karst development. Figure 4 provides comparative results of GIS-based analyses of LiDAR and CIR imagery data processing at two example sites (Fig. 2) within the study area, including: 1) CIR imagery; 2) NDVI analyses of CIR imagery; 3) DEM derived from LiDAR analyses; 4) geologic characterization based on intensity classification of LiDAR data; and 5) composite image of karst and structural data attained from GIS analyses compared with karst features physically documented during field mapping. CIR imagery clearly shows the sparse vegetation that is common throughout the Gypsum Plain with major sinkhole arroyos easily discernable as is expected from high-resolution orthoimagery. NDVI analyses, although designed as an index for vegetation, clearly shows spatial patterns that are different than what would be predicted based on visual interpretation of CIR imagery, but instead indicate regions where "healthy" vegetation is able to exploit fractures and karst development in order to access shallow groundwater. Therefore, regions that exhibit high values of NDVI (values close to 1) can be used in the arid, desert environment of the Gypsum Plain as a proxy for identification of probable karst development based on increased moisture content associated with these features; caution must be used, however, with NDVI analyses because areas with strong "shadow zones" in CIR imagery that occur along high angle scarps / sinkholes provide false positives as can be seen at Study Site 2 when compared with the DEM of the same location (Fig. 4).

Interpretation of Airborne-/ Satellite-Based, Remotely-Sensed Data
Digital elevation modeling from high-resolution LiDAR analyses within the study area clearly identifies topographic lows and highs (Fig. 4), both of which are related to karst development across the Gypsum Plain. Small depressions developed along an arroyo oriented southwest / northeast are visible at Study Site 1, while large-scale sinkholes with significant depth occur at Study Site 2, including both sinkhole development with well-defined, dendritic arroyos (eastern edge of Study Site 2) and near-vertical collapse structures (western edge of Study Site 2). Additionally, raised topographic ridges associated with near-surface hydration along fractures can be seen as roughly north / south linear topographic highs in Study Site 2. While not traditionally used as a method for geologic mapping, LiDAR intensity can be used for classified image analyses throughout the Gypsum Plain because of the high density of ground returns and low spatial density of vegetation. At both Study Site 1 and Study Site 2, intensity classification was used to delineate between gypsum bedrock, gypsic soil (gypsite), siliceous-rich alluvium (silica largely eolian derived) and paved road surfaces; however, as with NDVI analyses, caution must be used with LiDAR intensity classification because regions of lower density LiDAR ground returns associated with larger vegetation (i.e., trees) that create high absorption appear similar to values of paved road surfaces.
As a result of image analyses described above and field verification during physical land surveys throughout the study area, composite maps (Fig. 4) that delineate the extent of sinkhole development, spatial distribution of fractures, occurrence of surficially-expressed karst features and locations of probable karst development throughout the Gypsum Plain were constructed. High-resolution geomorphic maps, including surficial landscape composition, were derived from LiDAR intensity analyses and regions that likely exhibit additional near-surface karst developed were delineated based on moisture variability derived from NDVI analyses. Therefore, spatial imagery analyses throughout the study area indicate that both traditional (i.e., CIR and DEM analyses) and nontraditional (i.e., NDVI and LiDAR intensity analyses) remote sensing can provide powerful tools for karst and associated potential geohazard characterization in arid environments with low vegetation densities like those that commonly occur in evaporite karst terrains like the Gypsum Plain.

LAND-BASED, GEOPHYSICAL CHARACTERIZATION
Geophysical methods provide land-based, noninvasive mechanisms for characterization of geohazards within karstified terrains, but methods vary significantly based on depth and resolution of karst anomalies that are investigated. Brown et al. Images from top to bottom, in descending order, include: 1) CIR imagery; 2) NDVI analyses of CIR imagery used to asses probable karst development based on moisture variability; 3) DEM derived from LiDAR analyses; 4) geologic map derived from classification of LiDAR intensity; and 5) composite map of karst features identified through combined field surveys and image analyses. All images are 300 m wide and oriented with north towards the top of the page.
(2011) conducted a comparative study on geophysical techniques for subsurface porosity characterization in Ordovician carbonate strata in central Texas. In their study, they compared the effectiveness of microgravity, ground penetrating radar (GPR), direct current (DC) resistivity, capacitively coupled (CC) resistivity, induced polarization (IP) and ground conductivity (GC). Based on their analyses, GPR provided the highest resolution of karst phenomena at shallow depths (<40 cm) and microgravity provided the most suitable data for delineation of large-scale trends. However, resistivity (both DC and CC) proved the best source for imaging relatively-shallow karst phenomena at depths ranging from approximately one meter to ten meters (Brown et al., 2011), which falls within the average range of most karst geohazard concerns. Gypsum karst geohazards at depths of less than one meter are commonly breached by natural denudation processes and identifiable by traditional karst surveys or will likely be intercepted by construction activities during anthropogenic modification associated with infrastructure development; karst phenomena at depths greater than ten meters likely exhibit low probability of catastrophic failure in evaporite karst terrains. Therefore, for this study both DC and CC resistivity analyses were conducted to evaluate potential shallow karst phenomena.
Direct Current (DC) resistivity analyses were conducted using an AGI (Advanced Geoscience Inc.) SuperSting R8/IP multi-electrode resistivity meter with 56 electrodes spaced at 1, 2, or 4 m in a dipoledipole array. DC resistivity analyses were conducted at 16 sites adjacent to roads within the study area where shallow geohazards were suspected based International Journal of Speleology, 46 (2), 169-180. Tampa, FL (USA) May 2017 on road subsidence. Most sites were investigated with at least two different electrode spacings, where larger spacing increased depth of investigation but at reduced resolution. DC resistivity data was reduced using AGI's EarthImager 2D and a smoothed model inversion. Data with a high degree of misfit to predicted modeled parameters were edited and terrain corrections were applied. DC analyses were conducted in late spring 2016 when soil moisture was slightly elevated due to low-intensity spring precipitation events. Figure 5 provides representative DC resistivity inversions with 4 m electrode spacing and an approximate depth of resolution of 20 m for data collected at Study Site 1 and Study Site 2 (Fig. 2). These sites illustrate the most common karst phenomena occurrences documented through resistivity analyses; data was acquired parallel to roads in regions where visible road failure had occurred in the past. At Study Site 1, thin gypsic soil (~1-2 m thick) covers gypsum bedrock that exhibited high moisture content (saturated gypsum) in the upper regions at the time of survey that transitions into low moisture content (unsaturated gypsum) with increasing depth; small solution conduits were documented at shallow depths that appear largely filled with moisture-rich sediment although at least two resistivity anomalies appear to be either partially-filled or air-filled conduits. At Study Site 2, thin gypsic soils also occur over bedrock with variable moisture content that is more irregular than at Study Site 1. However, karst development at Study Site 2 largely consists of "leached bedrock" regions that are likely associated with preferential dissolution along calcite / gypsum laminae in the Castile Formation. At both Study Site 1 and Study Site 2, near-vertical fractures indicate regions of increased fluid migration and probable solutional widening.

Methodology for Land-Based, Geophysical Characterization
Capacitively Coupled (CC) resistivity analyses were conducted with a Geometrics TR5 OhmMapper in a dipole-dipole array configuration composed of five receivers connected by 2.5 m coaxial cables and a transmitter offset of 2.5 m. Data was collected at a transmission rate of once per second with a traverse velocity of approximately one meter per second. GPS data was simultaneously collected, using a Trimble Nomad 900 series logger connected to a Pathfinder Pro receiver and Zephyr antennae with a horizontal accuracy of less than 50 cm. CC resistivity analyses included more than 120 km of collected data along thoroughfares in the study area that was processed using AGI's EarthImager 2D software to produce smooth model inversions consistent with DC resistivity analyses conducted at proximal sites. CC analyses were conducted in mid-summer 2016 with relatively dry conditions and high temperatures. Figure 6 provides representative CC resistivity inversions with an approximate depth of resolution of four meters for data collected at Study Site 1 and Study Site 2 (Fig. 2); these data were cropped from continuous data collected and processed in multikilometer segments. While the data was collected proximal to data represented in DC resistivity analyses (Fig. 5) presented above, CC data was acquired directly within road paths for ease of continuous resistivity data collection, thus spatial variability is seen in comparison with DC data analyses that is further exacerbated by drier near-surface sediments / strata at time of CC data collection. However, results of CC data analyses indicate similar karst phenomena occurrences at Study Site 1 and Study Site 2, including shallow bedrock / soil contacts, solutional conduits and "leached bedrock." At Study Site 1, shallow solutional conduits appear as higher resistivity anomalies in CC resistivity data likely resulting from drying of sediment fills that appear moisture-rich in DC resistivity analyses. Leached zones at Study Site 1 appear better defined in CC analyses as a result of the higher resolution data presented here; however, the lack of depth of resolution precludes interpretation of deeper manifestations of karst development observed in DC resistivity analyses. DC resistivity tomography consistently correlates with CC resistivity analyses when karst heterogeneity is considered in relation to the location of DC and CC surveys that were conducted up to 10 m apart, although parallel to each other.

Interpretation of Land-Based, Geophysical Data
Resistivity studies within the Castile Formation proved an effective method for characterization of karst phenomena and identification of potential geohazards. However, data could not be accurately interpreted based on absolute resistivity values associated with specific karst or rock properties; instead, interpretations required analyses based on surrounding geologic characterization and environmental parameters at the time of data collection, including relative soil moisture and thermal fluctuations throughout the day. While extremely high resistivity values are generally associated with open void space in karst studies, they are also associated with low moisture regions in gypsic soils and poorlyfractured bedrock in this study. Similarly, extremely low resistivity regions are generally associated with water-filled voids or saturated media, including moist, soil-filled caves; however, it is difficult to differentiate these from regions of heavily-leached, moisture-rich bedrock and regions of increased capillary fluids within vadose soils in this study area. Therefore, it is critical that proximal geology and associated karst development assessed during traditional geomorphic and karst mapping be assimilated into interpretations, especially in arid, gypsic terrains that push instrumental limitations such as the environment of the Gypsum Plain. Interpretations made from resistivity tomography analyses at both Study Site 1 and Study Site 2 were confirmed through excavations at these sites (Fig. 7); consistently, field excavations confirmed resistivity tomography interpretations at sites throughout the study area, but detailed knowledge of the surrounding geology proved essential for proper interpretations of geophysical anomalies identified.

CONCLUSIONS
Evaporite karst within the Gypsum Plain has the potential to create significant geohazards; however, coupling of traditional karst studies with remotelysensed geomorphic / geologic analyses can provide an effective toolset to predict and proactively mitigate potential regions of thoroughfare failure. The two example sites used here illustrate the effectiveness of coupling traditional and remotely-sensed karst analyses. In the study area, the authors effectively coupled these methodologies over larger areas, including: 1) ~14 km 2 of traverse-based physical karst surveys; 2) 110 km of continuous, capacitively coupled resistivity analyses; 3) >4000 m of direct current resistivity analyses; 4) ~18 km 2 of LiDAR analyses; and 5) ~26 km 2 of CIR analyses.
The results of these analyses delineated 5121 sinkholes through LiDAR analyses, while only 396 open karst features (i.e., non-filled sinkholes) were identified with physical land surveys in the same area. LiDAR analyses proved to significantly better delineate shallow, filled sinkholes than traditional surface surveys; traditional surface surveys better delineated open cave and karst features because physical field observations enable delineation of features with orifices less than 30 cm in diameter which is below the accuracy of LiDAR data within this study. CIR imagery provides invaluable data on structural controls within evaporite terrains through lineament analyses, while Fig. 7. Images of representative field excavations of karst features: A) typical "leached bedrock" zone; B) suffosion cave filled with moisture-rich alluvium developed in gypsic soil; C) small, shallow cave detected adjacent to road (scale bar ~1 m); D) solutional conduit at base of excavation (note backhoe excavation was conducted during rain event which is why there is water / mud in trench); E) typical gypsum epikarst development beneath shallow soils (note that this site was excavated in the past as a borrow pit but illustrates well the irregular nature of soil / rock contacts within the study area).
NDVI analyses delineate regions where moisture flux is highest and thus regions where karst processes are likely to be currently active and additional nearsurface karst development is probable. Analyses of LiDAR return intensity provided promising results for mapping of geologic variability within evaporite karst terrains where vegetation densities are low and can be used to supplement existing geologic maps of study regions.
Resistivity surveys identified more than 400 regions of karst development, including 181 individual solutional conduits, 133 zones of solutionallywidened fractures, and 104 individual regions of bedrock leaching. DC resistivity data provided greater depth of resolution but proved time intensive, while CC resistivity data could be acquired relatively rapidly but at decreased depth of resolution. Interpretations of both DC and CC resistivity tomography proved consistent when temporal and spatial variability of data was considered. In this study, field excavations at DC resistivity sites were used to verify interpretations of proximal DC and CC resistivity data, which were then used as "training data" for interpretation of CC resistivity tomography collected throughout the study area.
Traditional karst studies under predicted total karst development within the study area when compared with remotely-sensed analyses; shallow, filled depressions were not readily discernable in the field. Geophysical methods could only be accurately interpreted with correlation to field investigations, both proximal cave studies and excavations, because resistivity signatures exhibited significant variability based on variations in soil composition, structural deformation, and nature of solutional development. Therefore, coupled airborne-/ satellite-based and geophysical remote sensing provided an effective means of delineating / characterizing geohazards to facilitate mitigation plans when coupled with traditional karst studies. Ultimately, this composite toolset proved essential for efficient and effective characterization of potential karst geohazards within the study area; excavations of predicted potential geohazards have confirmed these findings at multiple sites throughout the investigated region. While some of the remote-sensing techniques used in this study were novel (e.g., NDVI, LiDAR intensity classification) in their use for delineating karst phenomena, they proved effective throughout the sparsely-vegetated Gypsum Plain and provide promising results that are likely to produce similar results in other evaporite karst terrains.