Numerical models of groundwater flow in folded and faulted aquifers in mountainous areas: trade-offs between numerical error and runtime

Groundwater flow in folded and faulted terrains is governed by the geological structure, which exerts a significant influence on the flow directions and on the spring location.  Numerical modeling can be challenging because the complex geometry of model layers can result in steep inclination of aquifer bed and different thicknesses of saturated portions within the same numerical layer. Drying and rewetting of cells during model iterations leads to numerical instabilities, increasing numerical error and runtime. However, excessive simplification of the system, seeking for numerical stability, may lead to an unsatisfactory predictive capability of the model. Recent advances have addressed such issues, including the development of solvers that facilitate convergence and/or reduce computational errors due to model nonlinearities [1], [2]. The aim of this research was to develop and test a procedure for the simulation of groundwater flow in a complex karst, folded, multilayer aquifer, while minimizing numerical errors and runtime. We applied this procedure to a 3D model, previously developed in steady state conditions with an equivalent porous media approach using the Newton-Raphson formulation of MODFLOW-2005 (MODFLOW-NWT) [3]. The major impact of folded and faulted geological structures on controlling the flow dynamics in terms of flow direction, water heads, and spatial distribution of the outflows to the river and springs was accounted for in a numerical model where three aquifer layers and two semipermeable layers have been constructed respecting their true geometry as far as possible.  A transient simulation was performed on this model using monthly stress periods and variable pumping simulated by MODFLOW’s WEL package to test effects of withdrawals for water supply on the aquifer system. Initial runs showed a very high mass balance error (2% discrepancy over cumulative volume) and a runtime of 1 hour and 38 minutes. To reduce both mass balance error and runtime, the USGS software for the optimization of the MODFLOW-NWT solver inputs, NWTOPT [4], was used.  NWTOPT identified improved solver inputs, which gave a superior tradeoff between acceptable mass balance error (-0.39%) and much reduced runtime (40 minutes and 9 seconds). The added stability and shorter runtimes are particularly welcome if many iterations of the model were needed for automated calibration, application or uncertainty analysis. References [1]  Niswonger, R. G., Panday, S., & Ibaraki, M. (2011). MODFLOW-NWT, a Newton formulation for MODFLOW-2005. US Geological Survey Techniques and Methods, 6(A37), 44. [2] Hunt, R. J., & Feinstein, D. T. (2012). MODFLOW-NWT–Robust handling of dry cells using a Newton Formulation of MODFLOW-2005. Ground Water, 50(5), 659-663. [3] Preziosi, E., Guyennon, N., Petrangeli, A.B., Romano, E., Di Salvo, C. (2022) A stepwise modelling approach to identifying structural features that control groundwater flow in a folded carbonate aquifer system. Water, 14 (16), art. no. 2475,  DOI: 10.3390/w14162475 [4] Newcomer, M. W., & Hunt, R. J. (2022). NWTOPT–A hyperparameter optimization approach for selection of environmental model solver settings. Environmental Modelling & Software, 147, 105250.