Airborne geophysical surveys of the lower Mississippi Valley demonstrate system-scale mapping of subsurface architecture

A system-scale airborne geophysical survey

From 2018 through early 2020, we acquired more than 43,000 flight-line-kilometers (line-km) of airborne geophysical data over the MAP study area of ~140,000 km² (Fig. 2a, “Methods” section). Data collection included a high-resolution survey over ~1000 km² near Shellmound, Mississippi, regional surveys with 3–6 km line spacing across the entire study area, and over 3000 line-km of data acquired along streams and rivers to characterize potential surface water–groundwater connections beneath these important recharge pathways. Radiometric (Fig. 2b), magnetic (Fig. 2c), and inverted resistivity grids at multiple depth intervals (Fig. 2d–h) summarize the combined results from both regional survey phases. Together, this represents the first initiative to acquire system-scale airborne geophysical data over an entire US aquifer.

At the regional scale, radiometric data (Fig. 2b) correlate with mapped surficial geology³² and sediment age, with Holocene deposits clearly delineated as strong returns of multiple elements (light-colored areas in Fig. 2b) compared with Pleistocene sediments. Magnetic data gridded at this scale largely corroborated previously mapped structures, such as the line of southwest–northeast-trending magnetic highs (Fig. 2c) associated with mapped intrusive plutons along the Commerce geophysical lineament (CGL) adjacent to the RR in northeast Arkansas and southeast Missouri^45,46.

Inverted resistivity models in the uppermost 5 m (Fig. 2d and Fig. S1b) correspond closely with mapped surficial units³² (Fig. S1a). For example, Wisconsin-age valley train deposits, as well as point bars and meander scrolls of modern river networks, appear as resistive features in the uppermost 5 m, whereas fine-grained units such as backswamp deposits are characterized by low resistivity. At 20–25 m depth (Fig. 2e and Fig. S1c), intermediate to high-resistivity values are consistent with the coarse-grained lithology found throughout the MRVA. Lower resistivity can be found at this depth in sedimentary units outside the MAP region as well as over structural highs within the MAP region where Tertiary units are close to the surface. Low-resistivity values show the Vicksburg–Jackson confining unit (VKBG) in southeast Arkansas and northeast Louisiana, the Midway confining unit (MDWY) in northeast Arkansas, and an erosional remnant of Tertiary sediments beneath Crowleys Ridge. Beneath the Quaternary MRVA (~30–50 m depth in most of the region), resistivity values strongly correlate with Tertiary subcropping units (Fig. 2f–h and Fig. S1d, e). Most notable here is the low resistivity that corresponds with the known regional VKBG and MDWY confining units, both clay and shale rich. In contrast, the subcrop of the Middle Claiborne aquifer (MCAQ) sands are resistive (Fig. 2g, h and Fig. S1d), with a notable change in facies north of Memphis, Tennessee, associated with the coarse Memphis Sand of the Claiborne Group (dashed line, Fig. 2g, h).

Further correlation between resistivity and geologic structure is evident in cross-section view (Fig. 3a–c), where gridded resistivity models are compared with the top elevation of MERAS model surfaces⁴³. From west to east, prominent features in Fig. 3a include the highly resistive Paleozoic Ozark Plateaus aquifer that bounds the MAP region, the conductive east-dipping MDWY beneath and west of Crowleys Ridge, and the resistive MCAQ dipping to the east beneath the east side of Crowleys Ridge. This section highlights the variable degree of connectivity between the MRVA and underlying Tertiary aquifers. While the MCAQ sands appear connected to the shallow MRVA (veneer of moderate to high resistivity in the upper 30–50 m) immediately east of Crowleys Ridge, these aquifer layers become mostly separated by the Middle Claiborne confining unit (MCCU) farther east. To the south (Fig. 3b), the subcrop of the Claiborne Group aquifers (MCAQ and Upper Claiborne aquifer) suggests a direct connection to the MRVA beneath the Mississippi Delta, sandwiched between the VKBG confining unit in the shallow western half of the section and the Lower Claiborne confining unit (LCCU) that dips westward in the eastern part of the section. Towards the southern end of the AEM survey (Fig. 3c), the MRVA is disconnected from Tertiary aquifers by the MCCU and VKBG confining units in the western and eastern portions of the section, respectively. AEM-derived resistivity models largely corroborate the existing framework but also reveal greater detail in the overall structure and heterogeneity within units than could be previously determined through relatively sparse borehole observations. This study demonstrates that systematic mapping at high spatial resolution with AEM data illuminates model structural details expected to exist throughout the region but that cannot be fully understood with sparse observations.

