Surface and Groundwater Interactions: A Review of Coupling Strategies in Detailed Domain Models

Abstract

In groundwater numerical simulations, the interactions between surface and groundwater have received great attention due to difficulties related to their validation and calibration due to the dynamic exchange occurring at the soil–water interface. The interaction is complex at small scales. However, at larger scales, the interaction is even more complicated, and has never been fully addressed. A clear understanding of the coupling strategies between the surface and groundwater is essential in order to develop numerical models for successful simulations. In the present review, two of the most commonly used coupling strategies in detailed domain models—namely, fully-coupled and loosely-coupled techniques—are reviewed and compared. The advantages and limitations of each modelling scheme are discussed. This review highlights the strategies to be considered in the development of groundwater flow models that are representative of real-world conditions between surface and groundwater interactions at regional scales.

1. Introduction

Water security is an emergent issue in the world due to climate change and variability, and increasing water demands due to population growth and economic development [1]. Both surface water and groundwater are, inevitably, resources that are needed in order to meet future water demands. Groundwater is widely used in the national economies of many countries for different purposes, such as potable and industrial water supply, irrigation, and mineral water. Fresh groundwater is of particular importance as the main source of public water in the domestic and drinking water balance. To date, groundwater is the main source of potable water in most European countries and in the United States, which accounts for 75% of municipal water supply system [2]. Irrigation has grown rapidly over the past 50 years, and approximately 40% of the irrigated regions of the world are predominantly governed by groundwater resource [3]. However, in recent years, the depletion of groundwater has become a concern due to overexploitation. A greater understanding is needed to better manage the limited water resources, especially in drought-prone areas. There is a need to better understand the interactions between groundwater and surface water, which is fundamentally an integrated feature in the hydrologic cycle. These interactions could be influenced by both natural and anthropogenic processes [4]. Recently, research has focused on the proper selection of models that accounts for surface and groundwater interactions coupled into a single framework. However, there are many restrictions associated with this approach. There are considerable fluctuations of surface and groundwater interactions at the soil–water interface due to the heterogeneous nature of the porous medium with time, along with the variabilities due to climate change. The diversity of aquifer types, along with recharge and groundwater flow conditions, and the extensive range of river flow circumstances provide the ingredients for a complex system of interactions, which heightens the need for robust modelling and solution methodologies [5].

Surface water models are being used to understand surface water systems by including possible deviations due to both natural and human influences. The interrelation and coupling among various subsystems of the hydrological cycle has been noted by many pioneer researchers in the process of the physical and mathematical description of the flow processes [6]. The first physically-based and coupled model, Système Hydrologique Européen (SHE), was developed primarily based on fundamentals of surface water flow. The United States Geological Survey (USGS) PRMS (Precipitation-Runoff Modelling System) and National Oceanic and Atmospheric Administration (NOAA)’s National Weather Service models [7] are the recent developments of surface water models. The groundwater flow models have the ability to simulate predicted impacts from over-exploitation, to simulate the fate and transport of a contaminate plume, and to forecast the variation of the groundwater levels with time due to climatic variability. However, the traditional groundwater flow models neglect the incorporation of surface water flux into the subsurface environment. MikeSHE and SHETRAN were further developed and applied with both surface and groundwater flow [8,9,10].

From the literature, it has been noticed that the modelling and simulation of groundwater and surface water interaction have been conducted by many researchers [11,12,13,14]. Multiple factors have been reported that could influence the interaction between ground and surface water, such as hydro-climatic variables; physiographic structure; differences in the head between the catchment surface water and groundwater; and the surface and groundwater flow geometry [15,16]. The modeling of flows and transport phenomena as discrete groundwater systems has been achieved well by relatively simpler approaches [15,17]. However, the modelling process is complicated by the inclusion of surface water and groundwater interaction. Ever-present interactions between groundwater and surface water are a concern due to the continuous flow of groundwater into rivers [5]. The integration of groundwater and surface water interaction could be classified in many different ways. Figure 1 illustrates a schematic of the different available modelling approaches. It serves as the foundation and framework for the present review. The objectives of the present review are to compare the coupling strategies in detailed domain models, which are the fully-coupled and loosely-coupled schemes. The domain integrated modelling will be briefly reviewed and compared to detailed domain models.

