Multiple-point statistical modeling of three-dimensional glacial aquifer heterogeneity for improved groundwater management

Quaternary glacial aquifers are important water sources for irrigation in many agricultural regions, including eastern Nebraska, USA. Quaternary glacial aquifers are heterogeneous, with juxtaposed low-permeability and high-permeability hydrofacies. Managing groundwater in such aquifers requires a realistic groundwater-flow model parameterization, and characterization of the aquifer geometry, spatial distribution of aquifer properties, and local aquifer interconnectedness. Despite its importance in considering uncertainty during decision-making, hydrofacies probabilities generated from multiple-point statistics (MPS) are not widely applied for groundwater model parameterization and groundwater management zone delineation. This study used a combination of soft data, a cognitive training image, and hard data to generate 100 three-dimensional (3D) conditional aquifer heterogeneity realizations. The most probable model (probability of hydrofacies) was then computed at node spacing of 200 × 200 × 3 m and validated using groundwater-level hydrographs. The resulting hydrofacies probability grids revealed variations in aquifer geometry, locally disconnected aquifer systems, recharge pathways, and hydrologic barriers. The profiles from hydrofacies probability at various locations show spatial variability of the streambed and aquifer connectivity. Groundwater-level hydrographs show evidence of these aquifer characteristics, verifying the general structure of the model. Using the MPS-generated 3D hydrofacies probability and hydrologic data, a novel workflow was developed in order to better define high-resolution groundwater management zones and strategies. In general, the conditional probability of hydrofacies helps improve the understanding of glacial aquifer heterogeneity, the characterization of aquifer-to-aquifer and streambed-aquifer connections, and the delineation of groundwater management zones. This MPS workflow can be adapted to other areas for modeling 3D aquifer heterogeneity using multisource data.


