Modelling of shallow water table dynamics using conceptual and physically based integrated surface water-groundwater hydrologic models

We present a new conceptual scheme of the interaction between unsaturated and saturated zones of the MOBIDIC (Modello Bilancio Idrologico DIstributo e Continuo) hydrological model which is applicable to shallow water table conditions. First, MODFLOW was coupled to MOBIDIC as the physically based alternative to the conceptual groundwater component of the MOBIDIC-MODFLOW. FirstThen, assuming a hydrostatic equilibrium moisture profile in the unsaturated zone, a 10 dynamic water table level dependent specific yield was added to MOBIDIC-MODFLOW and calculation of the groundwater recharge in MOBIDIC was revisited using a power type equation based on infiltration rate, soil moisture deficit, and a calibration parameter linked to the initial water table depth, soil type, and rainfall intensity. a hydrostatic equilibrium moisture profile was assumed for simulating changes in water table levels. This resulted in a water table based expression of specific yield, which was included in the coupled 15 MOBIDIC-MODFLOW modelling framework for capturing shallow water tables fluctuations. Second, the groundwater recharge was defined using a power type equation based on infiltration rate, soil moisture deficit and a calibration parameter linked to initial water table level, soil type and rainfall intensity. Using the Water Table Fluctuation (WTF) method for a homogeneous soil column, the water table rise for a homogeneous soil column under a pulse of rain with different intensities (up to 30 mm/day) the parameter of the proposed groundwater recharge equation was determined for four 20 soil types i.e., sand, loamy sand, sandy loam and loam under a pulse of rain with different intensities. The fidelity of the introduced modifications in MOBIDIC-MODFLOW was assessed by comparison of the simulated water tables against those of MIKE SHEThe simulated water table levels were compared against those simulated by MIKE-SHE, a physically based integrated hydrological modelling system simulating surface and groundwater flow, in two numerical experiments. Two numerical experiments were carried out: a two-dimensional case of a hypothetical watershed in a vertical plane (constant slope) 25 under a 1cm/day uniform rainfall rate, and a quasi-real three dimensional watershed under one month of measured daily rainfall hyetograph. The comparative analysis confirmed that the simplified approach can mimic simple and complex groundwater systems with an acceptable level of accuracy. In addition, the computational efficiency of the proposed approach (MIKE -SHE took 180 times longer to solve the 3D case than the MOBIDIC-MODFLOW framework) demonstrates its applicability to real catchment case studies. 30

physically based integrated hydrologic models, the description of the flow in the unsaturated zone in externally coupled models is not based on the Richards equation and their simplified unsaturated-saturated coupling scheme can restrict their applicability in modelling of shallow water table fluctuations.
The challenging issue regarding the applicability of these externally coupled models in humid shallow water table regions is 25 related to the inconsistencies in the conceptualization of the interaction between unsaturated and saturated zones. Seibert et al. (2003) distinguishes three types of interactions between unsaturated and saturated zones based on the water table levels: Type 1. The water table is relatively deep and there is only one-directional interaction between UZ and SZ. This means the soil moisture state in the unsaturated soil is independent of the groundwater level and the role of groundwater in runoff generation process is not considered. 30 Type 2. The water table is about at the root zone level. The unsaturated soil can get water from the capillary rise in groundwater so the interaction becomes two-directional as groundwater recharge can be either positive or negative to create hydrostatic equilibrium with the water table. The unsaturated soil profile is assumed to have a constant vertical extension (unsaturated moisture capacity (the maximum amount of moisture the unsaturated layer can hold) is constant during the course of simulations). Such assumptions were made by different hydrological models such as TOPMODEL (Beven et al., 1995), SWAT (Neitsch et al., 2011) and MOBIDIC (Castillo et al., 2015).
Type 3. The water table is very close to the surface and unsaturated moisture capacity can no longer be assumed to be constant as the water table fluctuates. This is the case when water table rise results in a decrease in the unsaturated soil moisture capacity 5 and a small amount of the infiltration causes a significant rise in groundwater level (due to the significantly small moisture deficit in unsaturated zone). This is the main thrust of this paper as we introduce a conceptual approach compatible with the existing MOBIDIC framework for applications where strong interactions between the UZ and SZ exist.
The main problems in modelling of type 3 shallow water tables using externally coupled integrated models lies in 1) the assumption of constant specific yield in the saturated flow module of the integrated model is incompatible with the nonlinear 10 decrease in its value as water table rises up to the soil surface; and 2) the assumption of constant vertical extension for the unsaturated zone does not reflect the inverse relation between the unsaturated and saturated moisture storages discussed in (Duffy, 2013;Seibert et al., 2003).
With regard to these limitations, the objective of this study is to propose a series of modifications to the original conceptualization of the unsaturated and saturated flow process of MOBIDIC, in order to extend its applicability for modelling 15 of shallow water table fluctuations while retaining its computational efficiency. To this aim, the conceptual saturated flow scheme of MOBIDIC was replaced with MODFLOW as a physically based three-dimensional groundwater model using the sequential coupling approach (Guzha and Hardy, 2010). Then, a novel methodology for revisiting the calculation of the groundwater recharge in MOBIDIC, the specific yield in MODFLOW, and the interaction between the unsaturated and saturated zones in MOBIDIC-MODFLOW was developed. The fully coupled surface-subsurface model MIKE SHE is used as 20 a reference for comparison; hence the methodology is based on numerical benchmarking on hypothetical and realistic catchments. Using the Water Table Fluctuation (WTF) method (Healy and Cook, 2002) in MIKE SHE, the rises of a shallow water table were simulated under different sets of rainfall intensity, soil property and depth to water table. The simulated responses were then used to reformulate the groundwater recharge of MOBIDIC based on the assumption of a quasi-steady pressure profile in the unsaturated zone as water table fluctuates. The accuracy of the proposed modifications was first 25 evaluated in a two-dimensional case (constant slope), where the simulated water table rises of the two models under a uniform rainfall rate were compared. In a second experiment, the approach was tested at the catchment scale and under unsteady rainfall conditions. Comparison of the simulated water table responses of the MOBIDIC-MODFLOW against those of MIKE SHE allowed us to evaluate how the unsaturated-saturated interaction scheme of the externally coupled models can be adapted for applications in shallow water table regions. 30 More specifically, the calculation of groundwater recharge is revisited as the existing linear function in MOBIDIC (groundwater recharge is a linear function of moisture storage in gravity reservoir) is incompatible with the nonlinear behaviour of moisture fluxes between UZ and SZ typical of shallow water table regions. Such nonlinear behavior is associated with nonlinear decrease in magnitude of specific yields as the shallow water table rises to the soil surface.
In this paper, we applied the water table fluctuation method (WTF) (Healy and Cook, 2002) to investigate the behavior of groundwater recharge in shallow water tables (up to 1.5m deep) for four soil types i.e., sand, loamy sand, sandy loam and loam in a soil column under a single pulse of rainfall with different intensities. These simulations were carried out using MIKE-SHE as a reference model. MIKE-SHE assumes a constant specific yield (in saturated flow computation process) and the computed water table is corrected based on a specified mass balance threshold error in UZ-SZ coupling process. The 5 simulated water table levels of MIKE-SHE are then considered as the 'true' hydrologic response of the system and are utilized to reformulate the groundwater recharge model component of MOBIDIC. The accuracy of the proposed changes is first tested in a two-dimensional case where subsurface water is simulated in a vertical plane with constant slope. A constant rainfall rate is applied and the rise in groundwater levels is affected by groundwater recharge and by the lateral interaction between the saturated computational grids. In a second numerical experiment, the accuracy of the approach is further evaluated at the 10 catchment scale and under unsteady rainfall where the simulated water table levels of the two models (MIKE-SHE as the reference model and MOBIDIC-MODFLOW) are compared.

