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

The Mississippi Alluvial Plain hosts one of the most prolific shallow aquifer systems in the United States but is experiencing chronic groundwater decline. The Reelfoot rift and New Madrid seismic zone underlie the region and represent an important and poorly understood seismic hazard. Despite its societal and economic importance, the shallow subsurface architecture has not been mapped with the spatial resolution needed for effective management. Here, we present airborne electromagnetic, magnetic, and radiometric observations, measured over more than 43,000 flight-line-kilometers, which collectively provide a system-scale snapshot of the entire region. We develop detailed maps of aquifer connectivity and shallow geologic structure, infer relationships between structure and groundwater age, and identify previously unseen paleochannels and shallow fault structures. This dataset demonstrates how regional-scale airborne geophysics can close a scale gap in Earth observation by providing observational data at suitable scales and resolutions to improve our understanding of subsurface structures. Integration of large-scale airborne geophysical surveys can enable the system-scale assessment of aquifer characteristics and improved mapping of shallow subsurface structures, according to a comprehensive investigation of the Mississippi Alluvial Plain

G roundwater is the world's most extracted raw material, with annual withdrawals of 600-1100 km 3 /year 1,2 . Withdrawals from freshwater aquifers supply nearly half of the world's drinking water and 43% of water for irrigation 1,3 . Increasing withdrawals from finite groundwater stores driven by changes in the availability of clean surface water, climate, population growth, and management have led to significant stress on global groundwater systems [4][5][6] . Despite the critical importance of groundwater, aquifers remain a largely hidden resource that have not received the level of detailed exploration typically given to other natural resources such as oil and gas reservoirs 5 .
Groundwater withdrawals for irrigation made up~20% of the total freshwater use for 2015 in the United States, of which 58% was extracted from the High Plains (HP, 20.3%), Mississippi Alluvial Plain (MAP, 20.5%), and Central Valley (CV, 17.7%) aquifer systems ( Fig. 1a) 7,8 . Sustained withdrawals in the MAP region have led to significant groundwater-level declines in excess of 100 feet (30 m) from pre-development conditions before significant groundwater pumping in the mid-twentieth century through 2007, primarily in Arkansas and Mississippi 9 . The Mississippi River Valley alluvial aquifer (MRVA), the surficial aquifer within the MAP region, was among the most heavily withdrawn aquifers for irrigation in 2015 (11.7 Bg/day) 7 (Fig. 1a). Approximately 8 million acres (32% of the MRVA footprint) are irrigated 10 , with the overwhelming majority derived from the shallow, highly productive MRVA. The total revenue and economic impact of row crop agriculture were nearly $12 billion in 2017 11 .
Groundwater management over large regions often relies on scant data about the structure, extent, and connectivity of principal freshwater aquifers. The most detailed information about aquifer structure often comes from an ad hoc collection of borehole data acquired for disparate purposes over years to decades from various sources such as municipal, irrigation, or resource exploration wells, and is curated by multiple agencies, with no official central repository or common data standards. Regional-scale borehole information is highly variable in its spatial coverage, content, and data quality. Borehole data do not uniformly sample the subsurface either in land area or in depth, and information about data uncertainty is often missing or anecdotal. The need for consistent and spatially extensive data has led to rapid advancement in the use of remote sensing techniques. For example, advances in satellite and high-altitude airborne remote sensing methods have dramatically advanced our ability to understand water use and water budgets at regional to global scales [12][13][14][15] . Together with ground-based observations, remotely sensed data of changes at the Earth's surface are greatly improving the predictive capabilities of groundwater simulations 16 .
Airborne electromagnetic (AEM) methods are emerging as a highly cost-effective means for extending remote sensing capabilities to the characterization of the subsurface at the watershed to basin scale, filling a critical scale gap in Earth observation between ground-based and satellite or airborne remote sensing methods (Fig. 1b). AEM has been applied to groundwater studies for decades [17][18][19][20][21][22][23][24] , with numerous applications in studies of the cryosphere [25][26][27] , infrastructure 28 , and hazards 29,30 . However, most AEM studies have been applied to relatively limited areas (a few hundred to a few thousand square kilometers) that represent only part of a larger geological or hydrological system. In this study, we present the first AEM effort to map an entire major aquifer system in the United States, covering >140,000 km 2 within the MAP region.
The modern lower Mississippi River valley has primarily been shaped by the combined forces of glaciation and sea-level change [31][32][33] . Between 2.4 Ma and 250 ka, erosion removed a thick sequence of Pliocene gravels along with tens of meters of underlying earlier Tertiary sediments 34 , leaving <100 m of discontinuous Pleistocene sediments unconformably overlying a thick sequence of Cretaceous-Tertiary sediments of the Mississippi Embayment 35,36 . These Pliocene gravels are preserved only at higher elevations in a small area on Crowleys Ridge and outside the footprint of the MRVA 34,37 . Borehole data provide constraints on the thickness of Quaternary MRVA sediments 38,39 , and optically stimulated luminescence (OSL) measurements helped to refine the chronology of the valley over the last glacial cycle during the past~100 ka 33,40 . In the northern part of the MAP study area, tectonic activity associated with faults of the Cambrian Reelfoot rift (RR) has caused deformation of the Pliocene-Pleistocene unconformity that forms the base of the alluvial aquifer system. Fault activity and uplift have continued in the Quaternary with implications for modern-day hydrologic systems 34,41 as well as seismicity related to the New Madrid seismic zone (NMSZ) 42 .
The US Geological Survey's MAP Regional Water Availability Study began in 2016 to pilot a next-generation approach to integrated water availability assessments. Original development of the hydrogeologic framework for the earlier Mississippi Embayment Regional Aquifer Study (MERAS) model 43 incorporated >2600 geophysical logs and well data within the Mississippi Embayment 44 focused on the deeper Tertiary units that are primarily used for drinking water. Borehole constraints on the surficial MRVA have been sparse by comparison because of its primary use as an irrigation source. Out of the available 2652 boreholes, only 381 were able to define the base of the MRVA, or roughly one control point for every 200 km 2 . As part of the new MAP assessment, a regional AEM survey was designed as the foundation for updating the hydrogeologic framework of the MRVA together with updated borehole compilations 39,41 .
In this first-of-its-kind regional AEM study, we demonstrate the utility of airborne geophysical data for improving understanding of the structure of the surficial aquifer and surrounding geological formations that are of critical importance not only for effective water resource management in the region but also have broader implications for geologic mapping and hazards investigations. Our data refine system-scale detail of the MRVA, surficial confining material, and the configuration of deeper Tertiary units to improve understanding of regional hydrology. High-resolution flight-line data and grids reveal previously unseen features such as buried paleochannels in Mississippi and faults adjacent to the NMSZ in Arkansas and Missouri-structures that have potential for significant impact on the Earth system compared with their relatively small footprint (which makes detection difficult) within the larger study area. Beyond the results presented here, we expect these system-scale data and models will have broad community value for comprehensive and diverse investigations of the MAP region going forward, and demonstrate the value of taking a regional approach to understanding the subsurface that is foundational to many scientific and societal studies.