Introduction
Glacial aquifer systems supply water for various uses worldwide, particularly in the northern latitudes that experienced widespread glaciation during the Quaternary (Miller 1999;Russell et al. 2004;Ehlers and Gibbard 2007;Margat and van der Gun 2013).In the United States, glacial aquifers underlie parts of 26 states and supply water to about 42 million people, or 5% of the US population (Yager et al. 2018).Aquifer heterogeneity is defined as the spatial variation of aquifer hydraulic properties (Hiscock and Bense 2021).Glacial aquifer systems are heterogeneous and comprise sediment assemblages with a wide range of hydraulic properties (Brodzikowski and van Loon 1987;Anderson 1989;Whittaker and Teutsch 1999;Benn and Evans 2014).Aquifer heterogeneity results in complex patterns of groundwater flow, which impact head distribution and travel time (He et al. 2013), groundwater and surface interactions (Fleckenstein et al. 2006;Zhao and Illman 2017), limit infiltration rate and groundwater discharge (Åberg et al. 2021) and managed aquifer recharge (Maples et al. 2020).Aquifer heterogeneity produces tortuous groundwater flow paths (Comte et al. 2019) and challenges groundwater management (Kelly et al. 1 3 2013).Characterizing the spatial scales of heterogeneity in the subsurface volume of interest is essential to appropriately simulate flow and manage groundwater resources (Giudici 2010).The relevant data to characterize aquifer heterogeneity at a scale of interest are rarely available, and as a result, aquifer heterogeneity is simplified for groundwater simulation.Simplifying heterogeneity may lead to unrealistic aquifer characteristics (Lukjan and Chalermyanont 2017) and incorrect aquifer connections, resulting in overgeneralized groundwater flow paths.
Sustainable groundwater management requires an improved understanding of aquifer-to-aquifer connections, stream-aquifer connections, and vertical and horizontal water flow components and processes between or within an aquifer.Groundwater management activities such as regulating well pumping and well interferences, siting new pumping wells, protecting streamflow depletion, safeguarding recharge zones and enhancing aquifer recharge, necessitate a realistic aquifer framework and properties.Groundwater management decisions often rely on the results of groundwater models and monitoring systems.Well hydrographs are often used to map aquifer heterogeneity and properties during calibrations; however, even with good calibration results, uncertainty may persist.A unique set of parameters cannot be obtained during calibration-multiple hydraulic conductivity values may provide a good fit, and solutions are nonunique (Konikow and Bredehoeft 1992).Inadequate conceptual geological model and aquifer parameter estimations are sources of uncertainty in groundwater flow models (Rojas et al. 2010;Refsgaard et al. 2012;Enemark et al. 2019;Yin et al. 2021).Errors in aquifer heterogeneity and structure can result in biases related to parameterization, calibration and prediction (Christensen et al. 2017b).
Comparing the estimated aquifer parameters by the inverse calibration model to prior knowledge of aquifer structures or a 3D aquifer realization model can help validate the inverse modelbased estimated aquifer properties and model accuracy, reducing the potential impact of ill-posed calibrations on model-aided decision-making.Three-dimensional (3D) modeling of hydrofacies is fundamental to constructing a groundwater flow model for groundwater management.Hydrofacies probability is also a useful tool for understanding and assessing uncertainty in an aquifer structure and aquifer properties.However, it is not widely used in groundwater model parameterization and groundwater management zone delineation due to a lack of data and difficulties in implementing the concept in practical applications.
Aquifer permeability is influenced by lithofacies distribution (Bersezio et al. 1999;Zappa et al. 2006;Suriamin and Pranter 2018).Lithofacies from borehole data are frequently used to model hydrofacies.Different lithofacies can be grouped into single hydrofacies (Klingbeil et al. 1999;De Caro et al. 2020).Estimation of hydraulic conductivity is typically done by analyzing lithofacies data and limited pumping test wells to model groundwater flow.The groundwater flow model necessitates realistic parameterization of spatial variations in aquifer properties.Accurately estimating spatial variations of complex glacial aquifer hydraulic conductivities requires a large amount of data.The estimation of realistic aquifer heterogeneity and properties can be aided by modeling spatial variations in hydrofacies.The hydrofacies model can be used to estimate the hydraulic conductivity of an aquifer (Theel et al. 2020;Xue et al. 2022).Geophysical data are often used to enhance the horizontal resolution of hydrofacies modeling, but translating geophysical data to hydrofacies can be challenging due to the nonuniqueness of the information derived from geophysical data.Airborne electromagnetics (AEM) is a geophysical technique that provides spatially dense data between boreholes and can yield critical information on the aquifer extent and confining units (Hanson et al. 2012;Korus et al. 2013Korus et al. , 2017Korus et al. , 2021;;Korus 2018).To interpret AEM data, the AEM data must first be inverted and interpolated to generate a 3D resistivity grid.Then the resistivity-borehole lithological relationship is used to convert 3D AEM resistivity data into hydrofacies.Combining AEM-derived resistivity-depth models with borehole logs can improve geological heterogeneity simulation and reliability (Jørgensen et al. 2013;Foged et al. 2014;Korus et al. 2017).The integration of AEM with a borehole can be utilized to map the aquifer system (Knight et al. 2018), estimate the upper top of the saturated zone (Dewar and Knight 2020) and locate incised valley fill for recharging the aquifer (Knight et al. 2022).Combining AEM data with well hydrograph data is useful for mapping impermeable aquifer boundaries (Korus 2018).The information derived from AEM data is also valuable in bolstering the accuracy of groundwater models.Marker et al. (2015) employed AEM data in combination with borehole data to produce clay fraction estimates, categorized them into hydrostratigraphic zones, estimated hydraulic conductivity for each zone, used these values in groundwater calibration, and observed an enhancement in hydrological performance.Christensen et al. (2017a) used AEM data to create 3D distributions of hydraulic conductivity and reported that the accuracy of the groundwater predictive model improved, specifically, predictions of recharge area, head change and stream discharge.
Three-dimensional modeling of glacial aquifer heterogeneity is challenging because it requires spatially dense data and sophisticated methods.The cognitive approach produces only a single aquifer heterogeneity and involves significant effort and expertise to create multiple 3D realizations of aquifer heterogeneity in a glacial region.Geostatistical simulations such as two-point and multiple-point statistics (MPS) are also commonly used to model aquifer heterogeneity.The MPS approach has many advantages over cognitive and two-point geostatistical methods (Liu et al. 2005;Hu and Chugunova 2008).Unlike two-point geostatistics, which use a variogram approach to generate geological heterogeneity, MPS uses statistics from a training image (TI) and better simulates complex geological patterns than two-point techniques (Strebelle 2002).MPS generates multiple plausible realizations by combining TI, hard data (certain borehole data) and soft data (geophysical data); however, locational measurement and driller's borehole lithological logs descriptions are always subject to uncertainty and borehole data are not 100% certain.Despite such errors, hard data are treated as certain and used to constrain the MPS realization.MPS realizations are compelled to perfectly match the hard data during simulations, but MPS realizations are not forced to match soft data (Hansen et al. 2018;Vilhelmsen et al. 2019;Madsen et al. 2021).
In areas where there are no hard data, MPS fills simulation grids by borrowing geological features from a TI and soft data.Statistical information is derived from TI using a search tree template and stored in a tree (Liu 2006).The selection of an appropriate search template is vital to capture geological heterogeneity and information presented in the TI.A small search mask may exclude vital geological features in the TI, and too large a search mask prolongs MPS simulation time and leads to nonrepresentative statistics (Hu and Chugunova 2008).The TI pattern that appears in the search mask is recorded in the search tree (Ba et al. 2019).
Near the eastern margin of the High Plains Aquifer (HPA) in Nebraska, the geological system consists of complex glacial deposits and productive paleovalley (buried valley) aquifers (Divine et al. 2009;Korus et al. 2017).The boundary of the HPA in eastern Nebraska is concealed by glacial deposits and loess, and there are multiple generations of paleovalleys in the subsurface (Korus et al. 2017).These paleovalley aquifers are complex and heterogeneous, resulting in significant spatial variability in aquifer properties.The variability of aquifer properties, as well as the presence of low permeability hydrofacies, can limit water flow between and within the aquifers as well as recharge to the aquifer.Effective groundwater management in this area requires an understanding of paleovalley aquifer continuity, confining units, and features controlling groundwater flow.This study aims to answer the following questions by simulating detailed glacial aquifer heterogeneity: (1) Are the palaeovalley aquifers connected?(2) Is the streambed-aquifer system connected?(3) Can a 3D MPS-generated probability be used to define management boundaries?
Probabilistic simulations of glacial aquifer heterogeneity using MPS are not often utilized to analyze the connection between buried valleys aquifers, aquifers and streambeds and to define groundwater management zones.Simulations of various conditional realizations aid in determining plausible glacial aquifer heterogeneity and aquifer-to-aquifer connectivity, which is an important step in considering the uncertainty of aquifer characterization in groundwater modeling and management.Visual inspections and borehole logs are commonly utilized to validate hydrostratigraphic units produced by MPS (Barfod et al. 2018a).Although well hydrographs can provide information about the properties of aquifers and hydrostratigraphic units, they are rarely used to validate MPS-generated hydrofacies and aquifer heterogeneity.Unlike other examples, this study uses a novel application of groundwater head monitoring for qualitative validation of aquifer heterogeneity and computed 3D hydrofacies probability to account for uncertainty.To the authors' knowledge, hydrographs have not been applied to validate glacial aquifer heterogeneity.This study demonstrated that MPS can be a useful tool for modeling glacial-aquifer heterogeneity.Furthermore, a workflow for delineating management boundaries that take into account the spatially variable nature of these aquifers was provided, which is critical for tailored management solutions.