Combined resistivity models from the two phases of regional data collection (Fig. 2d–h and Fig. 3a–c) have been interpolated onto grids with 1 km × 1 km lateral and 5-m vertical resolution useful for investigation of regional-scale structure. However, the native resolution of resistivity models along flight paths is much higher, with spacing between sounding locations for the Resolve and Tempest AEM systems equal to 25 and 75 m, respectively, with finer-scale near-surface vertical resolution for the Resolve system. On a native-resolution cross-section of the Resolve system in northeast Arkansas (Fig. S2), significant detail reveals the internal structure and variability of the MRVA. For example, the Resolve models capture the topography of the aquifer base, including incised channels in the subcropping Tertiary unit that may have been formed during periods of glacial outwash (interpreted locations marked “C” on Fig. S2), and the thickness and extent of surficial low-resistivity material that may be local confining units and a barrier to recharge. Lateral transitions in the upper ~5–30 m correspond with mapped braid belts³³, suggesting that electrical resistivity can be an indicator of these distinct lithologic units. However, since the chronology of mapped braid belts is largely based on relatively shallow OSL dates and surficial mapping³³, there is little constraint on the age of deeper Quaternary sediments beyond the estimate of their deposition after 250 ka³⁴. Given the presence of shallow braid belt deposits older than the last sea-level lowstand ~20 ka³³, along with the observed internal structure of the MRVA with a discontinuous low-resistivity layer at ~30 m (Fig. S2), we hypothesize that the deeper Quaternary sediments may represent earlier filling of post-250 ka eroded channels of the ancestral Mississippi–Ohio River systems^34,38.

Derivative products: interpretations of hydrogeologic structure and properties

While lithology dominates AEM-derived resistivity values in the MAP region, porewater salinity is also known to influence resistivity^17,19,24, and groundwater salinity varies throughout the study area^47,48. Specific conductance (SC) measured in boreholes throughout the MERAS domain are generally low, with similar values in both Quaternary sediments (median log₁₀ SC 2.79 = 617 μS/cm) and deeper Tertiary units (median log₁₀ SC 2.67 = 471 μS/cm). Areas of high salinity are limited within the MRVA footprint, with only 6% of the area predicted to be >1000 μS/cm⁴⁸. Correlation between SC- and AEM-derived resistivity (Fig. 4) by MAP region (Fig. 2) suggests that SC has limited overall control on resistivity across the study domain, but can be important in specific regions where SC is high, generally in the Grand Prairie (Fig. 4a) and Boeuf (Fig. 4b) regions. Lithology appears to be a primary driver for resistivity, with Quaternary sediments typically higher in resistivity than Tertiary units (Fig. 4a–f) over the same SC range. Tertiary, and to a lesser degree Quaternary, units follow an SC-resistivity trend indicative of moderate-high surface conductivity caused by an increased fraction of clays or other fine-grained sediments that flattens the slope of this relationship (black curves, Fig. 4a–f), compared with Archie’s Law where surface conduction is absent^49,50. Notable exceptions where Quaternary and Tertiary resistivity values are similar occur in the Cache (Fig. 4c) and St. Francis (Fig.4d) regions, where coarse-grained MCAQ sand subcrop beneath the MRVA and appear similar geophysically (double-sided arrow in Fig. 3a just east of Crowleys Ridge). The sensitivity of AEM data to both model structure and porewater salinity makes it a powerful tool in hydrologic studies, and also highlights the need for borehole and other geologic observations to calibrate against to reduce uncertainty in the extrapolation of AEM interpretations across the entire survey domain.