2. Detailed Domain Models

Detailed domain models have been developed on surface and groundwater simulation. For the regional scale modelling, the fully-coupled and loosely-coupled schemes are the two widely used techniques. The traditional methodologies and the hierarchical modelling approach will be addressed within the local scale modelling discussion.

2.1. Regional Scale Modelling

A regional scale is defined as a catchment with an area between 1000 and 10,000 km². It can include variables such as climate, geology, geomorphology, landscape categories, and biological and human factors in the same region [18,19,20]. The surface and groundwater interaction on the regional scale incorporates the entire terrestrial hydrological cycle using observations at larger distances and simplified interaction processes. Fully- and loosely-coupled schemes are the two widely used techniques in regional scale modelling.

2.1.1. Fully-Coupled Scheme

Fully-coupled schemes are physically-based models that integrate the surface and groundwater components of the hydrologic cycle into a single software package. A fully-coupled scheme software package has the ability to solve governing equations for both surface and groundwater flow simultaneously. Fully-coupled models can investigate the interaction over a range of spatial and temporal scales. Stable results can be acquired by these models [7], resulting in a reduction in the number of ancillary software packages needed to address the various subsystems [21,22]. In fully-coupled systems, the governing equations for groundwater flow are based on the modified Richards’ equation for both saturated and unsaturated zones. The governing equation for surface water is based on the depth-averaged Saint-Venant equation given below.

Governing Equations

Richards’ equation describes the conservation of the mass of the water phase in the void space in the porous medium, along with Darcy’s equation for the volumetric flux given below.

• Porous Medium (3D) Richards’ Equation:

$$w_{m }\frac{\partial {\theta }_{s }{S}_w}{\partial t}+\nabla {.w}_m\overrightarrow{q}\pm {{\displaystyle \sum }}^{}{\Gamma }_{ex}=\pm Q$$

(1)

• Darcy–Buckingham Equation:

$$\overrightarrow{q}=-\stackrel{\rightrightarrows }{K}. k_r\nabla \left(\psi +z\right)$$

(2)

where, $w_m$ (dimensionless) is the volumetric division of the total porosity of the porous medium (or the primary continuum); $Q$ = source/sink rate, which defines the volumetric fluid flux per unit volume demonstrating a source (positive) or a sink (negative) to the porous medium system$; {\Gamma }_{ex}$ = exchange fluxes term$; S_w$= degree of water saturation; ${\theta }_{s }$= porosity; $\psi$ = pressure head; z = elevation head; $\stackrel{\rightrightarrows }{k}$ is the permeability tensor; and $k_r$ (dimensionless) is the relative permeability (function of saturation). The hydraulic conductivity tensor, $\stackrel{\rightrightarrows }{K}$ [$L{T}^{-1}$] = $\frac{\rho g\stackrel{\rightrightarrows }{k}}{\mu },$ where $\mu$ is the water viscosity [$M{L}^{-1 }{T}^{-1}$], g is the acceleration due to gravity [$L{T}^{-2}$], $\stackrel{\rightrightarrows }{k}$ is the permeability tensor of the porous medium [L²], and $\rho$ is the water density [$M{L}^{-3}$], which is a function of the concentration C [$M{L}^{-3}$] of any particular solute, such that $\rho =\rho \left(C\right)$. The water saturation is interconnected to the moisture content $\theta$ (dimensionless) according to $S_w=\frac{\theta }{{\theta }_s}$. In Equation (2), $\Gamma$_ex represents the exchange rate of the volumetric fluid $\left[T^{-1}\right]$ between the subsurface domain and all of the other categories of domains supported by the model, and it is expressed per unit volume of the other domain types. Generally, these additional domains are surface, tile drains, wells, discrete fractures, and dual continuum. The definition of $\Gamma$_ex (positive for the flow into the porous medium) relies on the conceptualization of fluid exchange between the domains [23].