Results and discussion
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 2 (Fig. 2a, "Methods" section). Data collection included a high-resolution survey over~1000 km 2 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 Fig. 2 Airborne geophysical survey coverage and summary of regional datasets. a Primary management regions in the MAP study area, with flight lines for each of the three phases of data collection completed through early 2020 (CR Crowleys Ridge). Results from the combined regional surveys gridded onto the 1 km National Hydrogeologic Grid 45 include b radiometric data presented as a ternary diagram that indicates the relative abundance of K, U, and Th in surficial sediments with areas of Holocene (H)-and Pleistocene (P)-aged sediments indicated 32 ; c the residual magnetic intensity map (in nanoTeslas, nT) shows faults 46  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 32 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. 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 43 . 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 33 , 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 33 , there is little constraint on the age of deeper Quaternary sediments beyond the estimate of their deposition after 250 ka 34 . Given the presence of shallow braid belt deposits older than the last sea-level lowstand~20 ka 33 , 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 10 SC 2.79 = 617 μS/cm) and deeper Tertiary units (median log 10 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 48 . Correlation between SC-and AEM-derived resistivity (Fig. 4) by Fig. 4 Relationships between groundwater salinity and AEM-derived resistivity. a-f Measured groundwater specific conductance (SC) versus AEMderived resistivity at borehole locations, organized by region (Fig. 2). SC measurements from MRVA, terrace (TRRC), and loess layers are shown in blue, with measurements from deeper Tertiary units in yellow. Black dashed curves illustrate theoretical SC-resistivity relationships for varying amounts of surface conductivity (σ s ) caused by the increasing fraction of fine-grained sediment. g Map view of a high-SC cluster in White County, Arkansas, correlates with decreased shallow resistivity in cross-section view (h) and corresponds to the high-SC observations circled in (a).
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 39 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 51 (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 52 . 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 38 , 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 53 ; 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 54 and in understanding aquifer vulnerability to contamination from agricultural or industrial operations 55 . Here, the thickness and extent of confining material is defined by measuring the thickness of lowresistivity 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 56 , 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 nearsurface 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 watergroundwater conduits-may be connected to the aquifer below. A persistent feature in many river profiles is the intermediateresistivity 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 57 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 54 , 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 32 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 60 . Tritium samples collected from the MAP (n = 582) were categorized into modern, mixed, and premodern groups 61 based on the atmospheric tritium input for the well location, sample date, and the measured tritium concentration 60 . Modern ages are expected for surficial alluvial aquifers in a humid region and 39.1% of MRVA samples fall in this grouping 61 ; however, 17.9% of samples were premodern 61 , 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 mthat 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 62 . 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 64 and the CGL 65 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 Fig. 7 Relationship between groundwater age and AEM-derived aquifer structure. Sample depth for tritium dates classified as premodern (a), mixed (b), or modern (c) is plotted against the degree of connectivity between the MRVA and deeper Tertiary units (Fig. 5f), with point colors that represent the degree of shallow aquifer connectivity (Fig. 5e).
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 33 ; 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 41 and support an argument for 53 m of Quaternary uplift 34 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 62 .
Multi-scale mapping. The multiple phases of AEM mapping (Fig. 2a), along with targeted ground-based 66 and waterborne 67 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 flightline 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 flightline 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 33 , 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 68 . 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 69 , 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-datapoint 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 nearsurface 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).

Methods
Airborne data collection. Airborne geophysical data were acquired in three phases through a competitive award to CGG Airborne, using two different AEM platforms that also included magnetic and radiometric sensors: (Phase 1) A high-resolution grid with flight-line spacing between 0.25 and 1 km covered an area of~1000 km 2 just northwest of Shellmound, Mississippi. (Phase 2) A regional grid covered over 95,000 km 2 with 6-12 km between flight lines over most of the area, plus an area of higher density (1.5 km line spacing) over part of southern Crowleys Ridge in Arkansas. This regional survey also included~2000 line-km of data acquired directly along several streams and rivers to characterize their potential for enhancing connectivity between surface water and groundwater. (Phase 3) A second regional survey extended coverage to an area of 140,000 km 2 with 6 km between flight lines, interwoven with the first regional grid such that the combined line spacing was 3 km over most of the study area. An additional 1400 line-km were acquired streamwise along the entire Mississippi River within the study area, as well as the Arkansas River from Little Fig. 9 Multi-scale mapping with airborne and ground-based methods. a Regional-scale structure mapped on 1 km grid cells interpolated from~3 to 6 kmspaced flight-line data over the entire MAP region. b Regional-scale structure enlarged to the~30 × 30 km Shellmound study area, compared with highresolution structure mapped on 100 m grid cells interpolated from 0.25 to 1 km-spaced flight-line data from the Shellmound study area (c) shows the ability to resolve detailed buried channel structure. d Comparison of near-surface high-resolution Shellmound AEM data (background) with very high-resolution ground-based electromagnetic data 66 acquired with a sensor towed over~4 km 2 on survey lines spaced by 25 m.
Rock, Arkansas, to its confluence with the Mississippi River. The first two phases of data collection were acquired with the helicopter-based Resolve frequencydomain AEM system, while the third phase was acquired with the fixed-wing time-domain Tempest AEM system. The Resolve AEM system was flown in the first two phases (Fig. 2a), first in March 2018 for a high-resolution study near Shellmound, Mississippi 53,70 , and again in the regional mapping effort from November 2018 to February 2019 that covered most of the MAP study area 71 . The Resolve AEM system is a frequencydomain electromagnetic sensor that operates with six independent transmitter-receiver coil pairs at discrete frequencies from 400 to 140,000 Hz 18 . Five of the coil pairs are oriented in a horizontal coplanar configuration (horizontal coils with vertical magnetic dipole) and separated by 7.9 m, with one coil pair in a vertical coaxial configuration (vertical coils with horizontal magnetic dipole in the flight direction) separated by 9.1 m. Data are recorded at a 10-Hz sampling interval with nominal flight speed 130 km/h and instrument height 30 m above the terrain. The Resolve instrument is slung~30 m beneath a Eurocopter AS350 B2 helicopter. The typical depth of investigation for the Resolve system is 60-100 m, with high vertical resolution in the uppermost 2-5 m below land surface. The Tempest AEM system was flown in the third phase of surveys, covering the entire study area (Fig. 2a) from November 2019 to March 2020 72 . Tempest is a time-domain electromagnetic sensor 73 acquired from a fixed-wing platform (Cessna 208B), with a transmitter loop attached to the airframe and a threecomponent receiver coil in a bird that is towed behind the aircraft. The Tempest transmitter is a horizontal loop with an area 155 m 2 and an average moment 43,400 Am 2 . The three-component receiver is towed nominally 108 m behind and 50 m below the transmitter, with actual position during flight recorded by GPS. Inertial measurement unit sensors are mounted in the aircraft and the receiver bird measured transmitter and receiver orientation, respectively. Recorded dB/dt data for each channel are transformed to B-field responses that would be observed for a perfect 100% duty cycle transmitter waveform 73 , and are provided in 15 channels with window center times from 11 to 13,400 µs. Data are recorded at a 5 Hz sampling interval with a nominal flight speed of 200 km/h and a flying height of 120 m above the terrain. The typical depth of investigation for the Tempest system is 200-300 m, with moderate vertical resolution in the uppermost 5-10 m below land surface. Complete details of system and acquisition parameters for each survey phase can be found in the contractor-provided reports [70][71][72] .
Inversion of AEM data and combined geophysical products. Resolve AEM data were inverted to recover subsurface electrical resistivity structure using the Aarhus Workbench software (Aarhus GeoSoftware, Aarhus, Denmark) that builds on the AarhusInv code 74 . Input data included 12 data channels of in-phase and quadrature measurements from each of the six coil pairs, along with laser altimeter measurements of the instrument height above ground. Data were manually edited to remove noisy measurements caused by manmade infrastructure before inversion. Inverted resistivity models were estimated on a 30-layer structure from depths of 1-125 m below ground, and models were output approximately every 25 m along flight lines. A laterally constrained inversion process added lateral constraints between neighboring models to promote spatial continuity across models. In addition to resistivity values, the algorithm also produced estimates of transmitter height above ground, parameter uncertainty, and depth of investigation 75 . Data and inverted model results are available in the online ScienceBase repository 70,71 .
Tempest AEM data were inverted to recover subsurface electrical resistivity structure using the GALEI code developed by Geoscience Australia, including specific functionality for the Tempest system 76 . Input data included 30 data channels of delivered B-field observations for both x-and z-component receiver coils (15 channels each). System geometry included transmitter altitude and attitude, as well as transmitter-receiver offset and receiver attitude, all of which were input into the inversion. Inverted resistivity models were estimated on a 30layer structure from depths of 4-400 m below ground, and models were output approximately every 75 m along flight lines. In addition to resistivity values, the algorithm also produced estimates of transmitter height above ground, parameter uncertainty, and depth of investigation 75 . Data and inverted model results are available in the online ScienceBase repository 72 .
We combined the Resolve and Tempest resistivity models in three steps. First, for each dataset, resistivity models were discretized onto a common 1 km resolution grid 45 for 5 m depth intervals, including model values down to the estimated depth of investigation for each grid point. Second, we kriged each dataset using the geostatspy library (https://github.com/GeostatsGuy/GeostatsPy) that incorporates GSLIB functions 77 , using empirical semivariograms to define the correlation structure for each layer. Lastly, the Resolve and Tempest kriged grids for each depth interval were combined using a weighted function that smoothly transitioned from a mix of 0.52/0.48 Resolve/Tempest in the shallowest layer to 0.00/1.00 Resolve/Tempest by a depth of 75 m 78 .
Derivative products Base aquifer. The elevation of the base of the MRVA was determined using a supervised machine learning algorithm 51 that uses 6130 existing borehole picks 39 along with 364 additional manual picks and AEM-derived resistivity data, as training. Based on these inputs, the algorithm predicts the base of MRVA elevation at all AEM data locations, with values interpolated to a regular grid with 1 km spatial resolution 79 .
Facies classification. The kriged resistivity grids described above were categorized into ten classes (Table 1) intended to group materials expected to have similar hydrologic and geologic properties based on their resistivity. In some cases, different juxtaposed geologic materials may have similar lithology, such as MRVA sands that overlie coarse-grained Tertiary sands, which cannot be distinguished based on their resistivity alone, and may, in fact, behave as a connected hydrostratigraphic unit. To provide additional geological interpretation, we also classified resistivity grids based on their position relative to the base of MRVA surface discussed above. For example, similar resistivity values above and below the base of MRVA surface would belong to the same class number, but are separated into MRVA and non-MRVA groups based on their vertical position. In this fashion, our classification system includes both hydrostratigraphic and geologic information. The non-MRVA class captures both Tertiary units beneath the MRVA and rocks and sediments outside the footprint of the MAP region, including Crowleys Ridge. Resistivity thresholds that separate the classes are chosen on a uniform logarithmic spacing between reasonable minimum and maximum values observed in the dataset. Classes qualitatively correspond with known lithologic transitions discussed above; however, sufficient data are not available to produce a quantitative calibration of class groupings against borehole and other geologic or hydrologic datasets for the entire region (Table 1).
Connectivity. We defined a vertically integrated connectivity metric as the vic across an interval of interest: For each layer (i) within an interval, the electrical conductance is defined as the product of layer thickness (t i ) and the inverse of the layer resistivity (ρ i ). The vic is the sum of layer conductances and is intended to highlight the presence and thickness of low-resistivity material that may be a barrier to vertical hydrologic connectivity across the interval.