Using 6130 published picks defining the depth to the base of the MRVA³⁹ within the AEM survey area, along with an additional 364 manual picks based on observation of resistivity cross-sections, we used a supervised machine learning algorithm⁵¹ (see “Methods” Section) to interpret the elevation of the base of the aquifer across the entire AEM dataset (Fig. 5a). An associated aquifer thickness map (Fig. 5b) is estimated by differencing the base of aquifer elevation surface from the land surface elevation, and saturated thickness (Fig. S3) is calculated by differencing the base of aquifer elevation from the 2018 potentiometric surface, assumed here to be the water table, derived from borehole observations⁵². By subregion (Fig. 2a and Fig. S3), saturated thickness is greatest in the St. Francis region where deep Quaternary scour channels of the Mississippi–Ohio River system have been documented east of Crowleys Ridge³⁸, and thinnest in the Grand Prairie region because of the combined influence of deep water table, a thick confining layer that limits recharge, and shallow subcrop of the VKBG (Fig. 2e). The difference in saturated thickness between AEM and borehole interpretations of the base of the aquifer surface is on the order of ±10–15 m (Fig. S3b). While not a large difference in absolute value, this can represent a significant percentage of aquifer thickness, which is typically <50 m across the region. A similar approach was used to produce interpretive elevation surfaces for the top of the aquifer (bottom of surficial confining layer) and base of the aquifer for the Shellmound high-resolution study area⁵³; these surfaces were incorporated into a model-scale 3D print of the aquifer system (Fig. S4) that is being utilized as an education and outreach tool.

Building on the interpreted base of MRVA surface (Fig. 5a), we developed a continuous classification scheme throughout the entire AEM survey area intended to aid both hydrostratigraphic and geologic interpretation (Figs. 3d–f and 5c, d, “Methods” section). Numeric classes binned from low (1) to high (10) resistivity indicate groups of material expected to have similar lithology and hydraulic conductivity, thereby providing a classification with hydrostratigraphic importance. For example, where coarse-grained MRVA sediments are in direct contact with Tertiary sands of the Claiborne Group aquifers (e.g., MCAQ), resistivity values alone do not distinguish these units. As such, they may behave as a connected hydrologic unit and are therefore considered as part of the same hydrostratigraphic class. To distinguish geologic units, each class is split into two categories (MRVA or non-MRVA) based on whether or not they are within the MAP region and fall above the previously defined base of MRVA surface (Fig. 5a). In cross-section view (Fig. 3d–f) near the typical transition between the MRVA and deeper Tertiary units ~35 m depth (Fig. 5c), both MRVA and non-MRVA categories can be seen for the classes that span this interface where similar resistivity values are found above and below. One notable feature here includes the channel of Quaternary sediment (MRVA class 3–4) running southwest from southern Mississippi, and adjacent to the west margin of the present-day Mississippi River, which captures a subtle channel incised into the top of the more conductive VKBG (Quaternary channel, Fig. 5c). At depth, classifications generally reflect the geometry of mapped Tertiary units, with additional detail on their geometry and heterogeneity.

The permeability of surficial sediments is an important factor in parameterizing recharge for groundwater modeling⁵⁴ and in understanding aquifer vulnerability to contamination from agricultural or industrial operations⁵⁵. Here, the thickness and extent of confining material is defined by measuring the thickness of low-resistivity layers, classes 1–4 (see “Methods” section), within the uppermost 15 m (Fig. 5e). Estimates of recharge that consider soils and surficial geology layers only, without consideration of their thickness or vertical position⁵⁶, may not fully account for recharge potential. By considering the thickness of confining material, regardless of whether it is present at the land surface or buried beneath a surficial layer of more permeable material, AEM results enable a robust metric for recharge potential. For locations where low-resistivity material is absent in the near-surface, increasing potential for aquifer recharge is indicated by the average near-surface classification (Fig. 5e), with increasing classes (i.e., higher resistivity) expected to have a coarser texture and thus more connection potential between recharge and groundwater.

In addition to characterizing aquifer connectivity at the land surface (Fig. 5e), we also consider connectivity between the base of the aquifer and deeper Tertiary formations (Fig. 5f). Connectivity at the base of the aquifer is defined by considering the vertically integrated electrical conductance (vic, see “Methods” section) from 25 m above to 25 m below the base of the MRVA. Areas where the MRVA is underlain by confining Tertiary units, such as MDWY or VKBG, are characterized by low connectivity values, whereas regions of MCAQ subcrop indicate high connectivity (Fig. 5f). Between these end-members, the connectivity metric indicates a range of potential for hydrologic connection between the shallow and deep aquifer systems.