The depth-averaged Saint-Venant equation for surface water can be written as:

$$\frac{\partial {\varnothing }_0{h}_0}{\partial t}+\nabla .d_0{\overrightarrow{q}}_o+d_0{\Gamma }_{0 }=\pm Q_0$$

(3)

where $Q_0$= source/sink rate, $d_{0 }$= water depth; ${\Gamma }_{0 }$= exchange fluxes; h₀ = surface hydrostatic pressure. Incorporating the Manning equation relating the flux with the gradient of the surface water system, it can be written as:

$${\overrightarrow{q}}_0= -\frac{d_0^{2/3}}{{\eta }^{ }{\varphi }^{1/2}}{k}_{r0}\nabla (d_0+z_0)$$

(4)

where $\eta$ is roughness, $\varphi$ is the surface water gradient, and $k_{r0}$ = relative permeability. The surface and subsurface flows are coupled through interaction terms that represent the flux interchange between the two surface–subsurface compartments [23,24].

$$-\nabla .{\omega }_m\overrightarrow{q}+{\Gamma }_{0 }\pm Q=w_m\frac{\partial {\theta }_s{S}_w}{\partial t}$$

(5)

$$-\nabla .d_0{\overrightarrow{q}}_0-d_0{\Gamma }_{0 }\pm Q_0=\frac{\partial {\varnothing }_0{h}_0}{\partial t}$$

(6)

It has been hypothesized that the first-order exchange term can be calculated using

$$d_0{\Gamma }_{0 }=k_{rso}{K}_{so}\left(h-h_o\right)$$

(7)

where h is the subsurface pressure head, $K_{so}$ is the conductance at the surface/subsurface interface, and $k_{rso}$ is the coupling-relative permeability/rill storage.

Fully-coupled software packages are based on three separate coupling strategies for the integration of hydrostatic surface water (i.e., Saint-Venant) and subsurface flow: first-order exchange, continuity of pressure and boundary condition switching [25]. A first-order exchange formulation or conductance has been applied for the coupling of surface and groundwater flow equations in HydroGeoSphere (HGS) and InHM software packages. The flux continuity can be maintained, and the continuity of pressure can then be imposed across the surface and subsurface domains [26,27]. The first-order exchange approach is based on Darcy’s Law for flow through an exchange interface. The surface and groundwater exchange flux (positive/negative as exfiltration/infiltration) is linearly dependent on the difference between the subsurface head at the uppermost node, and the surface head, as well as the value of the first-order exchange or conductance coefficient [27]. This mechanism has been used in investigations pertaining to flow simulation between different continua, such as fractures and macropores, rock and soil, and river–aquifer interactions through streambeds.

In another coupling strategy, the continuity of pressure coupling formulation is introduced and incorporated into the ATS (Advanced Terrestrial Simulator), Cast3M, and ParFlow software packages. These packages are fully-coupled models with a free-surface overland boundary condition on the top, based on pressure and flux continuity at the surface [28,29,30]. The surface water equations were used for the interface conductance [31].

The boundary condition switching coupling strategy has been used in CATHY software to enforce flux and pressure continuity at the surface/subsurface interface. This approach acts as a special treatment in variably-saturated subsurface flow models to detect atmosphere-controlled and soil-limited infiltration and evaporation dynamics [28]. It translates potential atmospheric fluxes into actual fluxes across the land surface, and the resulting surface storage [32,33]. A common hypothesis in these approaches is to maintain pressure and flux continuity at the surface/subsurface interface while ignoring the continuity of momentum and forces.