Water table fluctuation method
The Water Table Fluctuation (WTF) method is a simplified approach for the determination of groundwater recharge of an unconfined aquifer based on groundwater level fluctuations. This method is based on the assumption that the rise in 15 groundwater levels is due to the groundwater recharge (Healy and Cook, 2002). Considering the groundwater budget for a representative element, (Fig. 1), any change in the water table level (groundwater storage) would be due to a combination of recharge to groundwater ( ), inflow from upstream cell ( ), outflow to the downstream cell ( ) and evapotranspiration from groundwater ( ) as follows: is the change in groundwater storage. Evapotranspiration from groundwater is taken as the direct root water uptake from the saturated zone . Unlike in deep water table conditions, groundwater evapotranspiration in shallow water   table regions can be much greater than evapotranspiration from the unsaturated zone (Shah et al., 2007).
Assuming that water table rise is solely due to the recharge of groundwater requires the sum of other fluxes in Eq.1 to be zero.
This means that the determination of the groundwater recharge using WTF is best applicable over short periods (hours to days) 25 after onset of rainfall (before any significant redistribution of groundwater recharge to the other fluxes) (Healy and Cook, 2002). Therefore: Where [−] is the specific yield and Δℎ Δ is the change in water table over Δ . Application of the WTF method requires an estimation of the specific yield which is defined as the volume of drained water per unit drop in water table and aquifer area 30 (Nachabe, 2002): Where [ 3 ] is the volume of drained water, [ 2 ] is the area of the aquifer and ∆ℎ[ ] is the change in water table level. The specific yield is also defined as the difference between water contents at saturation and at field capacity (the moisture level below which water cannot be drained by gravity (Nachabe, 2002)). Such constant value for specific yield, however, holds 5 only for deep water tables where changes of the soil moisture profile in the unsaturated zone due to water table drop are relatively small and the volume of drained water can be approximated as ( − ) × ∆ℎ (Fig. 2a). In shallow water tables however, the specific yield is small (in Fig.  2b, shaded area is smaller than ( − ) × Δℎ ) and approaches zero when the capillary fringe zone extends up to the soil surface (Nachabe, 2002).

MIKE SHE
MIKE SHE is one of the widely used physically based integrated surface-subsurface hydrological model for a wide range of spatiotemporal applications ranging from detailed theoretical (single soil column) to operational watershed scale studies (Graham and Butts, 2005). It has a modular structure for computation of the hydrological processes with different levels of complexity, which is advantageous particularly in large-scale watershed studies (Kollet et al., 2017). A detailed description of 15 the computation of the hydrological processes in MIKE SHE can be found in Storm (1991) andDHI (2014). We present only the computation of flow in the unsaturated and saturated zones and their coupling approach.

Unsatuarted flow
The unsaturated flow is described using the one dimensional Richards' equation as (Downer and Ogden, 2004):