Structure and connectivity beneath rivers

Through numerous glacial cycles and concomitant changes in sea level, river networks have shaped the lower Mississippi River valley^33,38 and are important present-day sources of aquifer recharge. To characterize the detailed subsurface structure beneath rivers, >3000 line-km of AEM data were acquired directly along the Mississippi and Arkansas River paths, as well several smaller tributaries (Fig. 2a). For example, native-resolution, ungridded resistivity profiles with ~25 m spacing along flight paths for several of the tributaries (Fig. S5) illustrate similar detail as the main block flight lines in aquifer structure and geometry of the aquifer base and subcropping unit contact (Fig. S2). Following structure along the river paths enables a detailed view of the discontinuous nature of shallow confining materials beneath river systems, as well as how rivers—which are often important surface water-groundwater conduits—may be connected to the aquifer below. A persistent feature in many river profiles is the intermediate-resistivity layer at ~20–30 m depth that was more discontinuous along the main block flight lines and hypothesized earlier as a fine-grained unit separating younger and older Pleistocene sediments also indicated in Fig. S2.

In addition to streamwise surveys of river courses, gridded resistivity models from the complete regional dataset are intersected with the NHDPlus database of flowlines⁵⁷ to produce cross-sections along any river path (Fig. 6b–f). For example, flowlines intersected by the west–east flight-line block but not flown streamwise can be produced, such as the White River (Fig. 6c). Everywhere the resistivity grids intersected a stream, we use the same connectivity metric discussed previously to quantify the magnitude of connection potential between rivers and the underlying aquifer (see “Methods” section). Here, the vic connectivity metric is calculated on a version of the combined regional dataset made in 2-m depth intervals, integrated from 0 to 10 m beneath the river bottom to characterize the presence or absence of fine-grained material beneath rivers (Fig. 6a). In contrast to the MERAS model calibration for streambed conductance, which incorporated 43 streams with limited data for informing parameter values⁵⁴, the results here provide far greater coverage and finer granularity on the expected spatial patterns of streambed conductance that may be incorporated in groundwater model parameterization. To first order, surficial geology surrounding the streams³² is an important factor in the vic connectivity metric, with fine-grained units appearing less connected than coarse or intermediate units (Fig. S6).

Geological controls on groundwater age

Groundwater age can be used to assess groundwater availability by estimating recharge rates, delineating recharge areas, and characterizing aquifer susceptibility to surface contamination^58,59. Tritium is used to qualitatively differentiate groundwater age into modern (recharged in 1953 or later), premodern (recharged prior to 1953), or mixed (combination of modern and premodern water) groundwater categories⁶⁰. Tritium samples collected from the MAP (n = 582) were categorized into modern, mixed, and premodern groups⁶¹ based on the atmospheric tritium input for the well location, sample date, and the measured tritium concentration⁶⁰. Modern ages are expected for surficial alluvial aquifers in a humid region and 39.1% of MRVA samples fall in this grouping⁶¹; however, 17.9% of samples were premodern⁶¹, indicating significant variability in either recharge rates or sources of groundwater to the MRVA. The heterogeneity in MRVA groundwater age implies local-scale control, which also likely varies among the MAP regions depending on the relative importance of surface water recharge, aerial recharge through surficial confining units, or possible upwelling from deeper Tertiary units. Part of the challenge in understanding the heterogeneity of groundwater age in this system was not having regional-scale interpretations of aquifer architecture and connectivity between units.

To investigate the controls of shallow and deep connectivity on groundwater age, derivative product metrics for surface connectivity (Fig. 5e) and MRVA-Tertiary aquifer connectivity (Fig. 5f) are plotted together with tritium age categories (Fig. 7). Samples with premodern age (Fig. 7a) have characteristically high connectivity between the MRVA and subcropping Tertiary units, along with mixed surficial connectivity, suggesting that upwelling from deeper units may control these older groundwater measurements. Conversely, both mixed and modern tritium categories show a broader range of connectivity to deeper units (Fig. 7b,c), while the modern category has the largest fraction of points (85%) that indicate some degree of surface connectivity, suggesting that connection at the surface is a controlling factor for younger samples. While further data are needed on hydrologic gradients in the system to better predict actual flow paths, these insights suggest that geological structure provides at least partial control on groundwater age and is therefore also likely to be an important driver for vertical groundwater transport. The groundwater age observations—especially the preponderance of premodern groundwater—are difficult to understand without system-scale interpretations of aquifer architecture. The finding that geological structure controls groundwater transport and age is not itself a novel concept; however, we demonstrate that system-scale AEM data can map detailed model structure that facilitates a new understanding of hydrologic processes relevant to many applications and geologic settings.