A limited number of regional applications using fully-coupled models can be found in the literature, including its applications on small scales [34]. As one of the fully-coupled models, the Penn State Integrated Hydrologic Model (PIHM) was developed with a reservoir simulation module. PIHM can be converted into an optimization model under the Simulation-Optimization (S-O) modelling structure [35]. The Simulation-Optimization (S-O) modelling approach has been extensively recognized to resolve complications in water resources management [36,37,38,39,40,41,42,43,44]. The main idea of S-O modelling is to build interconnectivity between a simulation model and optimization scheme so that a trade-off between groundwater extraction rates and surface water reduction can be established. A fully-coupled model was also tested for lake-dominated hydrologic regimes with natural fluctuations, such as landscape disturbance. This integrated approach has been applied in spatially-distributed atmospheric flux information, including evapotranspiration, throughfall, and landscape diversity, in order to imitate field observations. However, it is recommended that the groundwater recharge rate and timing should be pre-defined between evapotranspiration and temporal throughfall fluid. Hence, groundwater recharge is being considered as a governing factor for the maintenance of lakes and wetlands by the use of the fully-coupled method [45].

The advantages of the fully-coupled scheme are also very clear. It has been found that the simulation results have a good agreement with the experimental/observed field data. The fully-coupled models can also be used for coastal and lake simulations. Additionally, this technique has been validated based on the large number of available software packages from previous studies. However, as noted, the fully-coupled approach should not be used in unsaturated zones with overbearing lateral flow [46]. In order to improve the accuracy of fully-coupled models, it is suggested that the quasi-adaptive meshing configuration should be adopted without considering the total number of grids. The quasi-adaptive meshing offers flexibilities to correct local discretization inaccuracies due to the topographic elevation grade [47]. In short, the fully-coupled models can serve as a strong watershed simulator for the examination of watershed hydrological characteristics with the desired accuracy.

Many pieces of software can be found using fully-coupled schemes. The dynamic interactions between the surface and subsurface are captured in a seamless way. However, the new parameters that appear in the flow equation with the coupling term should be determined experimentally or by fitting exercises.

Table 1. Fully-coupled software.

Available Software	Review Articles References
HydroGeosphere	[23,48]
ParFlow	[26,30,49,50]
OpenGeoSys	[51,52]
CATHY	[53,54]
InHM	[34]
MIKESHE-2003 version	[17]
HYDRUS	[17]
SUTRA	[16,17]
PAWS	[25,28]
PIHM	[25]
tRIBS + VEGGIE	[25]
IWFM	[55]
Cast3M	[28]
Community Land Model version 4.5	[56]

lists examples of software packages that are fully-coupled models.

With the application of the fully-coupled schemes, there are some disadvantages to the approach. There has been a lack of studies in the application of the fully-coupled approach at a regional scale. This limitation is likely due to the robustness of the software packages addressing the surface–groundwater interactions at the catchment scale. A direct comparison of the simulation results using different software is not practical. Most of the evaluations have been completed on local scales with limited testing and applicability to the real-world environment. The application of fully-coupled approaches on larger scales has not been extensively tested due to the complexity of the interactions between surface water and groundwater [34]. Another disadvantage of the fully-coupled models is the lengthy computational time. The calibration of the model is another challenge due to the large number of estimated parameters. The boundary and the initial conditions could also introduce larger uncertainties, resulting in high computational effort and cost. It is critical to scale up and down for simulations using the fully-coupled models [57].

2.1.2. Loosely-Coupled Scheme

The surface and groundwater can be solved individually by using a loosely-coupled scheme, either in succession without iteration, or iteratively within each time step. Two or more separate models could be used in this approach. There are many loosely-coupled models that were developed in the past that depended on these techniques [25].

Governing Equations

In the subsurface, the governing equations are derived from Darcy’s law, with the continuity equation [58]. The mass conservation equation for a confined aquifer can be written as (Mehdinejadiani et al., 2013):

$$\frac{\partial }{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)=S_s \frac{\partial h}{\partial t}$$

(8)

where $h$ is piezometric head; $t$ is time$; K_x$, $K_y$, and $K_z$ are the principal components of the hydraulic conductivity tensor; and $S_s$ is specific storage.

For an unconfined aquifer [58]:

$$\frac{\partial }{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)=0$$

(9)