Saturated flow
Saturated flow in MIKE -SHE is computed using the three-dimensional saturated flow equation as: Where , , [ −1 ] are the saturated hydraulic conductivity along the , and axes, respectively, ℎ[ ] is the groundwater head, [ −1 ] is the source/sink term and is the storativity coefficient [−]. The storativity coefficient is either 5 the specific yield, , for unconfined aquifers or the specific storage , , for confined ones (DHI, 2014). Equation (5) is solved numerically over computational grid squares using the finite difference method with the Preconditioned-Conjugate Gradient (PCG) solver, which is also included in MODFLOW (DHI, 2014).

Unsaturated-Saturated zone coupling
The explicit coupling approach implemented in MIKE SHE has the advantage of employing different times steps for each 10 zones (seconds to minutes for UZ and hours to day for SZ), which makes the system computationally less expensive compared to the implicit coupling approach (DHI, 2014). However, employing different time steps for the UZ and SZ may result in the generation of mass balance errors in calculating the water flux between the two zones due to 1) an incorrect value of (recommended value is − ) in computing the saturated flow and 2) a fixed water table level while the UZ progresses to the next time step (Storm, 1991). The coupling between UZ and SZ follows an iterative process by which the water table is 15 adjusted until the accumulated mass balance error for the entire soil layer (UZ+SZ) falls below a prescribed threshold. The process also calls for adjusting the soil moisture profile. The groundwater recharge is calculated once the iterative process has converged. It should be noted that is considered as constant in the computational process. The algorithm of iterative adjustments of the water table can be summarized as follows (Storm, 1991): Step 1) at the beginning of the UZ time step, the accumulated mass balance error for the entire soil layer (UZ+SZ) of a soil 20 column ( ) is calculated.
Step 2) if the accumulated error falls below the acceptable level, the water table does not require any adjustments and groundwater recharge is calculated as follows: Where is groundwater recharge, ∑ is the net inflow to the unsaturated soil. The calculated recharge is applied as the upper 25 boundary condition of the SZ module and simulation advances to the next time step.
Step3) if the calculated error is beyond the acceptable threshold, then the water table adjustment is required. Depending on the sign of the water table will be either raised or lowered.
3-a) if < 0, it means less water is stored in the profile and therefore water table has to be raised. Using the updated moisture profile in UZ (only the last three nodes of the UZ profile (above water table) are updated), go to 30 step 1 and repeat.

3-b) if
> 0, it means more water is stored in the profile and therefore the water table is lowered. Using the updated moisture profile in UZ (same as described above only for the three lowest nodes of the UZ profile), go to step 1 and repeat.
The iterative water table adjustment continues until the calculated falls below the acceptable limit. Then the new groundwater recharge is calculated as: 5 Where ℎ * and ℎ are the water table after and prior the adjustments, respectively. The calculated groundwater recharge will then be applied as the upper boundary condition of the SZ module and simulation advances to the next time step.

MOBIDIC
Modello Bilancio Idrologico DIstributo e Continuo (MOBIDIC) (Castelli et al.,2009) is a distributed continuous hydrologic 10 model in which the components of the hydrologic system are conceptualized as a system of inter-connected reservoirs. Such a conceptual formulation of the model makes it computationally more efficient especially in large scale watershed modeling (Castillo et al., 2015). The detailed description of the model is available in Castelli et al. (2009) and Castillo et al. (2015).

Unsaturated flow
In MOBIDIC, the unsaturated zone is described by two interconnected reservoirs i.e., the gravity and capillary reservoirs, to 15 account for corresponding gravity and capillary forces in the unsaturated soil. The water content at field capacity ( ) is the threshold below which moisture is entirely held in the capillary reservoir. The gravity reservoir interacts with the saturated zone via percolation to the groundwater and it may also redistribute the additional moisture to the downstream cells (lateral flow in the unsaturated zone) (Castillo et al., 2015). The capillary reservoir serves moisture to the plant roots (or evaporation if there is no plant in the computational grid) and can also take moisture from lower groundwater via capillary rise (Castillo et 20 al., 2015). The moisture capacities of the gravity reservoir and the capillary reservoir, [L] and [L], are defined as (Castillo et al., 2015): where [ ] is the unsaturated soil thickness and [−] are the water content at saturation, field capacity 25 and residual, respectively, which are determined based on soil texture classification (Rawls et al., 1982). At each time step, water available for infiltration (net precipitation plus the ponded water on the soil surface) is determined. The infiltration rate takes the minimum value among the saturated hydraulic conductivity, the allowable moisture capacity in gravity reservoir and the available water. The infiltration rate, however, might be underestimated for low permeable soils in dry conditions when the capillary force at the soil surface is not negligible Castelli (1996). The gravity 9 reservoir is replenished by the infiltrated water. Absorption flux (moisture that is extracted by from the gravity reservoir to the capillary reservoir) is calculated as (Castillo et al., 2015): Where [ −1 ](0 ≤ ≤ 1) is a linear coefficient that controls the rate of the moisture transfer between the two reservoirs.
Finer soils typically have higher value because of low downward moisture gradient. The groundwater recharge ( ) is 5 calculated as (Castillo et al., 2015) : Where [ ] is the depth to water table, , [ ] is updated moisture in gravity reservoir and [ −1 ](0 ≤ ≤ 1) is a coefficient that controls the rate of groundwater recharge. The available moisture storage in the gravity reservoir can also 10 contribute to lateral flux to the adjacent cell and is calculated as (Castillo et al., 2015): Where [ −1 ](0 ≤ ≤ 1) is a coefficient that determines the rate of lateral flow and , is the updated moisture state in gravity reservoir after reduction of groundwater recharge. It should be noted that, as the unsaturated flow is strictly vertical in MIKE SHE, we ignored the lateral redistribution of moisture in gravity reservoir ( =0) to make the structure of the two models 15 consistent.