Hidden faults along the CGL

Resistivity cross-sections west of the northeast-trending portion of Crowleys Ridge image significant up-to-the-east vertical offset—as much as ~50–75 m—that can be tracked over a distance of >100 km along two shallow fault strands (Fig. 8). These two faults closely follow the path of the ~10-km-wide CGL, just west of the RR and NMSZ⁶². The western fault splay shows a clear offset of the contact between an Upper Cretaceous high-resistivity layer at depth (McNairy Nacatosh aquifer) and the low-resistivity Tertiary MDWY within the uppermost 100–150 m (Fig. 8b–d), clearly extending to at least the base of the Quaternary aquifer. This western fault is evident on cross-sections between Fig. 8b and d (solid brown line Fig. 8a), with vertical offsets largest in the middle of this segment (up to 75 m) and maximum width ~500–1000 m. The eastern of the two faults is characterized in cross-section view by a dipping conductive layer that appears to terminate at the base of the MDWY, separating the uplifted resistive Cretaceous block to the west from the conductive MDWY to the east. Whether the dipping conductive layer is related to the MDWY or the fault structure itself is unclear. Unlike the western fault, offset is not seen at the top of the MDWY, possibly because of Quaternary erosion of this surface. The eastern of the two structures can be traced on resistivity cross-sections over a greater distance, following the CGL to the southwest before turning south towards the Western Margin fault (Fig. 8).

While generally coincident with the CGL, both hidden faults captured in the AEM data are previously unmapped, and thus provide new insight into the tectonic history of this area. The observed uplift on these faults along the CGL can plausibly help explain the tectonic origin of Crowleys Ridge given the proximity of these features to one another. Previous geophysical and geomorphological observations along and near the ridge margins provide important insight^63,64—especially in relation to the eastern fault indicated here—but have not revealed the source of uplift occurring 10–20 km farther west along the CGL. Existing geophysical data along the ridge margins⁶⁴ and the CGL⁶⁵ indicate faulting and support Quaternary fault activity, but do not have sufficient spatial scale and or resolution to delineate the discrete offset imaged by AEM data over the fault length of >100 km. The 50–75 m uplift (Fig. 8b–d) requires movement since the early Tertiary but is consistent with the Quaternary motion. The offset-wedge geometry suggests that the MDWY was lithified at the time of deformation. Given the lack of surface expression, the majority of offset may have occurred before the 32–45 ka Melville Ridge braid belt found immediately west of Crowleys Ridge³³; however, further age constraint of Quaternary activity is not available.

Dense borehole observations that quantify the overall thickness of Quaternary sediments identify the Pliocene–Pleistocene unconformity⁴¹ and support an argument for 53 m of Quaternary uplift³⁴ in the region (Fig. 8a). However, the typical minimum separation between boreholes is >2 km and borehole density declines west of the CGL, making them also insufficient to identify this narrow fault feature as a potential source of uplift. Although the CGL and identified faults are removed from the seismically active region of the NMSZ, the possibility of Quaternary tectonic activity and offset on the faults identified here is important to understand past variations in seismicity and future hazard in the region⁶².

Multi-scale mapping

The multiple phases of AEM mapping (Fig. 2a), along with targeted ground-based⁶⁶ and waterborne⁶⁷ geophysical data collection, provide an excellent case study in the value of subsurface mapping over multiple scales. From regional surveys aggregated to 3 km line spacing with 1 km grid cells to high-resolution AEM data in the Shellmound area with 0.25–1 km line spacing with 100 m grid cells, we are able to illustrate the value of an order of magnitude increase in flight-line spacing.