Solving the equations for unconfined groundwater flow is complicated due to the fluctuations of the aquifers’ thickness and variations, which could occur due to the withdrawal of groundwater [33,59]. The Dupuit principle and Boussinesq equation can be used to address this issue [60,61]. For loosely-coupled models, the coupling technique should be strong enough to transfer the data spatially and temporally [61]. Many of the loosely-coupled models have been established for specific catchments. It should be noted that some extra terms can be added to the governing equations, which should be solved by infiltration models. Infiltration models offer a systematic framework to understand the unsaturated zone and the groundwater flow. Many infiltration models had been established in the past. Green–Ampt, Horton’s, Philip’s, and Holtan’s are some of the well-known infiltration models.

A loosely-coupled scheme has been implemented with a direct linkage from hydrogeological data to hydrogeological numerical models in GIS, such as ArcArAz [55]. Unstructured meshes have been shown to be a unique characteristic in loosely-coupled models to simplify the boundary conditions. The efficiency of the loosely-coupled approach has been tested well in complex hydrological processes. The output of the surface water model has been applied as an imposed flux for groundwater recharge rates in MODFLOW [62]. It is also noted that some parameters, such as groundwater recharge rates, have the capacity to connect groundwater and surface water models. The importance of the transition zone between surface and groundwater flow has been noted from previous simulation results [62].

The loosely-coupled scheme has some advantages, such as the flexibility to apply individual tools to each environmental process, rather than one tool for all of the processes. This approach offers alternative software to reuse the hydrogeological conceptual model in GIS [55]. Similarly, the loosely-coupled scheme provides preference options in determining the platform to be used for any definite stage [63].

Several loosely-coupled frameworks have been developed for regional scales. However, there are some disadvantages, such as the difficulty in obtaining concrete decisions from the evaluation based on previous studies, because each structure was implemented for a particular catchment. The other disadvantages of the loosely-coupled scheme include the lack of a proper methodology to evaluate groundwater models. The groundwater model calibration is not sufficient for transient models with different recharge rates. The consideration of the transition zone is only possible by the construction of a fully-integrated dynamic model, which has two directional linkages between the surface and groundwater [62,64]. The geological properties of regional aquifers are normally simplified as heterogeneous in the estimation of spatial distribution of groundwater recharge. For this reason, the application of loosely-coupled models is extremely problematic, especially in the context of regional-scale modelling [65,66,67,68]. Convergence could be a challenging concern in the case of nonlinear models, including iteratively coupled models [7]. Inadequate statistics and information with respect to the interpolation of river water levels between gauges always creates difficulties when dealing with the measurement of the exchange term based on pressure differences [34]. For the estimation of exchange fluxes relying on pressure variances, the determination of river stages through the interpolation between the gauges is required, but loosely-coupled models are incapable of calculating river water levels at any point along the reach of a river. Apart from these above-mentioned complications, complexities in loosely-coupled models might also arise due to insufficient data of riverbed elevations with respect to a common datum (sea level) in potentially relevant stretches where the channel bottom would have to be derived using proxy data.

To date, MIKE SHE and FEFLOW, coupled with MIKE11 are the most commonly used loosely-coupled models for regional studies. Some loosely-coupled models are listed in

Table 2. Loosely-coupled software.

GW Model and Scheme	SW System	Unsaturated System	References
FEFLOW-3D Finite element	MIKE11-1D	2D; HELP	[69,70]
FHM-MODFLOW	HSPF-1D Stream Channel	Single storage reservoir-HSPF	[70]
MODBRANCH-MODFLOW	Saint-Venant-1D	Not Effective	[71]
IHMS-MODFLOW	DiCaSM	DiCaSM	[72,73]
GSFLOW-MODFLOW	PRMS-1D Stream Channel	1D Kinematic Wave to Richards’ equation	[74,75,76,77,78]
MIKE SHE-3D Finite Difference	Saint-Venant-2D	2-layer Water Balance Method	[79,80]
IWFM-Quasi 3D finite element/finite difference	Kinematic Wave Routing-1D	1D	[81,82]
MODCOU-SAM	Muskingum 1D/HECRAS-1D	Soil-water Budget	[83,84]
MODHMS-MODFLOW	Channel Flow Package (CHF1)-1D Saint Venant	Richards’ 3D	[85]
SWAT-MODFLOW	SWAT	Not mentioned	[86]
MODFLOW-MO-FLOW	MODFLOW stream package	MODFLOW stream package	[87]