Saturated flow
Saturated flow in MOBIDIC is described either by a simplified linear reservoir (conceptual scheme) or the Dupuit assumption (Dupuit, 1848). In this study, or by MODFLOW (Harbaugh et al.,2000) was coupled to MOBIDIC as a physically based threedimensional finite-difference based alternative groundwater flow model.(Three-dimensional Finite-Difference based 20 groundwater model). We used the latter to facilitate comparison between the simulated results of MOBIDIC and MIKE-SHE regarding their respective conceptualization of the saturated flow and its effect on the interaction between UZ and SZ. MOBIDIC and MODFLOW were coupled using the sequential coupling approach in which the groundwater recharge and water table depth act as the boundary conditions of the coupled MOBIDIC-MODFLOW model (Guzha, 2008). At each daily time step, the water table level determined from the solution of MODFLOW in the previous time step is used in MOBIDIC 25 for the calculation of groundwater recharge (equation 12). The calculated groundwater recharge is then applied as the upper boundary condition in MODFLOW and the water table level of the computation grid is updated. Note that, unlike in MIKE SHE, no fine spatiotemporal discretization of the system is required in MOBIDIC-MODFLOW, in which a daily time step has been used.

Coupling Unsaturated and Saturated zones
Unlike MIKE -SHE, coupling UZ-SZ in MOBIDIC-MODFLOW is not based on an iterative water table correction procedure.
Based on the calculated water table level in the previous time step, the recharge to groundwater can be positive (recharge to groundwater) or negative (extraction from groundwater) Castillo (2014). The latter occurs when the saturated storage is bigger than moisture storage in the gravity reservoir ( ). Therefore, the water table in the subsequent time step falls to establish a 5 hydrostatic equilibrium with moisture level in gravity reservoir Castillo (2014). However in the former, groundwater is recharged by a higher moisture level in the gravity reservoir. This results in a rise in water table in the subsequent time step Castillo (2014). During the dry periods, the capillary reservoir may also receive water from the capillary rise from water table which is calculated as (Castillo et al., 2015): where [ ] is the mean distance of the unsaturated layer to water table, 1 [ ] is the bubbling pressure, [−] is the product of pore size distribution index (Rawls et al., 1982) and pore size disconnected index (Brooks and Corey, 1964). The unsaturated soil water pressure ( ) is a function of saturation state of the layer and pore size distribution index ( ) and is calculated as (Castillo et al., 2015): It should be noted that the capillary rise is computed where the water table is within the soil profile ( ≤ ), otherwise the capillary rise flux will be zero.
In Table 1, the differences between MIKE SHE and MOBIDIC-MODFLOW with regards to the conceptualization of the unsaturated flow, saturated flow, UZ-SZ coupling process, and limitations associated with the application in very shallow 20 water tables are summarized. Unlike their disparity in describing the unsaturated zone, the two models share similar formulation and numerical solution technique for the saturated flow module. The iterative water table adjustments approach in MIKE SHE, takes into account the variations of specific yield as water table rises/falls. In MOBIDIC-MODFLOW, the specific yield is instead defined as a soil hydraulic parameter and remains unchanged during the course of simulation. This, however, leads to an underestimation of the water table rise as a result of a very small value of specific yield as the water table 25 approaches the ground surface (Abdul and Gillham, 1989).