Physical setting
The study area is located in eastern Nebraska, USA (Fig. 1).It includes the watersheds of Elm Creek, Taylor Creek, and Loseke Creek, which flow into Shell Creek, a tributary of the Platte River.The elevation of the study area ranges between 424 and 533 m.Groundwater is managed by the Lower Platte North Natural Resources District, which has established management rules for this area because of periodic groundwater shortages related to well interference and large-magnitude groundwater-level declines during irrigation pumping.The management area is approximately 643 km 2 and is about 40 km in the E-W direction and 17 km in the N-S direction.
Shell Creek watershed lies near the maximum western extent of Pleistocene glaciation in North America.The Laurentide Ice Sheet advanced over preglacial and proglacial deposits in this area, producing a complex arrangement of buried valleys, tills, and spatially variable sand and gravel bodies.Glacial sediments are mantled by multiple loess units.Bedrock consists of impermeable Cretaceous rocks, principally shale, chalky shale, and limestone, which serve as the lower boundary for groundwater flow.Unconsolidated sand and gravel units above bedrock serve as the primary aquifers for municipal and irrigation purposes.The glacial aquifers system consists of sand and gravel deposits in buried valleys atop bedrock (paleovalleys), and discontinuous sand and gravel bodies within unconsolidated glacial sediments (Divine et al. 2009).

Borehole data processing
Boreholes logs were obtained from the Nebraska Department of Natural Resources (NDNR) and test holes were drilled by the Nebraska Conservation and Survey Division (CSD) (Fig. 2).The borehole lithological logs were used to evaluate geophysical inversion, conduct MPS simulation, and support results and interpretation.Borehole uncertainty impacts the translation of AEM resistivity to sediment probability and is weighted to account for uncertainty.Borehole data can be weighted based on location accuracy, sampling frequency, drilling method, and borehole geophysical data criteria to account for borehole uncertainty (Barfod et al. 2016;He et al. 2014b;Høyer et al. 2015Høyer et al. , 2017)).Høyer et al. (2017) used a uniform uncertainty of 20% for all boreholes.In this study, data quality is rated using a subdivision into four groups based on location accuracy and borehole log resolution following the decision tree borehole quality assessment method (Korus and Hensen 2020).Location quality rating depends on the source of geographic coordinates.If a borehole location is determined using a global positioning system (GPS) or by measuring N-S and E-W distance from a nearby line (resulting in low location error, typically <10 m), it is considered to have good location quality.Otherwise, it is classified as having poor quality.For the quality of the lithological description, the borehole depth is divided by the borehole logs intervals and if the average is less than 6.1, the description is considered as detailed.If the average is 6.1 or higher, the description is considered not detailed.Using these criteria, the 1,007 boreholes in the study area were classified into four distinct categories.Category 1 boreholes have good location quality and detailed lithological descriptions.Category 1 boreholes were used as hard data because their locations and lithological logs were meticulously recorded by skilled geologists.Category 2 boreholes (a total of 394 boreholes) have detailed borehole lithological logs and good location quality (coordinates measured using GPS or measured footage N-S and E-W nearest sections line).Category 3 boreholes (a total of 387 boreholes) lack detailed borehole lithological logs but have good location quality.Both category 2 and 3 boreholes were not logged by a trained geologist or were unknown, and thus considered soft data (uncertain data).For category 2 and 3 boreholes, a 75 and 65% probability of certainty was assigned to category 2 and 3 boreholes, respectively.The lithohydrofacies from category 2 and 3 boreholes were merged and normalized and integrated with hydrogeophysical facies probabilities obtained from AEM resistivity data, which were used as soft data.A total of 218 boreholes were excluded in category 4 due to their inaccurate and ambiguous location quality to prevent errors from being introduced into the soft data when merging these boreholes with AEM data and utilizing them in MPS simulation.

AEM data processing and inversions
An AEM survey was conducted in 2016 using the SkyTEM304M system (Sørense and Auken 2004) to generate data for hydrostratigraphic units.SkyTEM is a helicopter-borne transient electromagnetic (TEM) method that consists of a transmitter coil and a receiver loop.The TEM method is based on the time-varying magnetic field caused by the rapid turnoff of current sent by the transmitter coil, and it provides subsurface electrical properties via Faraday's law of induction (Christiansen et al. 2006).The induction generates a secondary magnetic field measured at the receiver and contains information about the subsurface resistivity.The data are recorded in time windows called gates and consist of a single sounding at one location.The flight speed and the measurement interval determine the distance between soundings.SkyTEM has two moments: a low moment (3,000 A・m 2 ) enables measurement of early time data and provides details about shallow subsurface structures, and a high moment (160,000 A・m 2 ) enables measurement of late time data, penetrates deeper and offers details about deep subsurface structures (Sørense and Auken 2004).
There are 26 flight lines in the study area (Fig. 2), and the distance between flights is about 500 m.The AEM survey lines run east to west, and one line runs north to south.The AEM survey instrumentation and data acquisition for the study area were presented in AGF (2017).AEM data were processed using Aarhus Workbench 6.2 software which has geophysical and GIS platforms.Manual and automatic methods were applied to adjust altitude and voltage data, as described by Auken et al. (2009).A spatially constrained inversion (SCI) (Viezzoli et al. 2008) was used that forms a quasi-3D model by constraining laterally along and across the flight lines to generate resistivity depth for 30 layers.A 30-layer model was used, with thicknesses increasing logarithmically with depth.The thickness of the first and 29th layers is 3 and 25.9 m, respectively.The depth to the base of the 29th layer is 312 m below the land surface.The Kriging technique with a node spacing of 200 m × 200 m × 3 m was used to interpolate one-dimensional (1D) inverted resistivity to a 3D resistivity grid.The vertical resolution (3 m) was selected to capture vertical heterogeneity in resistivity.