.

2.2. Local Scale Modelling

Local scale modelling is applied for smaller spatial areas for surface and groundwater interaction. The cross section of the river can be incorporated in a particular small region with simplifications of other processes outside the research area. However, the influence of these outside processes could vary, and may play a significant role in determining the interaction. Thus, boundary conditions should be carefully defined, such as the groundwater recharge, runoff formation, and regional groundwater flow [34].

2.2.1. Traditional Methodologies

Several traditional techniques have been used in local scale modelling, including local grid refinement, local analytical correction, and local numerical correction [88]. In the following, these three techniques will be briefly reviewed.

Local grid refinement is one of the most common methods for well dynamics prediction. The resolution of the grid is dependent on the domain and knowledge of the local geology/hydrogeology [89,90,91]. Two or more distinctly sized grids have been associated with local grid refinement, in which a coarse grid occupying a huge area with the inclusion of regional boundary settings and one or more finer grids enclose the local areas. A one-way coupling methodology can be used to establish the link between the coarse and fine grids. In this approach, the conditions simulated by the coarse grid are applied into the fine boundaries, or through the process of two-way feedback between the two grids [92].

Table 3. Advantages and disadvantages of the local grid refinement technique.

Advantages	Disadvantages
Refinement or subdivisions of large grid cell into smaller spatial dimensions for the area of interest yields a more precise approximation of hydraulic head or drawdown on the well scale.	The refinement process is not well suited for large-scale, regional groundwater models because of remarkable increment in computational resources [93].
The solution could be obtained swiftly along with consistency between the regional and local zone around the wells, particularly for simple or small-scale.	Multiple sources and sinks, complex aquifer structure, transient flow conditions, multiple scaled of interest and strong anisotropy and heterogeneity associated with large problems would make this solution procedure very steep and might often lead to deficiency of convergence or numerical undulations [88].
The regular arrangement within the child grids and the resulting reduced computational burden due to large declination in the total number of cells can be avoided [93].	An uneven configuration across the parent-child boundary and the improvement and maintenance of separate model files are the main obstacles of this method [88].

highlights some of the advantages and disadvantages of this methodology.

The local analytical correction method was developed using a finite difference scheme. This technique is characterized by better enhancement for the forecasting of drawdown in the pumping well, which can be compared to the exact solution [94,95]. However, this corrected drawdown method would not be applicable in several realistic scenarios, such as for multiple sources and sinks within the well block with variable hydraulic conductivity, anisotropic media with detailed coordination, variable recharge, fractional well penetration, and transient flow or pumping circumstances [96,97,98,99,100].

The local numerical correction is based on the subdivisions of moderately large finite-difference grid cells from a regional model into multiple cells with gradually smaller spatial dimensions. This technique can be applied continuously until the preferred resolution of the well scale is achieved. It can be performed as an independent model after the creation of the local model. Restrictive assumptions of the local analytical correction and potential complications associated with local grid refinement can be avoided in this method. The conversion of facts and figures from a regional to a local model, and the choice of the initial and final boundary conditions for the generated local model have been identified as the most significant issues in this approach. As such, the use of multiple smaller-scale local models instead of resolving very large, complex matrices has been implemented under this process. This approach has the capability to minimize a huge and difficult matrix system into multiple smaller and improved conditioned matrices [88]. However, this approach could be very time consuming with the introduction of new boundary conditions or scales.

2.2.2. Hierarchical Modelling Approach