Water Table Fluctuation method for a soil column using MIKE -SHE
We use the WTF method over for a soil column to understand how the rise in groundwater is affected by rainfall intensity, soil types and depth to water table level and how much of the infiltrated rain will percolate to the groundwater. The step-by-step procedure is shown in Fig.4. Using a single soil column with closed boundaries on the sides and bottom, the rise in groundwater level will be only due to the groundwater recharge (see eqs. 1 and 2). In shallow water table regions, as the water table rises, the specific yield decreases nonlinearlynonlinearly decreases and, therefore, the actual rise in water table might be greater than what would be expected by assuming = − . This is defined as ''reference water table rise'' (see Fig. 3) and is calculated as: where [L/T] is the infiltrated water and ∆ℎ[ ] is the reference water table rise which is the expected rise in water table if the specific yield remains at its ultimate value ( − ) and the infiltrated water completely reaches to the water table.
However, depending on the initial water table level and rainfall intensity, it is possible for the water table rise to be smaller than ∆ℎ (e.g. for relatively deeper water table and/or low rainfall rate) or that the groundwater recharge takes a value greater 10 than the infiltration rate. The actual water table rise was determined through the coupled unsaturated-saturated flow simulation in MIKE SHE and knowledge of ultimate specific yield allowed the calculation of groundwater recharge using equation 17.
Therefore, by defining ( ) we can evaluate how the rise of the water table is related to the precipitation rate (assuming no infiltration excess runoff), soil hydraulic property, and depth to the water table. >1, =, < 1 means that the actual rise of the water table is larger, equal, or smaller than the reference rise, respectively. Once again, in MIKE SHE, is a constant input 15 of the saturated flow module and changes regarding its value (decrease in its value when depth to water table decreases) is captured by the iterative unsaturated-saturated coupling procedure explained in section 3.1.3. The procedure was repeated for four soil types (see Table 2 for the hydraulic properties) and different rainfall intensities (up to the minimum of saturated hydraulic conductivity of the four soils to prevent infiltration excess runoff generation on the soil surface) to further investigate the significance of the rise in the water table of soil types with different characteristics in response to a given pulse of rain. In 20 shallow water tables, loamy soils with larger capillary fringe range have substantially smaller specific yield values which result in a higher rise of water table compared to the sandy soils Gillham (1984). Therefore, its actual rise in the water table is expected to be larger than the reference rise. However, as loamy soils have smaller hydraulic conductivities compared to sand, their water table response will be delayed especially for deep water tables. To avoid this delayed response, we focused our analysis to cases where the initial depths to the water table are at a maximum of 1.5m. Computational time steps were one 25 second and one minute for UZ and SZ, respectively to avoid numerical instabilities in the simulated water tables by MIKE SHE.
its value approaches zero when capillary fringe is extended up to the soil surface. This means that, in such conditions a small amount of rain may result in a significant rise in water table (up to the soil surface if the water table is close enough to the surface such that the capillary fringe reaches the surface) as investigated by Gillham, 1984, 1989;Buttle and Sami, 30 1992;Sklash and Farvolden, 1979;Waswa and Lorentz, 2015). The schematic illustration of the method is given in Fig. 5. If the infiltrated water is completely transferred to the groundwater, the rise in water table will be: Where [L/T] is the infiltrated water and ∆ℎ[L] is the rise in water table level. ∆ℎ is referred as the 'reference rise'. However, depending on the initial water table level and rainfall intensity, it is possible that the water table rise be smaller than ∆ℎ (e.g. for relatively deeper water table and/or low rainfall rate) or that the groundwater recharge takes a value greater than the infiltration rate. The 'actual rise' of the water table then can be calculated same as in Eq.2. Therefore, by defining ( ) we can 5 evaluate how the rise of the water table is related to the precipitation rate (assuming no infiltration excess runoff)and depth to the water table. If >1, =, < 1 means that the actual rise of the water table is larger, equal, or smaller than the reference rise, respectively. Once again, in MIKESHE, the is a constant input of the saturated module and the changes regarding its value (decrease in its value when depth to water table decreases) is being captured by the coupling procedure explained in section 3.1.3. 10 The procedure was repeated for different soil types i.e., sand, loamy sand, sandy loam and loam. The hydraulic properties of the soils are given in Table 1. This was done to further investigate the significance of the rise in the water table of soil types with different characteristics in response to a given pulse of rain. In shallow water tables, loamy soils with larger capillary fringe range have substantially smaller specific yield values which result in a higher rise of water table compared to the sandy soils Gillham (1984). Therefore, its actual rise in the water table is expected to be larger than the reference rise. However, as 15 loamy soils have smaller hydraulic conductivities compared to sand, their water table response will be delayed especially for deep water tables. To avoid this delayed response, we focused our analysis to cases where the initial depths to the water table are at a maximum of 1.5m. Computational time steps were one second and one minute for UZ and SZ, respectively to avoid numerical instabilities in the simulated water tables by MIKESHE.

Simulation results 20
Simulation results of the water table rises are shown in Fig. 45. Each dotted curve in the plots is associated with a specific initial depth to water table (ranging from 0.3 m to 1.5 m) and rainfall rate. The plots in Fig.4 5 are divided into two zones i.e., >1 (water table rise exceeds the reference rise, red dots) and <1 (water table rise is less than the reference rise, blue dots).
Whereas the lower bound of the chosen initial depth to water table in Fig. 4 5 is 1.5m, the upper bound changes from 0.3 m for sand to 0.9 m for loamy soil (see Table 1). This is due to the differences in water retention characteristics of the soils and 25 initial moisture deficit in unsaturated zone. For example, in sandy soil with 30 cm initial depth to water table (capillary fringe is 15.98 cm), the moisture deficit is 14.98 mm. So, for the infiltration rate bigger than this value, the unsaturated soil becomes completely saturated and water table rises up to the soil surface as it can be seen in Fig. 4 5 in which( the upper curve ends at a rainfall rate of about 15 mm/day). The same argument applies to the other soil types. For loamy soil with capillary fringe value of 40.12 cm, the initial depth to water table is of 90cm corresponds to the moisture deficit of 23.97 mm in unsaturated 30 zone. Therefore, the upper curve for loamy soil ends at this depth to water table and rainfall intensity.
For a given initial water table level, moving from low to high rainfall intensities results in increasing the ratio of recharge to infiltration, shifting from a situation where <1 (blue dots) to >1 (red dots). Also, while the case where <1 is more common in sandy soils, the opposite is found as soil texture becomes finer (e.g. loamy soil). This means that the actual water table rise in fine textured soils (loamy sand to loam) is almost always larger than would be expected based on using a constant value for specific yield = ( − ). 5 Note that for all soil types analyzed, the > 1 values show small sensitivity to increases in rainfall intensity as the depth to water table gets larger, e.g. in sandy soils and loamy soils when the water table is deeper than 50cm and 1 m, respectively.
This can be attributed to the existence of significant initial soil moisture deficits found at larger initial water table depths. On the other hand, as the water table gets closer to the soil surface, an increase in the rainfall rate generates a significant water table rise, yielding to a substantial increase in . This demonstrates the importance of the soil moisture deficit (and its ratio 10 to the precipitated rainfall) in significance of water table rise. magnitude and dynamics.

