Capturing soil-water and groundwater interactions with an iterative feedback coupling scheme: New HYDRYS package for MODFLOW

Abstract

Accurately capturing complex soil-water and groundwater interactions is vital for describing the coupling between subsurface/surface/atmospheric systems in regional-scale models. The non-linearity of the Richards' equation for water flow, however, introduces numerical complexity to large unsaturated-saturated modeling systems. An alternative is to use quasi-3D methods with a feedback coupling scheme to join practically sub-models with different properties, such as governing equations, numerical scales, and dimensionalities. In this work, to reduce the non-linearity in the coupling system, two different forms of 10 the Richards' equation are switched according to the soil-water content at each numerical node. A rigorous multi-scale water balance analysis is carried out at the phreatic interface to link the soil water and groundwater models at separated spatial and temporal scales. With a moving-boundary approach at the coupling interface, the non-trivial coupling errors introduced by the saturated lateral fluxes are minimized for problems with dynamic groundwater flow. It is shown that the developed iterative feedback coupling scheme results in significant error reduction, and is numerically efficient for capturing drastic flow 15 interactions at the water table, especially with dynamic local groundwater flow. The coupling scheme is developed into a new HYDRUS package for MODFLOW, which is applicable for regional-scale problems.quasi-3D methods are categorized into (1) the fully coupling scheme, which simultaneously builds the nodal hydraulic 50 connections of models at both sides and implicitly solves the assembled matrices;(2) the one-way coupling scheme, which delivers the soil-water model solutions onto the upper boundary of the groundwater model without feedback mechanism; and (3) the feedback (or two-way) coupling scheme, which explicitly exchanges the head/flux solutions in vicinity of the interface nodes.The fully coupling scheme (Gunduz and Aral, 2005; Zhu et al., 2012) is numerically rigorous but tends to increase the 55 computational burden for practical conditions. For example, the potential conditional diagonal dominance causes nonconvergence for the iterative solvers (Edwards, 1996). Owing to high non-linearity in the soil-water sub-models, the assembled matrices can only be solved with unified small time-steps, which adds to the computational expense. The one-way coupling scheme, as adopted by the MODFLOW-UZF1 package (Grygoruk et al., 2014; Niswonger et al., 2006), as well as the free drainage mode of MODFLOW-SWAP model (Xu et al., 2012), assumes that the water table depth is of minor influence on 60 flow interactions at the phreatic interface, and is thus problem specific.The feedback coupling method, in contrast, is widely used (Kuznetsov et al., 2012;Seo et al., 2007;Shen and Phanikumar, 2010;Stoppelenbrug et al., 2005;Xie et al., 2012;Xu et al., 2012) as a compromise between numerical accuracy and computational cost. In a feedback coupling scheme, the soil-water and groundwater sub-models can ...

Coupling saturated and unsaturated flow: comparing the iterative and the non-iterative approach

“…Usually, they are linked by the exchange of information on groundwater recharge and the position of the groundwater table. This can be done either in an iterative manner to account for the feedback between the two models (Shen and Phanikumar, 2010;Zhu et al, 2011;Kuznetsov et al, 2012;Xie et al, 2012;Xu et al, 2012;Mao et al, 2019;Zeng et al, 2019) or in a non-iterative manner without feedback (Pikul et al, 1974;Yakirevich et al, 1998;Beegum et al, 2018;. Both methods have been applied successfully.…”

“…Both methods have been applied successfully. While iterative coupling in general yields a higher accuracy, non-iterative coupling requires less computational ressources (Zeng et al, 2019).…”

“…Differences among the coupling strategies exist also in terms of the spatial extent of the unsaturated zone model. Mostly, the unsaturated zone is modeled to span from the groundwater table to the land surface (Pikul et al, 1974;Yakirevich et al, 1998;Zhu et al, 2011;Kuznetsov et al, 2012;Xie et al, 2012;Zhu et al, 2012;Renard and Tognelli, 2016;Zeng et al, 2019). As the groundwater table is moving with time, a remesh is needed after each coupling time step.…”

“…With an iterative coupling scheme, consistency can be achieved, but there still remains an error due to neglection of the temporal changes in the exchange information. Zeng et al (2019) use linear extrapolation of the groundwater table position to reduce this error. Zhu et al (2012) solve the two models jointly linked via the pressure head values at the interface.…”

“…Beegum et al (2018) performed a small sensitivity analysis on this parameter and detected a strong influence on the dynamics of the groundwater table position. Some research was dedicated to getting a better estimation of this value in the context of groundwater modeling (Hilberts et al, 2005) and coupled saturated-unsaturated zone modeling (Pikul et al, 1974;Xu et al, 2012;Zeng et al, 2019). Hilberts et al (2005) presented an analytical expression for the specific yield depending on the groundwater table position and the soil hydraulic parameters assuming hydrostatic equilibrium in the unsaturated zone.…”