The hierarchical approach is an interactive and periodic technique, which is primarily based on the nested grid method [88]. It has the capability to address some of the disadvantages in the traditional methodologies by using dynamically coupled and fully integrated settings [94]. A dynamically integrated hierarchical patch dynamics paradigm (HPDP) has been established for groundwater modelling. This approach has some unique features, such as reduced assumptions, data input and post-processing, which significantly improve its calculating efficiency. The simulation results for regional model and groundwater systems are satisfactory without solving large matrix structures [67,90]. This approach also raises confidence in estimating exchange fluxes [100]. The dynamic fusion policy in the model has made real-time integration and subscale modelling possible. Computation steering application within the approach can be applied to manipulate the simulation process.

Furthermore, a computational steering application within this approach is able to provide the platform for human intelligence, allowing for human interruption during the time of execution, with modifications in the final output, such as the addition of patches in particular areas and layers of concern [66]. In short, the permission to focus indeterminately into a specific zone as long as the governing equations can maintain a certain scale is one of the vital advantages of this modelling strategy. The hierarchical approach is a realistic and time-efficient way to carry out complicated hydrological field settings. However, extensive data in a well-recognized structure is required for this process. The users may need a physical illustration of the database. Hierarchical structures have the tendency to adapt slowly to changing requirements. This system includes data to be frequently stored in many different entities. The limitations of this approach include deficiencies of structural independence, as well as a complex navigation scheme [101,102].

3. Domain Integrated Modelling

Domain integrated modelling is used for a simplification of the water cycle. Both the water quantity and water quality can be forecast. Two types of domain integrated models can be named: process-based models and stochastic models.

3.1. Process-Based Models

Process-based models can be referred to as deterministic hydrology models which incorporate all of the physical processes from the real world. This type of model uses streamflow, surface runoff, subsurface flow, etc. in individual-event models and continuous simulation models. Some of the highlights from previous studies using process-based models can be found in the following section.

Process-based models have been applied to represent the hydrological and ecological cycle in prairie wetlands [103]. The impacts of different water budget constituents were identified. Lateral contribution as runoff from the individual wetlands and their respective catchments were considered to be a fully integrated hydrological system. The prairie wetland region is hydraulically interconnected to the shallow aquifer systems, with continual flow exchange between precipitation, surface water and groundwater. The permanence of the wetland ponds has been proposed as another important criterion, which usually relies on factors such as runoff inputs and the hydrological position within a landscape, etc. The function of the artificial drainage of wetlands is critical to flooding, stream-water quality, and groundwater recharge.

A physically-based and distributed interaction model has been applied for greenhouse gas emission scenarios. The temporal dynamics of surface and groundwater interaction under climate change have been numerically investigated. From the simulation results, the impact of groundwater on stream flow is significant under climate change scenarios [104]. In order to improve the representation of surface and groundwater interaction, a groundwater model known as ‘CLASS’ has been developed. The interaction of surface and groundwater has been further established, including water table dynamics and land surface parameters. Other experimentally-based simulations, such as CK, CK-GW and CK-GWL have been included in the CLASS model in order to identify the impacts of regional hydrology variations. Among these, CK-GWL was significantly modified in the past few years in order to improve the accountability of subsurface lateral flow in addition to recharge and discharge processes. Many studies have used the CK-GWL model to simulate surface runoff, stream flows, wetter soil moistures, and a high water table [24]. However, the application and results obtained from statistical principles recommended further investigations using the regional climate model for better simulations of soil moisture and stream flows through the inclusion of the groundwater component in land surface modelling. Some other research has been conducted to identify the main reasons for groundwater fluctuations [105]. It was noted that groundwater level fluctuation is a combined consequence of many parameters, including direct and induced recharge, the extraction rate, and lateral groundwater flow, etc.

3.2. Stochastic Models

Stochastic models are developed like a black box by using mathematical and statistical theories to create a linkage between input data and the model’s output. Regression, neural networks, transfer functions, system identification are some of the methodologies that have been used in stochastic modelling. The application of these techniques is reviewed in the following.