Changes in conceptualization of the UZ-SZ interactions in MOBIDIC-MODFLOW
As discussed in section 3, MOBIDIC's unsaturated soil depth ( in Eqs. (7) and (8)) remains constant during the course of a simulation (both in rainy and subsequent draining periods). This is due to its conceptualization of UZ and SZ as their interaction is not aimed to address the reverse relationship between them in very shallow water table cases. The proposed changes in this 15 paper are intendedaim to make the model applicable for such cases. The changes are as follows: 1-It is assumed that the unsaturated soil layer thickness, , is no longer a constant input of the model and changes with water table fluctuations. Hence, the total moisture capacity of the reservoirs ( and ) are determined similarly as in Eqs.
(7) and (8), but with replacing by the depth to water table ( ). Therefore, water table rise/fall results in a decrease/increase in moisture capacities of the reservoirs. Such assumption is valid when groundwater level is treated as a moving boundary and 20 there is a continuous transfer of moisture between the unsaturated and saturated zones.
2-The specific yield in calculation of groundwater head by MODFLOW, is determined based on a soil water retention model (e.g. Brooks and Corey, 1964) and hydrostatic equilibrium assumption in unsaturated zone (suction profile in unsaturated zone changes from steady state to another over the changes in water table) (Hilberts et al., 2005) Duke (1972: 25 Where is depth to water table and other variables are as defined above. Validity of hydrostatic equilibrium assumption is justified in shallow water table regions where redistribution of the infiltrated water is rapid and an equilibrium moisture profile occur immediately Bierkens (1998). Therefore, with any changes in water table level, specific yield is adjusted for the calculation of groundwater dynamics in the next time step. This is especially important in shallow water table cases as specific yield decreases significantly as water table gets closerrises up to the soil surface.
3-The groundwater recharge [ −1 ] is a function of available moisture in the gravity reservoir, unsaturated soil moisture deficit and infiltration rate and is calculated as: Where [ −1 ] is the infiltration rate, and [−] is a function of soil type, infiltration rate and water table level, which determines how much water will be transferred to the saturated zone. It should be noted that, similar to the original 10 conceptualization of the unsaturated flow in MOBIDIC, only gravity reservoir contributes to the groundwater recharge. The is an increase in these quantities. It should be noted that, unlike MIKESHE, the unsaturated and saturated modules in MOBIDIC-MODFLOW run with a daily time step which means any change in the water table level causes an instantaneous adjustment of unsaturated moisture profile.

Interpretation of the new groundwater recharge equation in MOBIDIC-MODFLOW
To interpret the proposed groundwater recharge equation (Eqs. (19) and (20)), let us first analyze its limit values. 20 As 1+( ) < 1, a decrease in results in a decrease in (and vice versa). The coefficient can be either positive or negative. When positive, recharge to the groundwater reservoir is larger than the infiltration rate and when negative, groundwater recharge is a portion of the infiltration. The specific case for which = 0 occurs when R = I. Such flexibility in value range of is important as it makes the model applicable for cases where adding small rainfall amount can cause a significant water table rise, which may occur in humid regions. However, its precise determination can be problematic as it is 25 a complex function of water table level, soil type, and rainfall (infiltration) intensity. In general, for a known water table level, an increase in rainfall rate will require a larger values to match the water   relationship is nonlinear. This further shows the complexity of the interactions between the UZ and SZ under shallow water table conditions. Note that in sandy soils, most of the values are negative, while the opposite is found for the fine textured soils such as loamy soil. This is because the water table rise for loamy soils are larger than for sandy soils (smaller specific yield and available moisture deficit in unsaturated profile). Comparing the scatterplots in Fig. 67, we realizedreveals that that unlike loamy soils, in sandy soils with relatively deeper initial water tables, the changes in (with changes in rainfall rates) is 10 quite small. This is related to the initial moisture deficit in unsaturated profile as discussed in section 4.1.

Application of the modified unsaturated-saturated scheme of MOBIDIC-MODFLOW scheme in two test cases
The appropriateness of the proposed changes in the conceptualization of UZ-SZ interactions of MOBIDIC-MODFLOWin shallow water tables are was tested by comparing its predicted simulated water table results against those of MIKE SHE in two test cases. The description of the test cases follows.

Two-dimensional case 25
So far, our analysis focused on a single soil column in which the coefficient of the proposed groundwater recharge, (Fig.   67), were determined based on the simulated water table rise by MIKE SHE (Fig. 45). In order to analyze the performance of the proposed approach in cases where there is lateral interaction between saturated computational grids, simulations of the water table rise were performed for a small soil box measuring 14m long, 1.5m high and 1m wide filled by sandy soil identical to the one used in the soil column analysis. The soil box has closed boundary conditions on = 0, = 14m and bottom side and it is assumed that initial water table is ℎ( , 0) = 0.2m. Since the rainfall is continuous during the entire simulation period, the configuration of the box such boundary conditions results in higher expected water table rise in the toe of the slope (shallower water table) and generation of saturation excess runoff starting from downhill moving to the uphill. The hillslope is discretized into 28 columns ( = 0.5m) and the vertical discretization of the unsaturated layer in MIKE SHE is = 1cm.
A 10 mm/day constant rainfall is applied for 20 days and the simulated water table rises of the two models are compared. 5

Simulation results
The simulated water tables of the two models are shown in Fig. 78. It can be seen that the predicted water tables of the two models are very closeThe two models in general show similar behavior in predicted water tables. The slight mismatch between the predicted water table heads of the two models can be attributed to the fact that in MOBIDIC-MODFLOW, there is an immediate response of the the response of groundwater to the precipitated water is immediate, that is, within the same time 10 step (one day). This is not the case in MIKE SHE where the simulation time step is much shorter at one second. The saturated length (the length where water table at the soil surface) predicted by the two models are closely matched which shows that the. This means simplifications in UZ-SZ interaction of MOBIDIC-MODFLOW can mimic the complex dynamical interaction between the two zones. Note that theThe generated saturation excess runoff by two models were removed from the soil surface as the flow routing module was not included in the simulations. 15

Borden catchment
In order to further evaluate the suitability of the proposed conceptualization of UZ-SZ in MOBIDIC-MODFLOW, the water table fluctuations of the two models at the Borden catchment for the month of May 2015 were compared. The Borden catchment is located approximately 70 km of Northwest of Toronto, Ontario (Jones et al., 2006). The site is about 20m wide 20 and 80m long and it was subjected to several experimental (Abdul and Gillham, 1989) and numerical rainfall-runoff studies (Jones et al., 2006;Kollet et al., 2017;VanderKwaak and Loague, 2001). It is assumed that the catchment consists of a single homogeneous sandy soil identical as the one used in the previous simulations (see Table 12 for hydraulic properties) without any vegetation cover. The motivation for simplifying the watershed physiographic characteristics here was to evaluate the MOBIDIC-MODFLOW in a 'real' watershed, while emphasizing on differences in UZ-SZ dynamics simulated by the 25 conceptual approach as compared to a physically based numerical model. The digital elevation model (DEM) of the Borden catchment has a spatial resolution of 0.5m and is available at www.hpscterrsys.de and is shown in Fig. 8a9a. The initial water table level was assumed to be at 1m below the catchment's outlet (at elevation 1.98m). The site has closed boundaries on the sides and it has a horizontal bedrock at elevation 0 m. The measured rainfall at the Borden site for May 2015 used in the numerical experiment (data available at Environment Canada, 2015). The 30 rainfall hyetograph (Fig. 1011) consists of 10 events with total amount of 77.2 mm and maximum of 24.8 mm at day 25 where a significant rise in water table is expected. The catchment was subdivided into two zones based on the initial depth to water table as shown in Fig. 8b9b. Zone 1 represents the initial depth to water table less than 1.5m (red zone in Fig. 8b9b) and cells in zone 2 have the initial depth to water table larger than 1.5m (blue zone in Fig. 8b9b). The calculation of groundwater recharge and interaction between UZ-SZ followed the modifications presented in section 5 (the parameter in Eq. (19) was made variable for water table levels less than 1.5m) and changes in groundwater head in zone 2 was simulated using the 5 original conceptualization of MOBIDIC explained in section 3.2. Therefore, these two structures together cover the dynamics of the water table for the whole catchment. Note that, for water table depths greater than 1.5m, the moisture capacities of the reservoirs ( , ) remains constant with a rise or fall of the water table (section 3.2). Such distinction based on depth to the water table was made to limit the extrapolation errors associated with the determination of for water table depths larger than 1.5m. Besides, for these zones the assumption of an inverse relationship between the unsaturated and saturated 10 moisture storages is questionable as discussed in the introduction.

Simulation results
The difference in simulated water table level of the two models is shown in Fig. 9 10 for the selected simulation days after the rainy events. The predicted groundwater head of the models at the outlet of the catchment are also compared in Fig. 1011. It can be seen that the two models generally compare well. 15 However MOBIDIC-MODFLOW is slightly predicting a higher water table rise as compared to MIKE SHE following rainfall events (e.g. events on days 11, 18 and 25). This is attributed to the interaction of saturated flow between the two zones (the saturated flow moves from grids in zone 2 to zone 1). The magnitude of the overestimation decreases over the course of simulations (for example the differences in water table levels are smaller at day 23 than day 6). This further shows that how the interaction between the zones with different characteristic of the interaction between the UZ and SZ (zone 1 with inverse 20 relation between the unsaturated and saturated storage and zone 2 with direct relationship between them as classified in the introduction) can affect the overall behavior of the water table fluctuations at catchment scale. Therefore, that the rise in water table in a computational grid is also the result of the net incoming of saturated flow from surrounding grids for both models.
That is why the predicted water able in MOBIDIC-MODFLOW is still slightly rising between two rainy days (from day 12 to day 17 and day 19 to day 24) as shown in Fig. 1011. Consequently, the simulated water tables of the two models gradually 25 converge in the days following a rainfall event (for example between the days 14 to 17 and 19 to 24) as the transient behavior simulated by MIKE SHE dissipates following a long redistribution period.

Discussion
The quick rise of the shallow water table in response to the precipitation was also observed in the experimental work of (Abdul and Gillham 1984). In their study, the water table response of a sandy soil packed in a 160 ×120 × 8 cm box with a 12° surface 30 slope with different initial water table levels under a uniform rainfall rate was investigated. The objective of the experiment was to evaluate the effect of the capillary fringe on the rise in water table and on streamflow generation. The results revealed models at the catchment scale can be cumbersome due to the required data and computational burden to run such models. On the other hand, the conceptual shallow water table unsaturated-saturated interaction schemes in the literature e.g., (Seibert et al., 2003) do not allow for capturing of the different water table responses of different soil types to the fallen pulses of precipitation shown in Figure 5. The coefficients of such conceptual models are often determined through the calibration process against few water table observations which does not represent the disparity in water Investigation of such cases using conceptual hydrologic models have not received much attention in the literature (e.g. Seibert et al., 2003). (Seibert et al., 2003)  It should be noted that, While while the  coefficient parameter of the proposed groundwater recharge in Eq.(19) was determined for sand, loamy sand, sandy loam and loam, the model was tested for sandy soils where the moisture retention in unsaturated zone is small and the assumption of quasi steady equilibrium moisture profile in unsaturated zone is legitimate, where it was assumed that the unsaturated soil instantaneously adjusts to changes in water table level as discussed by (Bierkens, 1998). For finer textured soils, however, such the assumption of instantaneous equilibrium state between unsaturated and saturated zones may be questionable and, therefore,. The approach the approach presented in this study here should thereforemust be fully tested over a variety of soil textures to further confirm evaluate its exactitude to broader watershed conditions. The proposed changes were intended to keep the MOBIDIC's number of calibration parameters low which is advantageous in model calibration as it does not raise the concern on equifinality (Beven, 2001). The coupling approach 5 implemented in MOBIDIC-MODFLOW was tested in two-dimensional (single slope) and three dimensional cases (Borden catchment), both consisting of a homogeneous sandy soil. In addition to testing the proposed approach over a variety of soil textures, a next step will consist in extending the applicability of the approach against observations in real catchments where other hydrological processes, such as summer-fall evapotranspiration and spring snowmelt, affect groundwater recharge in shallow aquifers. 10

Conclusion
In this paper, we have presented a series of modifications to the conceptualization of the unsaturated and -saturated zones inflows of MOBIDIC a conceptual hydrologic model for application to extend its applications in to shallow water tables in which there iswhere an inverse relationship between the unsaturated and saturated moisture storages exist. The modifications were: 1) replacement of the model's conceptual saturated flow scheme with MODFLOW; 2) revisiting the conceptualization 15 of the interaction between the unsaturated and saturated zones and calculation of groundwater recharge; and 3) inclusion of a dynamic specific yield based on quasi-steady state pressure profile in the unsaturated zone. The key variables in modelling of such casesshallow water table regions are 1) groundwater recharge and 2) specific yield. Groundwater recharge is usually considered as a function of unsaturated moisture state and is a fraction of infiltrated moisture into the unsaturated soil layer.
However in shallow water table cases,s its value can be larger than infiltration due to the capillary fringe above the water table, 20 which causes a quick and significant water table rise (Jayatilaka and Gillham,1996). Specific yield in the shallow water table regions is no longer a constant soil property and significantly decreases as water table reaches to the soil surface. Investigation of such cases are mostly addressed using physical based models (e.g. Cloke et al., 2006) or using fully conceptual approaches (e.g. Seibert et al., 2003). Through the coupling of MODFLOW to MOBIDIC, We presented MOBIDIC-MODFLOW as a conceptual-numericalwe have modified the unsaturated-saturated flow interaction of the model based on the assumption of 25 quasi-steady pressure profile in the unsaturated zone as water table changes. In turn, the static specific yield in MODFLOW could be replaced with the water table dependent dynamic specific yield such as derived by Duke (1972). Such modifications allowed the rapid and dynamic water table responses of shallow water table regions to be adequately captured. model in which changes in water table are assumed to establish an equilibrium moisture profile in the unsaturated soil zone. Employing the analytical expression of specific yield derived by Duke (1972) in the coupled model, its dynamical behavior particularly in 30 shallow water table regions was captured.Also, The the groundwater recharge in the model was defined with a power type equation whose parameter ( ) was determined using the Water under variable rainfall intensityover a month of measured daily rainfall hyetograph. In the light of simulation results, the 5 following conclusions can be derived: 1-Inclusion of a dynamic specific yield in investigation of water table behavior in shallow water tableof shallow water table responses is importantessential. Indeed, when the water table is close to the soil surface, the significance of rise is much greater than would be expected based on the ultimate value of specific yield e.g., − , as investigated by Jayatilaka and Gillham (1996). While in MIKE SHE, this is accounted for by using the iterative water table 10 correction method, in MOBIDIC-MODFLOW this is addressed by assuming a hydrostatic equilibrium state between UZ-SZ. To the authors' knowledge, this has not been taken into account in simplified coupling of a hydrologic and a groundwater model.

2-
The presented analysis showed a complex relation between the water table rise and rainfall rate, soil type, and depth to water tables.The analysis of the water table rise with different rainfall intensity, soil type and water table level 15 shows a complex relation between the groundwater recharge and these factors. The simulation over different ranges of rainfall rates and water table levels showed that fine textured soils with large capillary fringe (consequently small unsaturated soil moisture deficit) could have water table at soil surface even with a low rainfall intensity. This is important for the investigation of saturation excess runoff in lowland near river zones of the catchments. Therefore, it is important for the model to distinguish such differences in the expected water table responses of different soil 20 types for its suitability at the catchment scale.
3-The comparison of the simulated water tables of the two models at Borden case (Mean Absolute Error equals 1.3 cm) along with their execution times (30 hours with MIKE SHE against 10 minutes with MOBIDIC-MODFLOW) clearly demonstrates the applicability of simplified approach implemented in MOBIDIC-MODFLOW in investigation of the groundwater-surface water interaction. Further improvements of the approach with inclusion of other hydrological 25 processes such as the effect of evapotranspiration on parameterization of the () in future work will enable its applicability in real, practical applications.