Quaternary glacial aquifers are important water sources for irrigation in many agricultural regions, including eastern Nebraska, USA. Quaternary glacial aquifers are heterogeneous and comprise sediment assemblages with a wide range of hydraulic properties. Effective management and sustainable utilization of these heterogeneous glacial aquifers necessitate the development of realistic groundwater-flow models and characterization of aquifer geometry. However, hydrofacies probabilities predicted using multiple-point statistics (MPS) and machine learning (ML) techniques are rarely used for parameterizing groundwater models and identifying management zones. This study used MPS to simulate 100 three-dimensional conditional aquifer heterogeneity realizations by combining soft data, a cognitive training image, and hard data. The most probable hydrofacies model (sand and clay probability) was then calculated at a 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. A new workflow was established using the three-dimensional hydrofacies probability generated by MPS and hydrologic data to define high-resolution groundwater management zones and enhance strategies. Subsequently, ML techniques such as Random Forest, Gradient Boosting Classifier, Extreme Gradient Boosting, Multilayer Perceptron, and Stacking Classifier were used to model three-dimensional probabilistic distributions of hydrofacies at a grid size of 200 m x 200 m x 3 m. The models were compared in terms of their capability to identify thin, permeable hydrofacies, lateral continuity, and vertical contrast between hydrofacies units. ML predicted thin, permeable, and laterally continuous hydrofacies better compared to hydrofacies predicted by MPS. Multilayer Perceptron and Stacking Classifier models show sharper vertical contrasts between fine and coarse hydrofacies compared to MPS and other ML models. Finally, three groundwater models were constructed using MODFLOW 6 with unstructured grids. The first model was parameterized using hydraulic conductivity derived from pumping and geological data, while the second and third models were parameterized using hydraulic conductivity estimated from MPS and stacking machine learning hydrofacies models, respectively. K-means clustering was used to translate the predicted hydrofacies probability into hydraulic conductivity values. While the entire water budgets of all models show minimal variation, zonal water budget analyses reveal significant differences in storage change, stream-groundwater interactions, and total inflow and outflow. This study effectively demonstrates the influence of three-dimensional aquifer heterogeneity modeling approaches on the outcomes of groundwater models. Such insights can prove invaluable for groundwater managers and policymakers in assessing the implications of groundwater model parameterizations on local groundwater management.