Time series modelling is one of the most well-known strategies to estimate hydrogeological properties, e.g., flow and storage issues. This modelling approach has been applied due to its accuracy and cost-efficiency [106,107,108,109]. Multi-annual groundwater records can be used to describe aquifer properties. However, delayed gravity yield can generate erroneous groundwater hydrographs in time series modelling. The errors are evident in well pumping events at the initial time, which can be minimized by incorporating this effect into a transfer function model. Future research should be conducted to investigate the efficiency of this modelling approach [110]. The regression method is another technique that has been widely used in previous studies to investigate the groundwater table and groundwater management [111,112,113]. For instance, short-term groundwater table fluctuations have been generated using a triple linear regression method. The interrelationship among all of the key parameters in the determination of groundwater fluctuations can be established. It was found, from this regression and sensitivity analysis, that precipitation, evaporation, and river stages have the most significant role [114]. It is recommended that geographical locations and hydraulic conductivity should be considered in the determination of groundwater table variations using this method.

4. Discussion

The groundwater and surface water interaction is complex. A comprehensive framework is required in order to consider the climate, geological, and hydrological conditions. The interaction is critical for the successful simulation of groundwater research due to its direct influences on recharge–discharge approaches [63]. The interaction is desired to improve the modelling of groundwater and surface water. From the present review, it can be seen that the scale and uncertainty arise in the groundwater and surface water interactions. Temporal and spatial scales on both groundwater and surface water quantity and quality have not been addressed clearly [115]. Local and regional scale alterations and interchange have not been measured precisely. Therefore, the simulation results may vary significantly.

Various coupling techniques have been compared in the present review. Fully-coupled models can deliver an instantaneous clarification to groundwater and surface water interactions by combining process-based equations with 3D subsurface demonstrations. The fully-coupled model works best when the groundwater and surface water interactions are needed at a larger or regional scale, with the consideration of all of the hydrological processes. Fully-coupled models are most effective compared to other schemes in integrated systems. This approach also has the capacity to handle complex physics in regional scale modelling. However, the lack of field data in large-scale modelling has been a major hindrance for the successful application of fully-coupled models [116]. A large number of computational resources is needed for fully-coupled model simulation, due to the large number of estimated parameters. The calibration could be complicated, due to the high uncertainty in the boundary and initial conditions. Thus, a loosely-coupled approach has been used in the past with relatively simple coupling schemes.

Many loosely-coupled structures have been developed for local-scale modelling, especially for regional studies, such as MIKE SHE, FEFLOW coupled with MIKE 11. The loosely-coupled models can communicate one-way information to sub-models with mean coupled equations and feedbacks [117]. In this approach, the selection and integration of sub-models, such as groundwater, surface water and the unsaturated zone, results in the uncertainty of the modeled results due to variables in the input parameters within each of the sub-models. Some of the models are dependent on surface water models, while others are not. Additionally, the boundary placement between the coupling models is another challenge in loosely-coupled models.

5. Conclusions

Groundwater and surface water models can be integrated, but the terminology for the model coupling has yet to be defined well from the literature. Diverse coupling approaches need to be compared for future modelling applications. Fully-coupled models have the capabilities to resolve the Richards’ and shallow water equations simultaneously, and thereby can provide proper explanations of any complex physics scenario, as well as estimating the parameters for regional scale models. However, the large number of required parameters and a limited knowledge of the boundary conditions may increase the complexity. On the other hand, loosely-coupled methodologies based on empirical relationships can be used in a mixture of process-based models, or process-based models and other algorithms. The results from single models can be utilized as inputs for other models. Simplification in calibration and uncertainty evaluations are the main benefits of loosely-coupled models. Convergence difficulties for nonlinear models as well as disparities among the individual models during the time of their independent calibrations are some of the drawbacks of loosely-coupled models. Examples of software that were developed based on both fully-coupled and loosely-coupled models have been listed, which might help researchers choose the appropriate models when interaction flux needs to be considered. The models could be further classified by the boundaries at the interface between the groundwater, the vadose zone, and the surface water. The assumptions of all of the models could be simplified by reducing some of these uncertainties as well as the variable spatial and temporal features of the groundwater and surface water flows. A separate suitable framework is desired in order to address the uncertainties for future development.