Hydrologic data
Groundwater levels are monitored in a network of two observation wells and 20 irrigation wells in the study area.Monitoring of water levels in irrigation wells is achieved through the use of pipes installed in the filter pack of the well between the well casing and the borehole wall.Automatic groundwater level monitoring systems were installed in July of 2020 to collect daily real-time groundwater-level data.The groundwater level monitoring was aimed to assess drawdown and potential impacts of well interferences during irrigation pumping.Well hydrographs were used to evaluate and verify the interconnections between the aquifer systems and streams.Water-table measurements from spring 2017 were interpolated using the Ordinary kriging method to determine the position of the water table in the aquifer system.There were numerous water level data available from a measurement campaign conducted in the spring of 2017.Water levels data in the spring of 2017 were used to map static water levels because spring is the best time to measure static water levels (no irrigation pumping).

MPS modeling
The MPS work is carried out in this study using the 3D geological modeling software GeoScene3D.The software allows for visualization and integration of different geo-data types and consists of all the necessary tools for a full MPS workflow, from generating the 3D training image through simulation to computing the statistics/hydrofacies probabilities.In this study, the borehole lithological logs were categorized into hydrofacies.Hydrofacies 1 consists of poorly permeable sediments (silt, clay, soil, clay and silt, silt and clay, clay and sand, clay and gravel, shale, till and siltstone) and hydrofacies 2 consists of medium to highly permeable sediments (sand, gravel, silt and sand, sand and gravel, fine sand, coarse sand, gravel and silt).Hydrofacies 2 forms the glacial aquifer, while hydrofacies 1 forms the confining units (aquitard) and bedrock units.Figure 3 depicts the methodology used to simulate aquifer heterogeneity.The simulations are performed using the preferential simulation path technique (Hansen et al. 2018) and the SNESIM algorithm of the open-source software MPSLIB (Hansen et al. 2016) which is embedded in GeoScene3D.

Soft data conditioning
The soft data were prepared in three steps.First, the AEM data were processed, inverted, and interpolated to a 3D grid.Second, borehole data were grouped, ranked, and assigned to the 3D grid.The borehole lithological logs were then classified as lithohydrofacies.Lithohydrofacies 1 consists of fine-grained sediments and lithohydrofacies 2 consists of coarse-grained sediments.This study used lithohydrofacies generated from boreholes with detailed lithological logs and good location quality to determine the corresponding resistivity values for each lithohydrofacies.The lithohydrofacies were assigned to a 200 m × 200 m × 3 m grid and the distribution of resistivity values for each lithohydrofacies was evaluated by examining the correlation between corresponding lithohydrofacies and resistivity values.The correlation between resistivity and lithohydrofacies can be affected by factors such as subsurface water quality, clay minerals, and porosity.In the study area, the effect of porewater conductivity on the relationship is negligible and clay minerals exert a dominant control on resistivity (Korus et al. 2017).The built-in automatic transfer function in GeoScene3D was used to establish the correlation between resistivity values and lithohydrofacies.The transfer function highlights resistivity values of 22 Ω・m as cutoffs for hydrogeoelectrical facies (Fig. 4).A hydrogeoelectrical facies was then generated from AEM resistivity data using this cutoff value.A hydrogeoelectrical facies is a group of sediments that have similar lithological, hydraulic, and electrical properties (Cattaneo 2014).For the successive phases of the workflow, it would be very useful to distinguish hydrogeoelectrical facies based on their resistivity value.Figure 4 shows the probability of the two hydrogeoelectrical facies for the whole range of electrical resistivities in the study area.As resistivity values increase, the percentage of hydrogeoelectrical facies 1 decreases while the percentage of hydrogeoelectrical facies 2 increases.The hydrogeoelectrical facies probabilities were merged with lithohydrofacies from uncertain borehole data and used as soft data.Using information derived from AEM data as soft data leads to a more detailed and representative subsurface realization than relying solely on borehole data (Barfod et al. 2018b).The soft data were used together with TI to determine the probability of hydrofacies at a simulated grid.Using the tau model (Journel 2002), the SNESIM algorithm adjusts the influence of soft data (He et al. 2014a, b;Ma and Jafarpour 2019).A higher value of tau means that the soft data will have a stronger influence, while a lower value of tau means that the soft data will have a weaker influence (Ma and Jafarpour 2019).The value of tau can be estimated from the data.The probability of finding hydrofacies at a specific simulation grid node is calculated based on conditional probability by combining the training image and soft data.One hundred realizations are carried out to ensure that reliable hydrofacies events occur at the simulated grid and to calculate the probability at the simulated location.The number of realizations should be large enough to ensure the stability of the results while also considering the intended processing time of the realizations (Remy et al. 2009).