The horizontal and vertical resolution of airborne geophysical surveys defines the scale of investigation, which depends on survey design parameters including instrument type, flight-line spacing, and total kilometers flown. Spatially extensive data capturing individual features several kilometers in size commensurate with mapping scales of ~1:1,000,000 (Figs. 2, 5, and  9a) are mapped with flight-line spacing on the order of 1–10 km and effectively inform regional-scale hydrologic models and decision support. High-resolution airborne data with sub-kilometer flight-line spacing greatly improve the resolution of smaller features appropriate for local ~1:100,000-scale mapping (Fig. 9c and Fig. S7).

The high-resolution data near Shellmound, Mississippi, map coarse-grained (high-resistivity) sediments associated with river meanders that also tend to have relatively high potassium (K) abundance (Fig. S7a, b, east), in contrast with backswamp deposits and other fine-grained overbank sediments that have relatively low-resistivity and greater thorium (Th) abundance. At 50–55 m depth (Fig. 9c and Fig. S7d), resistivity models map a buried paleochannel in the southeast corner of the grid that appears incised to depths of at least 75 m. This channel was previously unknown and may be a relict of the ancestral Mississippi–Ohio river system that flowed east of Crowleys Ridge during the Pleistocene³³, creating conduits for local groundwater flow. Deep channels such as this are likely present throughout the region, but are under-sampled by borehole observations and even the regional geophysical flight lines (Figs. 2d–h, 3a–c).

The combination of spatially extensive and high-resolution data enables detection and mapping of small-scale or hidden features that would otherwise be missed (i.e., the known unknowns in the subsurface). Our results demonstrate the value in AEM for imaging localized pathways for groundwater connection above and below the aquifer (Fig. 5e, f), buried paleochannels (Fig. 9c), and hidden faults with narrow zones of uplift previously unmapped (Fig. 8). These discrete and relatively small features may have an outsized and nonlinear impact on model predictions, hazard assessments, and decision making, but are impossible to recognize with sparse or limited spatial coverage.

Outlook and future directions

An emerging example of “big data” in geoscience, AEM is expanding the breadth of information that can be applied to near-surface investigation and modeling where detail about subsurface complexity is lacking. Many process-based hydrologic models have evolved to function in the absence of detailed subsurface structural information, largely because these data are typically not available and because incorrect assumptions about subsurface structure and unwarranted complexity lead to modeling errors⁶⁸. Technological advances in remote sensing and airborne geophysics present new opportunities to advance the level of geological detail that can be incorporated in hydrologic models through the entire vertical extent of an aquifer system. As “big data” have become available to constrain details of model structure, commensurate advances in modeling, such as machine learning and uncertainty quantification⁶⁹, will also be needed to realize the full potential of airborne geophysical datasets.

As the dimension and complexity of questions that may be asked of a model or array of models (i.e., the “decision-space”) continue to increase, the need for detailed data to support these questions also increases. When allocating resources for data collection, a balance should be sought between the optimality of data collection to address a narrow decision-space (such as a single, focused question) and robustness of the investment to address future questions that expand the decision-space. The variety of applications demonstrated in this study highlights the robustness of large-scale AEM data in this respect. Although AEM surveys can involve high absolute cost, at large scales they may be 3–4 orders of magnitude less expensive on a per-data-point or per-square-kilometer basis compared with traditional ground-based surveys or drilling. Because large surveys can cover parts of multiple counties and states and support the interests of multiple scientific disciplines or stakeholder interests, a community-driven approach to the acquisition of these foundational geoscientific datasets is advantageous. Benefits of a coordinated approach to acquiring system-scale AEM data include reduced costs to individual participants, leveraging resources to acquire more data than any individual group could achieve, and ensuring data consistency across the study area that maximizes their value.

Much as lidar has transformed our understanding of Earth’s surface, airborne geophysical data extend our view into the subsurface, transforming our ability to inform three-dimensional mapping from catchment to basin scales (Fig. 1b). Here, we demonstrated that system-scale AEM data provide a robust platform from which to address a host of subsurface questions. Airborne geophysical datasets represent the next generation in subsurface mapping, capable of filling in the gaps between existing boreholes with an order of magnitude or greater increase in data density. This will provide nationally consistent near-surface geologic datasets as a foundation for three-dimensional geologic interpretations, hydrologic models, and other subsurface studies that rely on detailed measurements of difficult to access belowground properties within 200–300 m of Earth’s surface (i.e., lidar for the subsurface).

Source: Resources - nature.com