Hard data conditioning
Hard data conditioning forces MPS realizations to match hard data exactly.In this study, lithohydrofacies generated from eight boreholes that overlap AEM flight lines (Fig. 2) were used as hard data to condition the MPS realizations and reduce uncertainty by forcing the models to match the hard data.

Training image (TI)
In hard-data-scarce areas, MPS realizations rely more on the statistics derived from a TI, which is a cognitive model (Høyer et al. 2017), and from soft data.The ability of the MPS approach to generate realistic representations of geological heterogeneity is heavily dependent on the user-generated TI.TI can be generated from a wide range of data sources; geophysical data such as airborne electromagnetic data (AEM) (He et al. 2013(He et al. , 2014a;;Jørgensen et al. 2015), field data in combination with borehole data (Dall'Alba et al. 2020), and geological background knowledge of the study area.A 3D TI was generated for this study using the same cell size as the 3D AEM resistivity grid but with a smaller footprint (Fig. 5).The aquifer heterogeneity incorporation into the TI was based on a conceptual understanding of the geology of the study area, its aquifer system, glacial processes, interpretations of hydrostratigraphic borehole units, surface geological information, and the AEM model (Korus et al. 2021).Then, expected aquifer heterogeneity in the vertical and horizontal directions was incorporated into the TI.The TI image consists of 106,796 nodes for hydrofacies units.The TI consists of a buried bedrock paleovalley trending W-E, a buried channel trending N-S, and glaciotectonic thrust wedges.The structure and geometry of the glacial deposits in the study area were incorporated into the TI to guide MPS simulation.

MPS parameters and hydrofacies probability calculation
Choosing appropriate MPS parameters, such as search templates and multiple grids, is critical for simulating heterogeneous and discontinuous glacial aquifer systems.The MPS algorithm and the parameters used significantly impact spatial structures derived from the TI (Liu 2006;Høyer et al. 2017).Larger search templates of multiple grids are used to detect large-scale structures and patterns in the TI, while finer templates are used to detect smaller structures within the TI (Mariethoz and Caers 2014).Different MPS model parameters were evaluated to simulate aquifer realizations using the preferential simulation path implemented in Geo-Scene3D.An advantage of the preferential simulation path is that more informed model parameters are visited preferentially to less informed ones (Hansen et al. 2018).The model was simulated with a cell size of 200 m × 200 m × 3 m.
A minimum search template of 3 × 3 × 3 and a maximum search template of 10 × 10 × 10 were used to scan and derive information from the TI.Different numbers of multiple grids were tested and of the multiple grids, six provided good results.The impact of multigrids on MPS realization has been well documented by Liu (2006).
The hydrofacies assigned to the simulation grid are determined using the TI, hard, and soft data.During simulation, hard data provides the first conditioning nodes and patterns to be matched (Madsen et al. 2021).If a grid cell lacks hard data, grid cell-based conditional hydrofacies are calculated by scanning information stored in a search tree and soft data.

Model evaluation
Qualitative evaluations of 3D hydrofacies probability were performed using high-quality borehole logs and heterogeneity features presented in the TI.In addition, well hydrographs were used to determine the connection between aquifers and possible well interferences.The connections between aquifers were assessed following the approach of Butler et al. (2013) and Butler et al. (2021) using several characteristics of the hydrographs, including (1) duration of drawdown during pumping, (2) length of the recovery period after cessation of pumping, (3) year-to-year changes in peak groundwater levels, and (4) slope of hydrograph during nonirrigation months (after drawdown recovery and before the onset of pumping).

Model evaluation
Pseudo wells were placed in the model space at various locations (Fig. 6) to extract attributes from the hydrofacies probability grid and soft data (AEM and boreholes) grids for comparison (Fig. 7).Pseudo wells (Fig. 7) show that the probability of hydrofacies 2 is low at shallow depth (elevation 490-520 m), which corresponds to loess and till (confining unit).The probability of hydrofacies 2 is also low below the bedrock surface (elevation 420 m), which corresponds to shale.The highest hydrofacies 2 probability exists between the elevation of 420-495 m.In general, probabilities of soft data match well with hydrofacies probabilities from MPS simulations at the locations of pseudo wells.High probability (> 0.7) from soft data corresponds to high hydrofacies 2 probability (> 0.6) from MPS simulation in all pseudo wells.

Glacial aquifer system and heterogeneity
Spatial patterns hydrofacies 2 probability varies with depth (Fig. 8).The horizontal slices reveal different zones of hydrofacies 2 probability: low (<0.2),low-medium (from 0.2 to 0.4), medium (from 0.4 to 0.6) and medium to high probability (>0.6).Low hydrofacies 2 probability at elevation 486 m indicates a shallow, clay-rich confining unit (clay, silt, and till units; Fig. 8a).The deepest slice at 386 m shows shale-dominant bedrock units (Fig. 8f).Low-medium hydrofacies 2 probability values indicate heterogeneous units, whereas high hydrofacies 2 probability values indicate coarse sediments deposited in a buried valley and other glacial sand bodies.High hydrofacies 2 probabilities are oriented N-S at elevations of 473, 444 and 432 m (Fig. 8b-d).The interconnectedness of high hydrofacies 2 probability varies spatially: it is localized and discontinuous in the eastern and western parts of the study area.High hydrofacies 2 probability is a productive aquifer composed of sand and gravel deposits in channel deposits and buried valleys.
The horizontal slices (Fig. 8), W-E profiles (Fig. 9), and N-S profiles (Fig. 10) show three types of aquifer connectivity: (1) connected N-S and W-E trending paleovalley   Hydrofacies 2 probability Fig. 6 Locations of pseudo wells, sources of well hydrographs, and profiles between monitoring wells aquifers (Fig. 8c), (2) medium to poorly connected valley aquifer type where heterogeneous units (medium hydrofacies 2 probability) are found between high hydrofacies 2 probability (Fig. 8c-e) and (3) isolated valley aquifer type where the aquifer is bounded by aquitard units (Fig. 8a, b).The depths of confining units (low hydrofacies 2 probability at shallow depth), aquifer zones (high hydrofacies 2 probability), and bedrock (low hydrofacies 2 probability) in the profile correspond well to the lithological borehole logs.The top layer is the main confining unit in the study area, which is underlain by the aquifer materials, as seen in the W-E and N-S profiles (see Figs. 9, 10 and 11).The paleovalley aquifers have a high probability of hydrofacies 2 (composed of permeable coarse-grained sediments).The aquifer materials are not continuous, and there are heterogeneous zones (hydrofacies 2 probability from 0.4 to 0.6).The thickness of the aquifer varies spatially.The heterogeneous zones consist of clay, silt, and coarse-grained sediments such as sand and gravel.The units below the aquifer materials are till and shale units with a low hydrofacies 2 probability (<0.2).
In general, the profiles (Figs. 9, 10 and 11) show the presence of heterogeneous, clay-rich sediments and complex hydrostratigraphic units.The aquifer zone has a variable thickness, limited spatial extent, and is discontinuous.The

Hydraulic conditions and connections between aquifers
All monitoring wells in the study area show groundwaterlevel fluctuations in response to irrigation pumping.During summer months, irrigation pumping occurs intermittently for a period lasting less than 80 days.The magnitude and rate of drawdown during this period, as well as the subsequent recovery of groundwater levels, relate to variable aquifer conditions (confined and unconfined), lateral extents, aquifer-to-aquifer connectivity, and other aquifer characteristics (Butler and Liu 1991;Butler et al. 2013Butler et al. , 2021;;Korus 2018).Differences in the patterns, magnitudes, and rates of drawdown and recovery highlight important spatial variations in the aquifer framework.As such, evidence from well hydrographs is used to qualitatively assess the hydrofacies probability model (Fig. 12).
The hydrograph for well G-032023 is typical of wells completed in the N-S paleovalley aquifer (Fig. 13).Water-level drawdown of 5-8 m occurs over several months.This steady drawdown is punctuated by 3-5-day periods of rapid (about 1 h) drawdown of more than 10 m, followed by similarly rapid recovery.The long-term declines in water level could be a result of the combined drawdown from multiple nearby irrigation wells.
The shorter, rapid periods of drawdown record the development of a seepage face in the filter pack, which is a common occurrence in wells that are screened in unconfined aquifers (Rushton 2006;Chenaf and Chapuis 2007;Behrooz-Koohenjani et al. 2011;Houben 2015a, b).Following the irrigation season, the length of the groundwater-level recovery is only a small fraction of the pumping period (Fig. 13).This rapid recovery is evidence of a laterally constrained aquifer (Butler et al. 2021).Indeed, this well is bounded by low-permeability units in a W-E direction and is located very near the western boundary (Figs.6, and 12a).
Well G-161226 is located on the eastern margin of the same N-S paleovalley aquifer, where it intersects the W-E paleovalley aquifer (Fig. 6).The hydrograph for this well is similar to that of G-032023, except it lacks evidence of filterpack dewatering and its recovery period is slightly longer.G-161226 is bounded by medium (from 0.4 to 0.6) hydrofacies 2 probability, so the longer recovery period could reflect lower transmissivity or lack of boundary influence.The lack of dewatering suggests that the aquifer is confined, which agrees with the hydrofacies model (Fig. 12a).
The well hydrographs from G-147342 and G-161227 (Fig. 13) are typical of wells screened in the W-E paleovalley aquifer.The hydrographs show drawdown varying between 8 and 15 m and recovery times of about 50 days.The hydrofacies probability shows a high hydrofacies 2 probability (>0.7) between these two wells and throughout the 6-km-wide aquifers, which matches with coarse sediments shown in borehole lithological logs (Fig. 12b).It also shows a confining layer above the aquifer.The groundwaterlevel fluctuations are consistent with the interpretation of confined conditions, interconnections between wells in the paleovalley, and minimal influence of aquifer boundaries.
Well G-133348 lies in a zone of significant heterogeneity east of the W-E paleovalley aquifer (Fig. 12b).Groundwater level declines and recovery times differ significantly from those described previously.Recovery is incomplete following each irrigation season.Annual peak water levels in the well have declined by about 6 m over the past 6 years.Hydrofacies 2 probability is low (0.2-0.4) and the aquifers in this area are poorly connected to the W-E paleovalley aquifer (Figs. 8 and 12b).Borehole logs also show heterogeneous and thin coarse sediments in the vicinity of G-133348.
Monitoring wells G-149812 and G-051317 are located in the middle part of the N-S paleovalley aquifer (Fig. 6).The  hydrographs show evidence for filter-pack dewatering, indicating they tapped an unconfined aquifer.The hydrographs (Fig. 14) show that the aquifer has fully recovered before the next irrigation season.Eastward, near Loseke Creek (Fig. 12c), the water-level surface drops abruptly, indicating a hydraulic boundary.The location of this boundary corresponds to a heterogeneous zone separating the E-W paleovalley from deposits farther east.On the easternmost side of this profile, well G-173012 shows incomplete recovery between irrigation seasons, similar to other wells in this heterogeneous zone.
The well hydrograph from G-126483 is typical of wells in the heterogeneous zone west of the N-S paleovalley (Fig. 14).It shows a high drawdown (~12 m) and rapid recovery (~5 days; Fig. 14).It also shows rapid dewatering of the filter pack.This well hydrograph is characteristic of a thin unconfined aquifer with low transmissivity and enclosed by hydraulic boundaries.The hydrofacies model (Fig. 12d) shows a thin high hydrofacies 2 probability (>0.6) bounded by heterogeneous hydrofacies 2 probability (probability ranges between 0.2 and 1.0).It also shows a poor connection with the N-S paleovalley aquifer (Fig. 12d).
The hydrographs for wells G-126483 and G-157921, which are located near the stream outlets (Fig. 12e, f), are unlike any hydrographs in the upstream reaches (Figs. 13 and 14).Both wells show evidence of pumping-induced drawdown, but water levels recover rapidly and completely.There is little to no long-term drawdown of the water table in the vicinity of these wells during the irrigation season, despite evidence of large-magnitude drawdown in wells G-164613 and 173012, which are located just a few kilometers away.These observations strongly suggest that the aquifers in the downstream area are compartmentalized and are affected by recharge boundaries.It also suggests that there may be localized pathways for rapid recharge in upland areas between streams.

Delineation of groundwater management zones
Seven groundwater management zones were developed based on hydrofacies probability and hydrograph characteristics (Fig. 15).Table 1 shows the detailed characteristics of each zone as well as the management issue.Zone 1 is an unconfined aquifer, while zone 2 is heterogeneous and has abrupt variations between confined and unconfined aquifers.The main management issue in zone 2 is local well interference.The horizontal slice (Fig. 15, zones 3-4) demonstrates the connection between N-S and E-W-oriented palaeovalley aquifers.The decline in groundwater level in connected valley aquifers is most likely the result of multiple well interferences caused by nearby irrigation well pumping (Fig. 15, zone 4).Well interferences increase groundwater-level drawdown and affect well yield for crop production.Zone 3 indicates a narrow, bounded, unconfined aquifer with high hydrofacies 2 probability.Excessive drawdown can be seen in the hydrographs of the wells in zone 3 located near the boundaries.Zones 3 and 4 can be managed as a single aquifer.Future well permits and siting should consider the effects of well interferences and sustainable yield in these zones.Variability in groundwater level declines and recovery in nearby wells indicate a weak connection between aquifers (Fig. 15, zones 5-6).These zones are characterized by extreme heterogeneity that affects direct aquifer recharge.Well permitting in zones 5 and 6 should consider the impact of high drawdown during      Legend are hydrologically connected only in one small area.This heterogeneity presents a major challenge for groundwater management because major changes in hydrologic behavior occur over short distances.MPS provides a framework for delineating proposed management boundaries (Fig. 15) and hydrographs reveal the hydrologic behaviors specific to each zone (Fig. 14).Thus, the workflow presented here can be used to create tailored management solutions for highly heterogeneous aquifers.Fluctuations and recovery in groundwater levels are influenced by groundwater pumping, well interference and the nature of the aquifer.The heterogeneity of the aquifer obstructs groundwater flow to wells and if the cone of depression intersects with impermeable layers, significant drops in water levels occur.Some well hydrographs in zone 5, which are typical of wells that extract water from aquifers surrounded by impermeable formations, exhibit abrupt decreases in groundwater levels during irrigation pumping (Korus 2018).The hydrofacies probability profiles (Figs. 11 and 12e) show the existence of a localized aquifer.This aquifer is exposed to the land surface that may serve as a favorable location for the implementation of managed aquifer recharge techniques, such as an artificial recharge basin (Knight et al. 2022;Pepin et al. 2022;Uhlemann et al. 2022).
The use of appropriate aquifer and streambed heterogeneity could help to reduce uncertainty in groundwater-stream interaction modeling.The discontinuity and existence of low hydraulic conductivity affect groundwater flow (Åberg et al. 2021).In a heterogeneous aquifer system, the use of a classical model that involves layering a confining-layer system may result in errors in the groundwater model.The results from hydrofacies probability improve hydrostratigraphic conceptualizations and aid in parameterizing groundwater flow models for model calibrations.The hydrofacies 2 probability can be clustered and the relationship between the clustered hydrofacies 2 probability and hydraulic conductivity can be assessed to parameterize the groundwater model.The K-means clustering method can be used to group hydrofacies probabilities, which can then be used to estimate the hydraulic conductivity for each group.Furthermore, K-means clustering can be used to generate hydrofacies probability zones, and inverse groundwater modeling can be used to estimate hydraulic conductivity for each clustered zone.Marker et al. (2015) used K-means clustering to categorize voxel models of electrical resistivity and clay fraction into hydrostratigraphic zones and estimated hydraulic conductivity for each zone using hydrological calibration.
The heterogeneity of streambed sediment controls the water exchange between streams and groundwater.The hydraulic conductivity of the streambed is more variable in meandering rivers than in straight channels (Abimbola et al. 2020).Fleckenstein et al. (2006) compared stream and groundwater interactions in homogeneous and heterogeneous aquifers and reported that facies distribution affects flux and groundwater levels.Kalbus et al. (2009) investigated the influence of aquifer and streambed heterogeneity on groundwater discharge distribution and reported that aquifer heterogeneity has a larger effect on stream-groundwater interactions than streambed heterogeneity.MPS hydrofacies simulation can help to assess streambed hydraulic conductivity and thickness which are vital for assessing groundwater-stream interaction modeling.The correlation between MPS-generated hydrofacies probability along streams and hydraulic conductivity can be used to assess the spatial variations of hydraulic conductivity of the streambed and streambank.This can help to reduce uncertainty in the results of stream-aquifer interactions modeling.In addition, MPS hydrofacies probability can be used in conjunction with stream coring data to estimate the spatial variation of streambed thickness for the stream-groundwater interactions model.In the upstream parts of the study area, thick confining units are found between the streambed and aquifer.The hydrofacies 2 probability indicates thin clay-rich units between the streambed and the aquifer, indicating no direct connection between the streams and the aquifer in the middle parts of Loseke and Taylor creeks.This complicates water infiltration through streambeds and causes a delay in stream and groundwater exchange during groundwater pumping, prolonging stream depletion.Near all stream outlets (see Figs. 12e and 15, zone 7), there is a high hydrofacies 2 probability that is exposed to the land surface.This permeable unit can lead to stream depletion as the groundwater level in the aquifer decreases during pumping, or it can provide groundwater contribution to streams.The effects of irrigation pumping on streamflow and ecosystems should be considered in areas where the streambed and aquifer are connected.In general, the streambed hydrofacies probability can be extracted and the relationship between streambed hydraulic conductivity and hydrofacies probability can be established to estimate spatial variability in streambed conductivity.

Conclusions
In this research, glacial aquifer heterogeneity was simulated to support the groundwater management system in the Shell Creek watershed of eastern Nebraska, USA.Borehole data, training images (TI), and airborne electromagnetics (AEM) were used to simulate aquifer heterogeneity realizations and compute hydrofacies probability at each simulated node of 200 m × 200 m × 3 m.Monitored groundwater heads were used to validate aquifer heterogeneity modeling.The following conclusions were drawn based on the expected value of hydrofacies 2 probability at the simulated node.
• The principal aquifers in the area are two buried valley aquifers that intersect at a right angle.These aquifers are marked by high hydrofacies 2 probability and comparatively high spatial continuity.• The heterogeneous aquifers consist of hydrofacies 2 probability between 0.4-0.6 and have a variable thickness.• The paleovalley aquifers are disconnected in some places and bounded by low hydrofacies 2 probability, which can cause significant drawdown when a cone of depression intersects this low hydrofacies 2 probability area.• Thin layers are resolved well in the vicinity of hard data, resulting in an abrupt change of aquifer framework near these grid nodes.• Hydrofacies 2 probability shows a poor connection between streambank, streambed and aquifer at upstream reaches of the streams.In downstream reaches near the stream outlets, however, there is medium to high (> 0.6) hydrofacies 2 probability between the aquifer and streambed.
• The well hydrographs show that groundwater responds differently to pumping in upstream and downstream reaches, particularly water levels recover faster near stream outlets in some wells than the monitoring wells located upstream.• The 3D hydrofacies 2 probability at the simulated nodes and well hydrographs are essential for delineating groundwater management zones, prioritizing future well-monitoring locations and identifying a potential site for managed aquifer recharge.

Fig. 1
Fig. 1 Locations of a the State of Nebraska (NE) and the glaciation ice limits in the USA, and b the study area in Nebraska.c Geographic features of the study area, including a digital surface elevation model (DEM, m above sea level)

Fig. 2
Fig. 2 Airborne electromagnetics (AEM) flight lines, Nebraska Conservation and Survey Division (CSD) test holes, Nebraska Department of Natural Resources (DNR) boreholes, and W-E and N-S

Fig. 3
Fig. 3 Multiple-point statistic (MPS) simulation workflow used in this study

Fig
Fig. 4 The transfer function illustrating the probability of hydrogeoelectrical facies

Fig. 5
Fig. 5 Training images: a three-dimensional b W-E and N-S slices c only hydrofacies 2 presented Fig.7a

Fig. 7
Fig. 7 Comparison of hydrofacies 2 probability from MPS and soft data probability; the borehole locations are shown in Fig. 6 as "Fig.7a-e"

Fig. 8
Fig. 8 Horizontal slice view of hydrofacies 2 probability from MPS a 486 m b 473 m, c 444 m, d 432 m, e 410, f 386 m

Fig. 9
Fig. 9 W-E profiles of hydrofacies 2 probability in the northern part of the study area at different locations.Test holes (th) are hard data; Borehole category 2 (ct2) and borehole category 3 (ct3) are soft data.
Fig. 10 W-E profiles of hydrofacies 2 probability in the middle and southern (downstream reaches) of the study area at diffrent locations; Legend presented in Fig. 9. Test holes (th) are hard data; Borehole

Fig. 12
Fig. 12 Profiles of probability of hydrofacies from MPS between groundwater monitoring wells, a) A-A' b) B-B' c) CC' d) D-D' e) E-E' f) F-F'.Location of the profile is shown in Fig. 6.Legend presented in Fig. 9.The water level refers to spring 2017

Fig. 13
Fig. 13 Hydrographs from observation wells showing drawdown and recovery 2018-2022 in the upstream reaches of the study area.See Fig. 12a, b for locations of wells

Fig. 14
Fig. 14 Hydrographs from observation wells showing drawdown and recovery during 2018-2021 in the middle and downstream reaches of the study area.See Fig. 12c-f for locations of wells

Fig
Fig. 15 Proposed management boundaries with respect to hydrofacies 2 probability, represented as a horizontal slice through the model at 444 m elevation

Table 1
Management zones, recovery, hydrogeologic characteristics and management issue(s)