Hillslope-storage Boussinesq model for simulating subsurface water storage dynamics in scantily-gauged catchments

Accepted Manuscript

Hillslope-storage Boussinesq model for simulating subsurface water storage dynamics in scantily-gauged catchments Soumyaranjan Sahoo , Bhabagrahi Sahoo , Sudhindra N. Panda PII: DOI: Reference:

S0309-1708(17)31020-5 https://doi.org/10.1016/j.advwatres.2018.08.016 ADWR 3189

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

6 November 2017 20 August 2018 26 August 2018

Please cite this article as: Soumyaranjan Sahoo , Bhabagrahi Sahoo , Sudhindra N. Panda , Hillslope-storage Boussinesq model for simulating subsurface water storage dynamics in scantily-gauged catchments, Advances in Water Resources (2018), doi: https://doi.org/10.1016/j.advwatres.2018.08.016

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights: 

We enhanced the hsB model considering surface ponding, saturation-unsaturation zone dynamics, and bedrock leakage flux.



We applied the enhanced hsB model in a semi-gauged real headwater catchment.



We advocated that all the hydrological fluxes must be considered for its field

AC

CE

PT

ED

M

AN US

CR IP T

application.

1

ACCEPTED MANUSCRIPT

Hillslope-storage Boussinesq model for simulating subsurface water storage dynamics in scantily-gauged catchments

Soumyaranjan Sahooa*, Bhabagrahi Sahoob, and Sudhindra N. Pandac

Research Scholar, School of Water Resources, Indian Institute of Technology Kharagpur,

Kharagpur-721302, India, E-mail: [email protected] b

CR IP T

a

Associate Professor, School of Water Resources, Indian Institute of Technology Kharagpur,

c

AN US

Kharagpur-721302, India, Tel: +91 3222 281884 (O), E-mail: [email protected]

Professor, Agricultural and Food Engineering Department, Indian Institute of Technology

Kharagpur, Kharagpur-721302, India, E-mail: [email protected]

ED

M

*Corresponding author: Soumyaranjan Sahoo ([email protected])

Abstract

PT

Limited geo-hydro-meteorological data availability in most of the subsurface water catchments world-wide has constrained the effective surface and subsurface water management. For

CE

improved subsurface water modelling in real semi-gauged catchments, the direct-rainfall-

AC

recharge (DRR)-based low-dimensional hillslope-storage Boussinesq (hsB) model is enhanced by incorporating saturated phase processes of surface ponding (SP) in agricultural lands, unsaturated zone processes during non-monsoon (dry) period, bedrock leakage (BL) flux between the hillslope-aquifer and underlying confined aquifer, and root-zone water balance (WB) on the basis of precipitation forcing, catchment topography, and land-use/land-cover dynamics. The developed approaches are field-tested by considering the catchment as both semi-gauged regarding limited water table elevation (WTE) data and fully ungauged (no WTE 2

ACCEPTED MANUSCRIPT data). The results reveal that the developed WB-SP-BL-based hsB model variant performed the best in reproducing the observed WTE in six experimental wells in the Kanjhari reservoir catchment (125.04 km2) in eastern India with the Nash-Sutcliffe model efficiency > 67%, root mean square error

0.28 m, and percentage bias of ±2.5% for both the gauged and ungauged

scenarios. Conclusively, it is endorsed that accounting for all the hydrological fluxes in the hsB

semi-gauged catchments.

CR IP T

framework would make it amenable for useful prediction of subsurface flow in ungauged and

KEYWORDS: groundwater; recharge; hillslope drainage; hillslope-storage Boussinesq model;

AC

CE

PT

ED

M

AN US

subsurface storage-discharge; ungauged catchments

3

ACCEPTED MANUSCRIPT 1. Introduction Assessment of subsurface storage and baseflow is essential for sustenance of riverine hyporheic zones (Merill and Tonjes, 2014), aquatic ecosystem by providing thermal (Kurylyk et al., 2014) and chemical (Broda et al., 2014) refuge, planning of domestic and agricultural water supply systems (Matonse and Kroll, 2009), and estimation of environmental flow

CR IP T

(Ouyang, 2012). To study the subsurface flow phenomenon at hillslope-scale, the existing methods are based on the Dupuit-Forchheimer (D-F) approach (Dupuit, 1863), Boussinesq equation (Boussinesq, 1877) or the numerical solution of complex three-dimensional (3-D) equations (e.g., Troch et al. 1993; Chen et al., 1994a, 1994b). However, due to the assumptions

AN US

of 1-D flow-domain, hillslope of unit-width with homogeneous soil, free-surface boundary overlying the saturated strata, and simple end-point boundary conditions, many of these theories cannot be considered realistic (Paniconi et al., 2003). Similarly, field application of 3-

M

D Richards equation is not always feasible in ungauged or scantily-gauged catchments due to the difficulty in model parameterization, unavailability of information on soil hydraulic

ED

properties at the required spatial scale, inaccuracy and numerical instability problems

PT

associated with coarser spatiotemporal model discretization, and high computational time (Paniconi et al., 2003; List and Radu, 2016). Consequently, in the hydrologic literature, for

CE

modeling the subsurface flow dynamics, an emphasis has been given formulating simple, lowdimensional and physically-based models accounting for the spatiotemporal variability in

AC

atmospheric, geomorphologic, and pedo-hydrologic features of the catchment (Duffy, 1996; Woods and Sivapalan, 1999; Grayson and Blöschl, 2000). Fan and Bras (1998) laid the first attempt to collapse a 3-D soil mantle to a 1-D soil

moisture storage capacity profile, which is simple, yet effectively 3-D. In the hillslope-storage Boussinesq (hsB) model, Troch et al. (2003) improved the approach of Fan and Bras (1998) and extended the Boussinesq equation by incorporating the continuity and Darcy’s equations.

4

ACCEPTED MANUSCRIPT The hsB model can serve as a potential tool for ungauged catchment modeling (Troch et al., 2004) as it is of reduced dimensionality with parsimonious model parameters and is simpler than the recent 3-D hillslope-based models (e.g., Hazenberg et al., 2015). However, in spite of the improvements in the theoretical framework and modelling domain of the hsB model and its variants, it is evident from the very few real-world application cases that proper estimation of

CR IP T

the component hydrological fluxes is necessary for predictions in ungauged/scantily-gauged catchments (e.g., Matonse and Kroll, 2009; Carrillo et al., 2011; Broda et al., 2012; Troch et al., 2013) . Furthermore, only two of them (e.g., Carrillo et al., 2011; Troch et al., 2013) consider some of the essential hydrologic processes of infiltration and evapotranspiration,

AN US

neglecting the vadose zone moisture transport dynamics. Moreover, the previous real-world applications of this model are restricted only in the well-studied, highly instrumented and experimental data-sufficient hillslopes (e.g., Matonse and Kroll, 2009; Carrillo et al., 2011;

M

Matonse and Kroll, 2013; Liu et al., 2016).

While simulating low streamflows by the hsB model, Matonse and Kroll (2009) considered

ED

recharge input into the hsB model by subtracting evapotranspiration and interception losses

PT

from precipitation, however, neglecting infiltration dynamics. This approach, henceforth termed as the direct-rainfall-recharge (DRR) approach, ignores the Hortonian runoff

CE

mechanism; and overestimates the recharge that can invite substantial error in predicting low streamflows constraining the model’s applicability to handle only the saturation excess runoff.

AC

For baseflow simulation, Matonse and Kroll (2013) improved the approach of Matonse and Kroll (2009) by coupling the hsB model with the lumped Sacramento Soil Moisture Accounting (SAC-SMA) rainfall-runoff model, requiring a number of a priori parameters to incorporate the components of evapotranspiration, land cover, interflow, baseflow, overland flow, and lateral flow. Due to the underlying complexity, this approach may not be suitable for scantily-gauged catchments.

5

ACCEPTED MANUSCRIPT The problem is further aggravated in case of the monsoon-type tropical catchments, in which the topsoil layer is subjected to precipitation, irrigation, and evapotranspiration. The water, percolated from the topsoil, first passes through the transmission zone based on its soilmoisture dependent unsaturated hydraulic conductivity, and finally reaches the water table as recharge. This phenomenon is observed during the non-monsoon (dry) period (see Fig. 1b).

CR IP T

Conversely, in many tropical monsoon-type catchments, paddy crop is the dominant land use during the monsoon season, wherein, to deal with the dry-spells through adequate water storage, the croplands are bounded with dikes (earthen embankments) of about 10–40 cm in height. These local controls store precipitation and surface runoff resulting in ponding and

AN US

surface saturation and form seasonal artificial wetlands (Chen and Liu, 2002). Note that infiltration during ponding condition is always higher than that during no ponding condition due to the additional head of ponded water (Chen et al., 2002; He et al., 2009; Jha et al., 2017).

M

When the percolation rate from the topsoil is higher than the hydraulic conductivity at any soil layer in the transmission zone, the upper layer of this zone gets saturated due to the

ED

addition of water. During this situation, an unsaturated zone is sandwiched between two saturated zones (Fig. 1c). The lower one is the saturated Boussinesq aquifer, where the flow is

PT

lateral according to Dupuit’s assumption. Conversely, in the top saturated zone, flow is vertical of the

CE

due to gravity and soil suction is developed due to the negative soil matric potential,

underlying vadose zone. The flow continues to be vertically downward as long as the negative

AC

soil matric potential is present. With continuous addition of moisture or water, the unsaturated moisture content,

changes to saturated moisture content,

moisture, positive pressure (

; and with the further addition of

) builds up. Hence, the above dynamics creates a saturated zone

with a moving (variable) lower boundary above the transmission zone. In this case, the root zone water balance approach (e.g., Carrillo et al., 2011) may not yield accurate results.

6

ACCEPTED MANUSCRIPT Furthermore, the bedrock leakage (BL) flux is a dominant component of the hillslope-scale overall water balance (Graham et al., 2010) that influences the subsurface stormflow response (Tromp-van Meerveld et al., 2007), improves the simulation of pore water pressure distribution (Ebel et al., 2008), and is essential for simulating the subsurface flow response at the hillslopeor catchment-scale (Meerveld and Weiler, 2008).

CR IP T

Therefore, the specific research gaps in the previous hsB-based studies are identified as: i) Non-consideration of the fluxes to and from the aquifer, ii) Assumption of impermeable bedrock, iii) Lack of hsB parameterization method for its application in ungauged catchments, and iv) Assumption of spatially uniform aquifer hydraulic properties. These limitations

AN US

constrain the applicability of the model in scantily-gauged catchments and necessitate for the development of a generalized approach accounting for all the hydrologic fluxes to the hillslope aquifer.

M

In light of the above viewpoints, this paper tries to answer the research questions: 1) How including more hydrological processes to estimate recharge could improve the classical DRR

ED

approach in the hsB-based model?; and 2) How to parameterize a low-dimensional model in an

PT

ungauged catchment? To address these issues, this study tries to advocate a novel recharge-flux estimation method embedded in the hsB model framework by considering all the physical

CE

hydrological processes in the vadose zone overlying the water table. As recharge varies with precipitation, topography, and land use and land cover (LULC) characteristic of the basin,

AC

these dynamics are accounted for in the enhanced hsB model framework incorporating the water balance modeling of topsoil along with the inclusion of unsaturated zone processes during non-monsoon (dry) period. Similarly, during monsoon (wet) period, recharge due to surface ponding (SP) in croplands with rainwater is estimated by considering the saturated phase modeling. Validating the low-dimensional hsB-based model in a data-scarce catchment would open avenues for understanding the subsurface dynamics in ungauged catchments since

7

ACCEPTED MANUSCRIPT in the hydrologic literature, an emphasis has been given for ‘prediction in ungauged basins’ (Sivapalan et al., 2003; Hrachowitz et al., 2013). Hence, the present approach augments the hsB model to make it a wholesome package for real-catchment modeling. The manuscript is structured as follows. Section 2 briefs the hsB model and the earlier efforts to advance the theoretical and model framework; and Section 3 details the model

CR IP T

enhancements brought out in this study. Study area and application of the enhanced model variants in the study area are described in Sections 4 and 5, respectively. Section 6 discusses

2. The hsB model: a brief description 2.1. A brief description

AN US

the results, and Section 7 concludes the study.

Estimation of subsurface storage-discharge along a hillslope unconfined aquifer of unit-width,

M

resting on sloping bedrock was given by Boussinesq (Childs, 1971). Combining it with the

domain is expressed as (

)

+

(1)

PT

*

ED

continuity equation, the modified form of the Boussinesq equation accounting for 1-D flow

where h = groundwater elevation perpendicular to the bedrock [L]; t = time [T]; Ks = saturated

CE

hydraulic conductivity [LT-1]; f = drainable porosity [-]; α = bedrock slope angle of underlying

AC

impermeable layer; x = flow distance from the catchment outlet [L]; and N = rainfall recharge rate [LT-1].

Fan and Bras (1998) introduced the hillslope width function (HWF) [L] to represent the 3-D

soil mantle into 1-D form by a soil moisture storage capacity function, Sc(x) = fw(x)D, where D(x) is the average soil depth at a flow distance x from the outlet. Incorporating HWF in the continuity form of Eq. (1) and combining with Darcy’s equation, Troch et al. (2003) presented the 1-D hsB model to study the subsurface storage dynamics, expressed as 8

ACCEPTED MANUSCRIPT *

(

)+

(2)

Equation (2) is applicable in hillslopes with shallow soils on a relatively impermeable bedrock (or lower permeable subsoil), where the D-F assumptions are valid. It computes groundwater stage hydrograph, h(x, t) at any distance x from the channel. The subsurface discharge hydrograph, Qgw(x=0, t) at the outlet can be obtained by a simple mass balance of the

CR IP T

saturated storage as (∫

∫

)

(3)

In Eq. (3), the first and second terms represent the input flux and aquifer saturated storage,

AN US

respectively. 2.2. Assumptions

In this study, the following assumptions undertaken for developing the hillslope-scale hsB model are also valid:

The D-F assumptions for groundwater flow hold.

ii.

The hillslope-storage varies along the longitudinal distance x from the outlet; however,

ED

M

i.

this variability is ignored in the transverse direction of x (i.e., in the y-direction).

PT

2.3. Model history

Since its inception, the hsB model has undergone several improvements by inclusion of

CE

complex hillslope geometry and curved bedrock (Hilberts et al., 2004), non-constant drainable

AC

porosity by intervening drainable porosity as a function of water table height and bedrock depth (Hilberts et al., 2005), unsaturated zone modeling by modifying the flow domain (Hilberts et al., 2007), bedrock leakage (Broda et al., 2011), and lateral flow through the unsaturated zone (Kong et al., 2016). Recently, the hsB-based models and other lowdimensional models (e.g., Hazenberg et al., 2015, 2016; Pan et al., 2015; Sahoo et al., 2017) are also coupled with the surface flow models. Furthermore, Broda et al. (2012) developed the

9

ACCEPTED MANUSCRIPT hsB analytic element model (hsB-AE) to improve the hsB model dealing with shallow subsurface water by coupling it with the analytical element (AE) approach (Haitjema, 1995).

3. Enhancement of the hsB model 3.1. The model domain

CR IP T

To incorporate the dominant catchment processes in the enhanced hsB model framework, the hillslopes can be considered as the basic unit of study, in which, the entire soil column resting above the bedrock forms the hillslope unconfined aquifer (Fig. 1a). The aquifer can be assumed as a three-layered system (Figs. 1b and 1c), in which the depth of the first layer was selected as

AN US

45 cm herein which is based on the average root zone depth of the dominant vegetation in the study area. The second soil layer is the transmission (vadose) zone that extends up to the zone of saturation or unconfined aquifer water table, representing the third layer. These three layers

M

from the top surface exchange the fluxes of infiltration and evapotranspiration, surface- and subsurface runoffs, and bedrock leakage with the boundaries of atmosphere, stream network,

AC

CE

PT

ED

and underlying confined aquifer, respectively.

10

ACCEPTED MANUSCRIPT P

ET

Tf C

b) Non-ponding period

c)

P

Ponding period

CR IP T

Nl

P

ETpot

ETadj

Tf

Tf

Drz

RWB Transmission zone Dtz

C

Dynamic water table

NWB

D

C

h

Db

Nl

Upper saturated zone

RSP Transmission zone

Dynamic water table

NSP

Unconfined aquifer

M

Unconfined aquifer Bedrock

Ponding water Root zone

PD

Root zone

AN US

I

Unconfined aquifer

Bedrock Nl

ED

Fig. 1. (a) Conceptual flow domain of a hillslope-unconfined aquifer overlying the bedrock

PT

with slope angle α; and (b, c) Schematic of hydrological processes considered along with water

CE

pressure distributions in the respective zones in the enhanced hsB framework.

3.2. Hydrological process components

AC

Since the existing hsB model variants do not explicitly account for all the hydrological process components as mentioned earlier, this model was enhanced by i) Considering the surface and unsaturated zone water balances for recharge estimation during the non-monsoon (dry) period; ii) Accounting for the surface ponding and vadose zone processes during the monsoon (wet) period; and iii) Revising the original assumption of impermeable bedrock by a leaky bedrock, which is more realistic.

11

ACCEPTED MANUSCRIPT Recharge to the hillslope unconfined-aquifer water table depends on the unsaturated zone soil moisture convection dynamics through the transmission zone, which can be modeled using the Richards equation as *

(

)+

(4)

= volumetric soil moisture content in the transmission zone [L3L-3]; Kv = vertical

where

CR IP T

hydraulic conductivity [LT-1]; ph = variable soil water pressure head [L]; z = elevation above the bedrock [L]; and Kv(ph) can be estimated by the Mualem-van Genuchten equations (Mualem, 1976; van Genuchten, 1980; and Ines and Mohanty, 2008) ⁄

( *

[

]

) +

+

where Se = relative saturation [-];

(6)

are the residual and saturated moisture contents

= shape parameter [L-1], expressed as inverse of bubbling pressure; n

are the shape parameters [-] that account for pore size distribution and tortuosity,

respectively;

ED

and

(5)

M

[L3L-3], respectively;

and

AN US

*

= 0.5 for most of the soils (Mualem, 1976); and

⁄ . Equation (4)

PT

was discretized by finite difference scheme in time and space (soil column depth), with a layer spacing of 10 cm.

CE

Similarly, the capillary rise flux [LT-1] can be estimated as (Gardner, 1958) ]

(7)

AC

[

where

= coefficient [-] that is a function of b; b is the parameter related to the Brooks-Corey

soil water retention parameter = 2+3B (Eagleson, 1978); B = Brooks-Corey pore size distribution index (Brooks and Corey, 1964);

= soil matric potential at effective saturation

[L]; Dtz(t) = dynamic depth of the transmission zone [L] at any time t, expressed as Dtz(t) = [Dh(t-1)]/cosα (Fig. 1c); and F = capillary fringe depth [L]. The Mualem-van Genuchten model was used to parameterize the Richards equation for routing soil moisture in the vadose zone. 12

ACCEPTED MANUSCRIPT However, the Brooks-Corey relationships were used to model the daily-scale capillary rise only, without any routing. 3.2.1. Inclusion of surface water balance during non-monsoon (dry) period (WB approach) The hillslope unconfined aquifer interacts with the overlying root zone and the intermediate transmission zone. Hence, the net recharge flux to the hillslope aquifer depends upon the flux

CR IP T

between these layers, which can be estimated by the water balance approach in the root zone as (8)

where Drz = depth of the root zone considered for water balance approach [L];

AN US

volumetric soil moisture content of the root zone [-] at Drz = 45 cm; I(x,t) = infiltration rate at the land surface [LT-1]; ETa(x,t) = actual evapotranspiration from the root zone specific to any land use class [LT-1]; and RWB(x,t) = percolation rate from the surface soil layer to the transmission zone [LT-1].

M

The catchment-scale ETa was estimated for two scenarios. For the ponding phase, when

ED

there is no moisture stress, ETa is at the rate of the potential evapotranspiration (ETpot) estimated as (Swain and Sahoo, 2015) ∑ ∑

)

(9)

PT

(

CE

where CCi = crop coefficient of the th LULC class [-]; ALi = area of the th LULC class contributing to the total area [L2]; nl = total number of LULC classes constituting the

AC

catchment; and ET0 = reference evapotranspiration [LT-1], estimated by the FAO-24 Radiation model (Doorenbos and Pruitt, 1977). For non-ponding phase, when there is moisture stress, and the root zone depletion is

higher than the readily available water, the adjusted ETa (ETadj) can be estimated as (10) where MS = moisture stress coefficient, estimated as (Allen et al., 1998)

13

ACCEPTED MANUSCRIPT

(11) where AWT = total available soil moisture in the root zone [L]; AWR = readily available soil moisture in the root zone under the unsaturated condition, expressed as a fraction of the total available soil moisture that vegetation can extract without any water stress [L]; and RZD = root zone depletion [L].

CR IP T

The infiltration rate can be obtained as

(12)

where Ic = variable infiltration capacity of the surface soil [LT-1]; and Tf = throughfall rate [LT], estimated as (Carrillo et al., 2011) { where P = precipitation rate [LT-1];

AN US

1

(13)

variable canopy storage [L], estimated by canopy canopy

M

water balance accounting for the precipitation and evaporation;

storage capacity [L] (Dickinson, 1984); and LAI = leaf area index [-]. The variable infiltration

ED

capacity of the soil can be estimated by Horton’s infiltration model as (14)

PT

where If = steady state (ultimate) infiltration capacity [LT-1]; I0 = initial infiltration capacity

CE

[LT-1]; and H = Horton’s decay coefficient [T-1]. Although many other infiltration models, viz., the two-parameter Philip’s, Kostiakov and Green-Ampt models are available in the literature,

AC

estimation of the Horton’s model parameters in the field condition is easier due to simplicity; and also there are accuracy issues with the Kostiakov model at large time scales (Philip, 1957; Parlange and Haverkamp, 1989; Kale and Sahoo, 2011). Rearranging Eq. (8), one obtains the percolation rate into the transmission zone using the water balance (WB) approach as (15)

14

ACCEPTED MANUSCRIPT To estimate RWB, Eq. (15) is solved numerically first for the current value of θrz, and then for the updated values of θrz at subsequent time steps as illustrated in Fig. (2). Now, the transmission zone water balance can be solved to obtain the recharge from the ith LULC class as [

]

where z = elevation above the bedrock [L]. Incorporating Dtz, Eq. (16) can be re-written as [

(∫

)

(16)

CR IP T

(∫

]

)

(17)

AN US

Note that the depth of the transmission zone is dynamic due to the changing water table. Under such a scenario, the recharge can be higher than the downward flux of water as the rising water table takes up water from the transmission zone (Hilberts et al., 2007). Now, the enhanced hsB model accounting for the WB approach of effective recharge

(

)+

(18)

PT

*

ED

M

estimation during dry period can be given by

where NWB is a function of the LULC.

CE

3.2.2. Consideration of surface hydrological processes with surface ponding during monsoon (wet) period (WB-SP approach)

AC

Combining the continuity and Darcy’s equations, the recharge during ponding (saturated) phase can be estimated by the Laplace equation as (Agrawal et al., 2004) (19) where

hydraulic head [L], expressed as

; and pz = piezometric head [L].

Equation (19) was solved numerically using the central difference scheme by dividing the root zone depth into different sub-layers, each of 10 cm thickness. Applying the finite difference 15

ACCEPTED MANUSCRIPT form of Eq. (19) at all the grid points yields ng–1 number of simultaneous equations (ng = number of grids). The solution of Eq. (19) involves the dynamic saturated zone as the lower boundary; and the hydraulic head at the land surface as the upper boundary with the root zone node,

, expressed as (20)

CR IP T

where PD = ponding water depth in paddy fields [L] that occurs when the moisture content of topsoil θts>θs. PD depends on rainfall, irrigation, surface runoff, and evapotranspiration at any time t, which can be computed by the generic mass balance equation as (Khepar et al., 2000;

AN US

Agrawal et al., 2004)

(21)

where Qr = surface runoff [L]; RSP = percolation from the root zone under ponded condition [L]; and θts = moisture content in the topsoil [L3L-3].

)

(22)

ED

(

M

The ponding water state prevails if

If on any day, the PD value exceeds the maximum allowable value of 40 cm (the maximum

PT

dike height), then PD is reset to 40 cm, and the exceeded amount is assigned as Qr. The simultaneous equations were solved by the Gauss elimination method in matrix form to

CE

obtain Ψp at each grid point. Substituting Ψp in Darcy’s equation, the percolation flux from the

AC

root zone with surface ponding (SP) can be obtained as *

where

+

(23)

= vertical grid spacing [L] of the problem space = 10 cm; and Ψp(nrz) = hydraulic

head at the discretization node in the root zone. Hence, recharge under surface ponding condition for the ith LULC class can be expressed as

16

ACCEPTED MANUSCRIPT [

(∫

]

)

(24)

Now, the enhanced hsB model for the wet period considering the WB-SP approach can be expressed as

(

[

)+

3.2.3. Consideration of bedrock leakage flux (BL approach)

]

(25)

CR IP T

*

In the earlier hsB-based modeling studies, the oversimplification of impermeable bedrock could add error in the hsB model solutions, specifically when applied to a real-world

AN US

catchment. The bedrock leakage flux can be estimated using Darcy’s equation as (Broda et al., 2011) (

)

(26)

M

where Nl = bedrock leakage flux [LT-1]; Kb = bedrock vertical saturated hydraulic conductivity [LT-1]; Db = bedrock thickness measured perpendicular to the bedrock slope [L]; and ht and hb

ED

are the hydraulic heads [L] at the top and bottom of the bedrock, respectively. In this study, ht and hb are the water table elevations obtained from the two piezometers drilled up to the top

PT

and bottom of the bedrock, respectively.

CE

3.2.4. Generalized versions of the enhanced hsB model for alternate wetting and drying phases and land uses

AC

During the non-monsoon (dry) period, the enhanced WB-BL approach-based hsB model considering the water balance and bedrock leakage flux can be expressed for any land use class other than agriculture as * [

(

)+ ]

(27)

17

ACCEPTED MANUSCRIPT Similarly, for the agricultural land use, where substantial ponding can be expected during the monsoon (wet) season, the enhanced WB-SP-BL-based hsB model can be expressed as *

(

)+

[

]

(28)

CR IP T

Note that NWB and NSP are mutually exclusive in the sense that only one of these can be greater than zero. The presence of these two terms in Eq. (28) makes the expression more generic for addressing the problem of alternate wetting and drying (AWD) as there could be a period of non-ponding phase in between two ponding phases.

AN US

3.2.5. Model discretization and boundary conditions

Equations (18), (27) and (28) were discretized by the finite difference scheme in space and a multistep ordinary differential equation (ODE45) solver in time with ∆x = 30 m and ∆t = 1 h, with the output print time of ∆t = 1 day. Similarly, although the depth of the transmission zone

M

is varying, the spatialization of the recharge sub-models, NWB (Eq. (17)) and NSP (Eq. (24)) was

ED

fixed as ∆z = 0.1 m and ∆x = 30 m along the vertical and lateral coordinates. For solving the hsB model in a real-hillslope, the following boundary conditions were used: i) no-flux

PT

boundary at the catchment-divide (x = L, Qgw = 0) ii) Dirichlet boundary condition at the downstream channel outlet (x = 0, h = 0); and a recharge boundary, NWB-BL(x, t) or NWB-SP-BL (x,

CE

t) at the top surface depending on the hsB model variant to be used. The recharge forcings were

AC

computed as separate functions in the MATLAB environment as detailed in Fig. 2. The assumptions of x = 0, h = 0 is unrealistic for estimating the baseflow in a hillslope

where a seepage face is usually expected, and this creates a problem for calculating discharge at the lower boundary using Darcy’s law. However, this assumption is used by several researchers (e.g., Beven, 1981; Matonse and Kroll, 2009) to decouple it from the effect of a stream, as incorporating the stream-aquifer interaction is beyond the scope of the present study. Furthermore, since a numerical technique was used to estimate the discharge ensuring mass 18

ACCEPTED MANUSCRIPT balance, the simulated flow rate could be marginally affected by such an assumption as studied by Beven (1981) and Matonse and Kroll (2009). The steps illustrated in Fig. 2 were executed in the MATLAB environment as follows: i) the estimated aquifer hydraulic parameters (Ks, Kb, f, θs, θr) and the parameters of the Brooks-Corey and Richards equations were directly read from the look-up table. ii) The

CR IP T

meteorological and piezometer datasets were read at each simulation time step to estimate the actual evapotranspiration and capillary rise, respectively. iii) The initial values of ponding depth (if any) were read. iv) Based on the ponding or non-ponding situation, percolation through the transmission zone was estimated by transmission zone water balance using root

AN US

zone percolation, capillary rise, and bedrock leakage processes. v) The hillslope saturated storage and discharge were estimated. vi) The ponding and the transmission zone depths were

AC

CE

PT

ED

M

re-estimated to repeat the steps (iii)-(v).

19

ACCEPTED MANUSCRIPT

Start Ks, f, α, D, θr , θs , parameters of Brooks-Corey and Richards equations

Meteorological data, θrz , Dtz , h t , h b at time t ET0

AN US

1-D Richards Eq.

LAI Throughfall, Tf [Eq. (13)]

ETa = ETadj [Eq. (10)]

Infiltration capacity, Ic [Eq. (14)]

Percolation, RWB [Eq. (15)]

Percolation, RSP [Eqs. (19, 23)]

Infiltration, I [Eq. (12)]

Capillary rise, C [Eq. (7)]

θtz (z, t)

Transmission zone WB [Eq. (17) or Eq. (24)]

M

Defining channel network (Geomorphic analysis)

NWB or NSP

PT

ED

Hillslope delineation

t= t + 1

MS

LULC

CR IP T

No

Yes

DEM pre-processing by TauDEM (pit removal, computation of flow direction, contributing area and slope)

CE

ETpot [Eq. (9)]

PD>0?

DEM

HWF, w(x)

Yes

Intitialize PD, θtz , h(x) Re-estimated PD [Eq. (21)]

No

θ=θs ?

No

Nl [Eq. (26)]

hsB module: Eqs. (27, 28) Boundary Conditions: h = 0 at x= 0; Qgw = 0 at x=L h(x, t), Qgw (0, t) Last time step?

Yes

Stop

AC

Fig. 2. Workflow of the methodology followed for modeling the subsurface storage dynamics using different variants of the enhanced hsB model.

4. Study area and pedo-hydrogeologic database The study area considered herein is the ungauged Kanjhari reservoir catchment, which is a tributary of the Baitarani River in Odisha State in eastern India having a monsoon-type tropical

20

ACCEPTED MANUSCRIPT climate. This study area lies in between 21˚33'-21˚41'N latitudes and 85˚38'-85˚48'E longitudes, with an area of 125.04 km2. The elevation map of the study area, the locations of the experimental wells, their well-logs, and the double-ring infiltration tests are illustrated in Fig. 3. This area receives an average annual rainfall of about 1332 mm, with an average 75 numbers of rainy days per year, and a range of maximum daily rainfall intensity of 65.8–234.0

CR IP T

mm. More than 80% of the annual rainfall is generally received during June to September from the south-west monsoon depicting the wet period, which is the primary forcing for groundwater recharge. Similarly, the daily temperature varies in the range of 8–44 °C during the winter and

AN US

summer seasons.

W1 W2 W3 W4 W5 W6

2.5 5.0 7.5

CE

0.0

AC

Drilling depth (m)

b)

PT

ED

M

a)

10.0

Weathered Granite Clay with gravel Surface Clay

Fig. 3. (a) The Kanjhari reservoir catchment showing the elevation map, locations of experimental wells (W1–W6) at distances of 2640, 4860, 840, 3060, 1650 and 1170 m from the channel, respectively, and locations of double-ring infiltration tests; and (b) Well-logs.

21

ACCEPTED MANUSCRIPT The hard rock terrain of the Kanjhari reservoir catchment is made up of typical geomorphic landforms of isolated residual hillocks (inselbergs), structural and denuded hills, pediment and buried pediment, valleys, and linear ridges (Das et al., 1997). The subsurface soil of the catchment contains the secondary porosity that facilitates for the subsurface storage and groundwater movement in different weathered formations. This catchment is geo-hydraulically

CR IP T

interconnected, in which the near-surface weathered zone forms the unconfined aquifers, the source for flow of subsurface and groundwater in the deeper soil zones. The predominant type of surface soil distributed in the study area can be sub-grouped into coarse-loamy, fine–loamy,

PT

ED

M

AN US

and loamy (Fig. 4).

CE

Fig. 4. (a) LULC; and (b) soil maps of the Kanjhari reservoir catchment.

The daily-scale meteorological data of precipitation, relative humidity, and maximum and

AC

minimum temperatures for the study period at the nearest meteorology station, located at 21.37°N latitude and 85.35°E longitude, were collected from the Meteorology Centre, India Meteorological Department, Bhubaneswar. For the estimation of reference evapotranspiration, the solar radiation and wind speed data available with the National Aeronautics and Space Administration, USA (NASA, 2015) were used. Similarly, the daily soil moisture data of two

22

ACCEPTED MANUSCRIPT grid cells covering the study area at a spatial scale of 0.25° × 0.25° was obtained from the satellite-based reanalysis data products of SMOS-BEC (2012-2017). Figure 5 illustrates the timeseries of daily rainfall, soil moisture, and estimated actual evapotranspiration that were used for the daily-scale recharge estimation using the approach as shown in Fig. 2. The soil moisture contents in the root zone and the transmission zone were

CR IP T

measured by installing the Time Domain Refractometry (TDR) soil moisture probes up to three depths of 0.3 m, 1 m, 2 m and 3 m near all the well locations. For calibration and validation of the hsB model against the water table height above the impermeable layer, the daily-scale water table data were collected from two experimental open wells (W1, W2) along with the 15-

AN US

daily data from other four observation wells (W3–W6). Note that the catchment was considered ungauged regarding baseflow measurement; however, at least limited timeseries data of groundwater table is available for model setup. Further, the range of values of aquifer

M

hydraulic properties, viz., saturated hydraulic conductivity and specific yield were obtained from the pumping test results as published by Dhiman (2012). The soil samples were collected

ED

from the well locations, and grain-size analysis was performed to determine the aquifer

PT

formation material. From the dominant type of formation material, the aquifer hydraulic properties were estimated following the chart by Heath (1983). The LULC and soil type maps

CE

were obtained from the Odisha Space Application Centre (ORSAC), Bhubaneswar. To study the infiltration characteristics of the surface soil, double ring infiltrometer test was carried out

AC

at four locations (Fig. 3) in the catchment representing the different soil types (Fig. 4).

23

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 5. Rainfall, actual evapotranspiration (ETa) and soil moisture (SM) variables considered for estimation of recharge flux to the hillslope unconfined-aquifer.

5. Application of enhanced hsB model variants to a real-world catchment

M

5.1 Derivation of the hillslope width function (HWF)

ED

The hsB model operates on hillslopes as elements building the catchment. In this study, the entire catchment was considered as a single hillslope for applying the hsB model (e.g., Troch et

PT

al., 1994). The hillslope geomorphology, as regards of the HWF, was derived as the probability

CE

density of flow distances from the channel network (Bogaart and Troch, 2006). Hence, to obtain the HWF, the streamlines were first delineated from the available DEM as illustrated in

AC

Fig. 2 (e.g., Sahoo and Sahoo, 2018). The methodology of streamline delineation is illustrated in Fig. 2. For the computation of

flow direction, contributing area, and land slope, the DEM pre-processing was performed using the ‘Terrain analysis using Digital Elevation Model’ (TauDEM) version 5.1.2 (Tarboton, 2014). Channel network extraction involves the determination of the threshold contributing area, Ac; where, the channel networks were identified as those pixels having the contributing area larger than Ac. In this study, the value of Ac was determined by trial-and-error until the 24

ACCEPTED MANUSCRIPT extracted channel network, for a given estimate of Ac, resembled the natural drainage network. Then the flow distance was estimated by the D8 pour point algorithm (O’Callaghan and Mark, 1984); and, subsequently, from the probability density function (PDF) of flow distance, the HWF was estimated. Figure 6a illustrates the flow distance of different locations in the catchment from the channel network. The probability density of flow distances is illustrated in

CR IP T

Fig. 6b. In general, the hillslope width gets decreased from the streamlines towards the catchment divide (Fig. 6b) representing the divergent planform of the study area. 4

7

14000 a) 12000

Hillslope width, w (m)

5

b)

AN US

Northing (m)

6

4000

10000

x 10

3000

8000 6000

2000

4000

4 3 2

1000

1

5000 10000 Easting (m)

0

0 0

1000 2000 3000 4000 5000 6000 7000 Flow distance, x (m)

ED

0 0

M

2000

Fig. 6. (a) Map of flow distances from the channel network; and (b) hillslope width function

CE

PT

for the whole catchment estimated from the frequency of flow distances.

5.2. Real catchment application considering it both as gauged and ungauged regarding of

AC

water table elevation (WTE) 5.2.1. Evaluation of model parameters The bedrock slope angle, α of the catchment was determined from the well-logs at different places along the transect from the catchment boundary to the channel. The model parameters common to both the gauged and ungauged scenarios, as evaluated for the study area, are given in Table 1. The subsequent sections describe the parameterization of aquifer hydraulic property for the gauged and ungauged scenarios. 25

ACCEPTED MANUSCRIPT Table 1 hsB Model Parameters used for the Kanjhari Reservoir Catchment. Parameters Soil depth, D (m) Bedrock thickness, Db (m) Bedrock slope (%)

Estimation method Well-log Well-log Well-log

Values 10.00 2.50 10.56

Bedrock hydraulic conductivity, Kb (m/day)

Soil class of aquifer (Heath, 1983)

1.5 × 10-3

(m3/ m3) (m3/ m3)

RETC code (van Genuchten et al., 1991)

5.2.2. Gauged and ungauged catchment description

0.13 0.41 0.56 0.027 1.23

CR IP T

Residual moisture content, Saturated moisture content, Brooks-Corey index, B (-) Shape parameter, (cm-1) Shape parameter, n (-)

AN US

As defined by Sivapalan et al. (2003), an ungauged catchment is one where sufficient record of hydrological observation (both regarding data quality and quantity) is not available to aid in modeling the desired hydrological variable. Thus, a catchment may be gauged at one period

M

and ungauged at other periods. The desired variable can be precipitation, streamflow, and

groundwater modeling.

ED

sediment yield for surface water modeling; and water table data and baseflow rate in case of

PT

The developed enhanced versions of the hsB model were applied to the study area considering it both as gauged and ungauged catchments. In the present study, the catchment

CE

was considered as gauged for the period of modeling and ungauged otherwise, because information on WTE is available only during the periods of model calibration and validation.

AC

In a more generalized way, it was evaluated, how would have the developed model performed had the study area been an ungauged one. 5.2.3. Gauged catchment case The direct rainfall recharge (DRR), water balance (WB), water balance–surface ponding (WBSP), and water balance–surface ponding–bedrock leakage (WB-SP-BL) based enhanced hsB model variants were calibrated at daily-scale using the WTE of wells ‘W1’ and ‘W2’,

26

ACCEPTED MANUSCRIPT measured during the water year of 1st June 2012 to 31st May 2013. Out of these four variants, the best performing WB-SP-BL-based model was calibrated at 15-daily scale using the WTE of the wells W3, W4, and W5. Note that, as the period from 1st June to 31st May of the next year is considered as one water-year in India, the data periods for calibration and validation purposes were selected accordingly. It was considered that the hillslope unconfined aquifer is ̅ ). First, to deal with equifinality problem, the

CR IP T

made up of drainable pore space (

possible range of the soil class-specific aquifer parameters, Ks and f were obtained using the well-log information at each well location based on the works of Heath (1983) and Dhiman (2012). Subsequently, manual calibration was carried out by varying the Ks and f values within

AN US

the pre-estimated range. The spatial variation of the parameters was obtained using the inverse distance weighted (IDW) interpolation technique. To estimate the unknown information on soil parameters at any location (x,y) from the known information available at the surrounding

M

locations, the IDW interpolation method assumes that the surrounding locations, close enough to (x,y), are more alike than the distant surrounding locations as (Lu and Wong, 2008)

∑

⁄

(

)

( ⁄

ED

∑

)

( ⁄

)

∑

( ⁄

)

(30)

PT

∑

(29)

CE

where di = Euclidian distance of the centroid (x,y) of the grid of interest from the centroid of

AC

the ith grid with the known parameters; and ng = number of grid points with the known information on soil parameters. As the 1-D hsB model uses the grid-averaged parameters Ks(x) and f(x) in the x-direction only, the effective parameter values at a distance x from the stream network was estimated as the geometric mean of all the

and

grid values

located on the isoline at a flow distance x as [∏ [∏

]

(31a)

]

(31b) 27

ACCEPTED MANUSCRIPT where nx = number of grid points at a flow distance x. With the spatially varied values of Ks(x) and f(x), the model simulations were performed until the hsB solutions were closer to the observed WTE, resulting in the maximum NSE estimates. Subsequently, with the spatially varied calibrated values of Ks(x) and f(x), the model was validated against the daily-scale WTE observations at W1 and W2; and against 15-daily WTE

CR IP T

at W1–W6 locations. 5.2.4. Ungauged catchment case

In case of ungauged case of groundwater modeling, the available data is only well logs, from which the hydraulic parameters can be estimated. Due to lack of observed water table data, the

AN US

modelers are left with two options: i) Use the model parameter values from similar modeling studies on identical catchments by regionalization of parameters, or ii) Use the values obtained from aquifer formations through geophysical investigation. However, the second case might be

M

more logical in real-world scenarios since the possibility of finding a gauged catchment nearby the ungauged catchment is very remote. Hence, the aquifer hydraulic properties were used

ED

directly for modeling and assessing the model performance. As the aquifer hydraulic

PT

conductivity is a function of grain size and grain distribution characteristics, the well-logs were analyzed to determine the principal aquifer forming material (stratigraphic formation) by

CE

sieving. Subsequently, for the principal material of the stratigraphic layer, viz., coarse sand, medium sand, fine sand, silt and clay, the range of the Ks and f were obtained from the

AC

literature (Heath, 1983); and the representative material-specific values for that soil layer was estimated using the geometric mean of the respective lower and upper limiting values as obtained from the literature. The equivalent location-specific Ks and f values for the entire soil column were computed as ∑

(32a)

∑

28

ACCEPTED MANUSCRIPT ∑

(32b)

∑

where

equivalent saturated hydraulic conductivity of the aquifer material forming the

mth layer [LT-1];

equivalent porosity of the aquifer material forming the mth layer [LT-1];

thickness of the mth soil layer [L]; and nsl = total number of soil layers [-].

CR IP T

5.2.5. Model spin-up In most of the data-scarce catchments, where information on spatial distribution of groundwater elevation and soil moisture are unavailable, specifying the initial conditions is challenging (Ajami et al., 2014a). In such a case, model spin-up is generally performed by

AN US

running the model with a single year of precipitation forcing till equilibrium in modeling results is achieved. Herein, for both the calibration and validation, a 14-year spin-up period was used to achieve the equilibrium in the simulated WTE and also in the subsurface discharges with the same recharge input, so that the differences between the simulated WTE and also

M

between the subsurface discharges during the nth and the (n+1)th years were < 1% (e.g., Ajami

ED

et al., 2014b).

5.2.6. Model performance evaluation measures

PT

The criteria of Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970), coefficient of

CE

determination (R2), percentage bias (Pbias), and root mean square error (RMSE) were used for verifying the efficacy of the enhanced hsB model variants in simulating the corresponding

AC

observed hillslope saturated storage in terms of water table elevation for both the ungauged and gauged catchment setups. To consider the performance of a model satisfactory, the acceptable value of NSE, R2 and Pbias are >50%, >0.5 and ±25%, respectively (Moriasi et al., 2007).

6. Results and discussion 6.1. Performance of the hsB model considering the study area as gauged regarding water table elevation 29

ACCEPTED MANUSCRIPT Figure 7 shows the hsB model-predicted WTE by different approaches during calibration at well W1 (Fig. 7a) and well W2 (Fig. 7c), respectively; along with the corresponding Box-andWhisker plots of WTE (Figs. 7b, d). From the estimates of different performance evaluation measures of NSE, Pbias, R2, and RMSE presented therein, in general, it is envisaged that the accuracy of the water tables simulated by the enhanced hsB model variants varies in the order

CR IP T

DRR < WB < WB-SP < WB-SP-BL. The calibrated values of Ks and f for different enhanced hsB variants are mentioned in Table 2. Table 2

Calibrated and well log-based (values in parentheses) aquifer hydraulic parameters for the

W1

W2

f (-)

Ks

0.35 0.24 0.34 0.27 (0.25)

2.5 2.85 2.42 1.95 (1.70)

W3

f

AC

CE

PT

ED

DRR WB WB-SP WB-SP-BL

Ks (m/day) 3.14 2.17 2.57 2.25 (2.30)

M

hsB variants

AN US

gauged and ungauged catchment cases at different well locations, respectively.

30

0.41 0.37 0.35 0.38 (0.32)

Ks

f

W4 Ks

W5 f

Ks

f

Not calibrated Not calibrated Not calibrated 2.77 0.30 3.41 0.32 3.10 0.26

ACCEPTED MANUSCRIPT

0 50

7

h (m)

100

5 150

3

1 200 15.5.12 13.8.12 11.11.12 9.2.13 10.5.13 6 c) Well W2

5

4

AN US

h (m)

3

CR IP T

a) Well W1

9

Rainfall (mm/day)

11

2 1

M

0 15.5.12 13.8.12 11.11.12 9.2.13 10.5.13 Date (dd-mm-yy)

Fig. 7. (a, c) Calibration of the DRR, WB, WB-SP, and WB-SP-BL based enhanced hsB model

ED

variants at experimental wells W1 and W2 using the daily observed water table elevation

PT

timeseries, considering the catchment as gauged. Here h is the water table elevation above the bedrock; and (c, d) the respective box-and-whisker plots. The red crosses in the box-and-

CE

whisker plot represent the outliers.

AC

It can be envisaged from Fig. 7 that forcing the model with the recharge as a function of P,

ETa and Tf only, as in the case of the DRR approach, resulted in over-estimation of the WTE. This may be attributed to the over-estimation of recharge flux due to non-consideration of the dynamism of interception, evapotranspiration, infiltration, saturation and infiltration excess surface runoffs, and vadose-zone saturated and unsaturated flows. As inquired by the first research question, the DRR-based hsB approach might be useful in the preliminary stage of

31

ACCEPTED MANUSCRIPT model development or testing the model at the hillslope-scale. During calibration, the DRR approach overpredicted the observed WTE considerably throughout the year characterized by outliers in the Box-and-Whisker plot (Fig. 7b). Hence, the DRR-based hsB model is not practically applicable ‘as such’ to the real-ungauged catchment considered in this study, as confirmed by the NSE, Pbias, R2, and RMSE estimates of < -407%,

-22%,

0.60, and ≥

CR IP T

1.38 m, respectively, for both the wells W1 and W2 during the calibration.

Similarly, the WB approach-based enhanced hsB model underestimated the observed water table considerably during the calibration period, resulting in poor NSE, Pbias, R2, and RMSE estimates of < 38%, ≥ 5%, <0.87, and ≥0.37 m for both the wells, respectively. This could be

AN US

due to the non-consideration of ponding in paddy fields that, in turn, underestimated the monsoon recharge. Note that the infiltration rate under ponded phase is higher than that under non-ponded phase, which is not accounted for by the WB approach. This under-estimation of

AC

CE

PT

ED

M

recharge rate resulted in the underestimation of the water table.

32

ACCEPTED MANUSCRIPT

0 a) Well W1

100

7

300 400

5

500

3

600

1 700 15.5.13 13.8.13 11.11.13 9.2.14 10.5.14 c) Well W2

6

AN US

h (m)

4

CR IP T

200

h (m)

9

Rainfall (mm/day)

11

2

M

0 15.5.13 13.8.13 11.11.13 9.2.14 10.5.14 Date (dd-mm-yy)

Fig. 8. (a, c) Validation of different variants of the enhanced hsB model at the experimental

ED

wells W1 and W2 using the daily observed water table elevation timeseries considering the

PT

catchment as gauged; and (b, d) the corresponding Box-and-Whisker plots.

CE

Although similar underestimation is noticed during model validation at well-W1, the simulated water table overestimated the observed value at well-W2 after cessation of monsoon

AC

rainfall. During the validation period, precipitation was temporally distributed almost evenly with lesser intensity (hanging bars in Fig. 8a), which might have allowed the water to infiltrate leading to reduced number of days with surface saturation and ponding. Under such a situation, the WB-based model could have accounted for the recharge flux more accurately. However, during the calibration period, precipitation occurred in spells (hanging bars in Fig. 7a) with single high precipitation events, leading to surface saturation and ponding that cannot be 33

ACCEPTED MANUSCRIPT accounted for by the WB approach; hence, underestimated the recharge. The effect of recharge was reflected in the temporal behaviour of the water table. As well-W2 is surrounded by agricultural fields, this effect is more prominent there, resulting in an overestimation of the WTE. This further justifies the need for improvement in the WB approach-based hsB model variant by incorporating the saturation process modeling.

CR IP T

The limitation in the WB approach-based hsB model was addressed by incorporating the saturated phase processes in the WB-SP approach that performed better than the WB approach during the calibration period (Fig. 7); although, there was still overestimation by the WB-SP approach-based hsB model. Furthermore, during the validation period, the negative NSE

AN US

estimates of < -69% at both the W1 and W2 wells confirmed this overestimation (Figs. 8a, c). The overestimation of the simulated WTE by the WB-SP-based hsB model variant could be due to the inherent assumption of impermeable bedrock undertaken in the development of the

M

hsB model, which restricts any exchange flux with the underlying confined aquifer. However, in reality, the bedrock does not behave as an aquifuge and allows exchange flux between the

ED

unconfined and underlying confined aquifers, although at much a smaller rate, due to the lesser

PT

vertical hydraulic conductivity of the relatively compact bedrock. Note that a large amount of data illustrated as outliers in the whiskers of Figs. 8b and 8d also indicate the inefficacy of the

CE

WB-SP-based hsB model in reproducing the observed WTE during validation. The incorporation of bedrock leakage flux for recharge estimation in the WB-SP-BL

AC

approach-based hsB model variant enhanced the model performance significantly as evident from the estimates of NSE ≥ 70%, Pbias of ±2%, R2≥ 0.74, and RMSE

0.36 m during both

the calibration and validation phases (Figs. 7, 8). Hence, it is surmised that incorporating the bedrock leakage flux into the underlying confined aquifer along with the fluxes estimated by the water balance approach and surface-ponding (WB-SP-BL approach) from the paddy fields resulted in the best model efficiency, which answers the first research question fully.

34

ACCEPTED MANUSCRIPT The WB-SP-BL-based hsB model was also calibrated and validated at 15-daily scale at three well locations of W3, W4, and W5 (Figs. 9, 10) due to the heterogeneity of aquifer parameters; and further validated at other two well locations of W1 and W6 using the observed 15-daily water table data (Fig. 10a). As seen in Fig. 10a, this WB-SP-BL-based hsB approach performed satisfactorily in reproducing the observed WTE with the estimates of NSE >69% and R2

CR IP T

0.70 for all the wells. Hence, it is surmised that the revised hsB model framework

makes it suitable for studying the subsurface dynamics in scantily-gauged catchments at 15daily time-scale, as in case of many real-world scenarios, the WTE data is generally measured

AN US

at larger time-scales.

W4

6

NSE(%) = 70.40; Pbias(%) = 1.05; RMSE(m) = 0.21

W5

PT

8

ED

4 1.6.12 18.12.12 6.7.13 22.1.14 10.8.14 26.2.15 10

NSE(%) = 72.50; Pbias(%) = -0.13; RMSE(m) = 0.25

CE

6 1.6.12 18.12.12 6.7.13 22.1.14 10.8.14 26.2.15 Date (dd-mm-yy)

Simulated h (m)

Simulated h (m)

NSE(%) = 72.99; Pbias(%) = 0.20; RMSE(m) = 0.24

7 1.6.12 18.12.12 6.7.13 22.1.14 10.8.14 26.2.15 8

M

Wter Table Elevation (WTE), h (m)

9

10

R2 = 0.74

9 8

8

8

9

10

R2 = 0.74

7

6 6

Simulated h (m)

W3

7

8

R2 = 0.71

9 8 7 7

8

9

Observed h (m)

AC

Fig. 9. Calibration of the WB-SP-BL-based hsB model at 15-daily temporal scale using water table data of wells W3, W4, and W5 owing to the heterogeneity of aquifer parameters at these locations.

35

ACCEPTED MANUSCRIPT

7

R2 = 0.71

5

W4

R2 = 0.73

6 NSE(%) = 69.25; Pbias(%) = 1.32; RMSE(m) = 0.24

NSE(%) = 71.43; Pbias(%) = -0.12; RMSE(m) = 0.25

4 4 1.6.12 17.5.13 2.5.14 17.4.15 1.4.16 17.3.17 26.5.15 23.10.15 21.3.16 18.8.16 15.1.17 14.6.17 10 5 R2 = 0.75 R2 = 0.73 W2 W5 4 8 3 NSE(%) = 71.96; Pbias(%) = -0.48; RMSE(m) = 0.26 NSE(%) = 74.91; Pbias(%) = 0.58; RMSE(m) = 0.28

2 6 1.6.12 17.5.13 2.5.14 17.4.15 1.4.16 17.3.17 26.5.15 23.10.15 21.3.16 18.8.16 15.1.17 14.6.17 R2 = 0.71

10

W3

9

8

CR IP T

Wter Table Elevation (WTE), h (m)

6

8

W1

W6

R2 = 0.74

NSE(%) = 72.89; IOA(%) = 91.41; RMSE(m) = 0.25

NSE(%) = 72.02; Pbias(%) = -0.60; RMSE(m) = 0.24

AN US

6 7 26.5.15 23.10.15 21.3.16 18.8.16 15.1.17 14.6.17 1.6.12 17.5.13 2.5.14 17.4.15 1.4.16 17.3.17 Date (dd-mm-yy) Date (dd-mm-yy)

Fig. 10. (a) Reproduction of observed 15-daily water table elevation (WTE) timeseries by the WB-SP-BL-based hsB model variant.

M

6.2. Performance of the hsB model considering the study area as ungauged regarding water

ED

table elevation

Considering the catchment as ungauged, similar to the gauged case, the simulated water table

PT

in well-W1 was overestimated in the ungauged case by the DRR approach (Fig. 11). The WB, WB-SP, and WB-SP-BL approach-based hsB model variants also performed similarly as the

CE

gauged catchment case. In this ungauged scenario, the observed WTE were best reproduced by the WB-SP-BL approach with NSE > 69%, Pbias > -2.5%, R2 ≥ 0.75, and RMSE = 0.22 m

AC

(Figs. 11a, c) at both the wells W1 and W2. Hence, it is advocated that all the fluxes to the aquifer need to be accounted for while applying the hsB model to natural catchments. Furthermore, this reveals that the hsB model performed satisfactorily for both the ungauged and gauged catchment conditions. A comparative study between Figs. 8a (gauged case) and 11a (ungauged case) revealed that, in reproducing the WTE at well-W1, the NSE estimates by the DRR approach-based hsB 36

ACCEPTED MANUSCRIPT model for gauged and ungauged cases were -1967.47% and -2210.29% (unacceptable), respectively; whereas the corresponding estimates were 22.45% and 21.83% (unacceptable), 69.60% and -114.29% (unacceptable), and 70.84% and 67.70% (acceptable) for the WB-based hsB, WB-SP-based hsB, and WB-SP-BL-based hsB variants, respectively. Similar results were obtained in case of well-W2, where the performance of the gauged case was better than that of

CR IP T

the ungauged case. This could be due to the difference in model parameterization method, which also confirms that the present approach could be instrumental for estimating the soil hydraulic parameters in the ungauged catchments. 0 9

100 200 300

h (m)

7

400

5

500

3

M

600

Rainfall (mm/day)

a) Well W1

AN US

11

c) Well W2

6

h (m)

PT

4

CE

2

100 200 300 400 500 600 700

Rainfall (mm/day)

ED

1 700 15.5.13 13.8.13 11.11.13 9.2.14 10.5.14 0

AC

0 800 15.5.13 13.8.13 11.11.13 9.2.14 10.5.14 Date (dd-mm-yy)

Fig. 11. Reproduction of the daily observed groundwater table elevation timeseries by the DRR, WB, WB-SP, and WB-SP-BL approach-based hsB model variants for the experimental wells: (a, b) W1; and (c, d) W2, considering the study area as fully ungauged. Here h is the water table elevation above the bedrock.

37

ACCEPTED MANUSCRIPT 6.3. Uncertainty in the prediction of both the gauged and ungauged approaches The Quantile Regression approach (Weerts et al., 2011) was used for assessing the uncertainties of the WB-SP-BL-based hsB model using the calibrated and validated results for the gauged catchment scenario, and only the validated results for the ungauged scenario (Fig. 12). The results in Fig. 12 reveal that in most of the cases, the simulated values are well within

CR IP T

the range of 50% confidence interval (CI). The peak of the WTE hydrograph at well-W1 for the ungauged scenario was beyond 50% CI and within 90% CI. This indicates that the uncertainty associated with the simulated values of WTE is more in case of the ungauged

7

7

6

6

5

AN US

h (m)

catchment scenario than that in the gauged scenario.

7

6

5

11.8.12

7.11.12

3.2.13

4 2.5.13 15.5.13 5

4

3

7.11.13

3.2.14

2.5.14

3

2 15.5.12

11.8.12 7.11.12 3.2.13 Date (dd-mm-yy)

2.5.13

2 15.5.13

11.8.13

7.11.13

3.2.14

2.5.14

3

Validation W2 (gauged)

PT

Calibration W2

4 15.5.13 5

Validation W1 (ungauged)

4

ED

h (m)

4

11.8.13

M

4 15.5.12 5

5

Validation W1 (gauged)

Calibration W1

11.8.13 7.11.13 3.2.14 Date (dd-mm-yy)

Validation W2 (Ungauged) 2.5.14

2 15.5.13

11.8.13 7.11.13 3.2.14 Date (dd-mm-yy)

2.5.14

CE

Fig. 12. Quintile regression based uncertainty analysis of the WB-SP-BL-based hsB model

AC

variant for the gauged and ungauged scenarios (CI = confidence interval).

7. Conclusions The key findings from this study are as follows. 1. Soil-atmosphere interface and unsaturated zone processes are the integral parts of hillslope-based aquifer dynamics study for partitioning rainfall into the runoff, evapotranspiration, soil moisture, unsaturated-zone storage, and aquifer recharge flux.

38

ACCEPTED MANUSCRIPT Ignoring or over-simplifying any of these processes in the hsB model framework may invite substantial modeling uncertainty. 2. Bedrock leakage flux has to be accounted for in the hsB model framework while applying the model to a real-world catchment since no bedrock ideally behaves as an aquifuge. While performing the catchment-scale groundwater study both in the

CR IP T

unconfined and confined aquifers, quantification of this flux is essential; and it would provide insight into the shallow-deep aquifer interaction mechanism.

3. Simulations considering all the water flux estimation approaches revealed that the WBSP-BL-based enhanced hsB model could predict the temporal variation of the

AN US

subsurface water table in the experimental wells very well both at daily and 15-daily scales. This highlights the efficacy of the simplified WB-SP-BL-based hsB model for modeling the subsurface water in both the gauged and ungauged catchments, making it

PT

Acknowledgements

ED

as a future scope of research.

M

amenable for coupling with various land-surface schemes of the climate change models

We thank the Associate Editor, the Editor and three anonymous reviewers whose comments

CE

helped to improve the manuscript substantially in its present form. We thank Prof. Peter A. Troch, Dept. of Hydrology and Atmospheric Sciences, Univ. of Arizona, USA for providing

AC

the MATLAB codes of the original hsB model. The well-log, aquifer parameters, and water table data are obtained from the ICAR funded National Initiative on Climate Resilient Agriculture (NICRA) project (Sanction No. NICRA/CG/95/2011) and Central Groundwater Board, Bhubaneswar. The providers of other datasets as cited in the text are duly acknowledged. Thanks are also due to the fellowship provided to the first author by the Ministry of Human Resources Development, Government of India.

39

ACCEPTED MANUSCRIPT References Agrawal, M.K., Panda, S.N., Panigrahi, B., 2004. Modeling water balance parameters for rainfed rice. J. Irrig. Drain. Eng. 130(2), 129–139. doi:10.1061/(ASCE)07339437(2004)130:2(129). Ajami, H., Evans, J.P., McCabe, M.F., Stisen, S., 2014a. Reducing the spin-up time of integrated surface water-groundwater models. Hydrol. Earth Syst. Sci. 18(12), 5169–5179. doi:10.5194/hess-18-5169-2014.

CR IP T

Ajami, H., McCabe, M.F., Evans, J.P., Stisen, S., 2014b. Assessing the impact of model spinup on surface water-groundwater interactions using an integrated hydrologic model. Water Resour. Res. 50(3), 2636–2656. doi:10.1002/2013WR014258. Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Chapter 8 - ETc under soil water stress conditions, Crop evapotranspiration-Guidelines for computing crop water requirementsFAO Irrigation and drainage paper 56, Irrigation and Drainage. doi:10.1016/j.eja.2010.12.001.

AN US

Beven, K., 1981. Kinematic subsurface stormflow. Water Resour. Res. 17, 1419–1424. doi:10.1029/WR017i005p01419. Bogaart, P.W., Troch, P.A., 2006. Curvature distribution within hillslopes and catchments and its effect on the hydrological response. Hydrol. Earth Syst. Sci. 10(6), 925–936. doi:10.5194/hess-10-925-2006.

M

Boussinesq, J., 1877. Essai sur la théorie des eaux courantes. Mem. Acad. Sci. Inst. Fr., 23(1): 252-260.

ED

Broda, S., Larocque, M., Paniconi, C., 2014. Simulation of distributed base flow contributions to streamflow using a hillslope-based catchment model coupled to a regional-scale groundwater model. J. Hydrol. Eng. 19(5), 907–917. doi:10.1061/(ASCE)HE.19435584.0000877.

PT

Broda, S., Larocque, M., Paniconi, C., Haitjema, H., 2012. A low-dimensional hillslope-based catchment model for layered groundwater flow. Hydrol. Process. 26(18), 2814–2826. doi:10.1002/hyp.8319.

CE

Broda, S., Paniconi, C., Larocque, M., 2011. Numerical investigation of leakage in sloping aquifers. J. Hydrol. 409(1–2), 49–61. doi:10.1016/j.jhydrol.2011.07.035.

AC

Brodie, R.S., Hostetler, S., Slatter, E., 2008. Comparison of daily percentiles of streamflow and rainfall to investigate stream-aquifer connectivity. J. Hydrol. 349(1–2), 56–67. doi:10.1016/j.jhydrol.2007.10.056. Brooks, R., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrol. Pap. 3. Color. State Univ., Fort Collins. Carrillo, G., Troch, P.A., Sivapalan, M., Wagener, T., Harman, C., Sawicz, K., 2011. Catchment classification: Hydrological analysis of catchment behavior through processbased modeling along a climate gradient. Hydrol. Earth Syst. Sci. 15(11), 3411–3430. doi:10.5194/hess-15-3411-2011. Chen, S., Liu, C., 2002. Analysis of water movement in paddy rice fields (I) experimental studies. J. Hydrol. 260(1– 4), 206–215. doi:10.1016/S0022-1694(01)00615-1. Chen, S., Liu, C.-W., Huang, H.-C., 2002. Analysis of water movement in paddy rice fields (II) simulation studies. J. Hydrol. 268(1–4), 259–271. doi:10.1016/S0022-1694(01)00615-1. 40

ACCEPTED MANUSCRIPT Chen, Z., Govindaraju, R.S., Kavvas, M.L., 1994a. Spatial averaging of unsaturated flow equations under infiltration conditions over areally heterogeneous fields: 1. Development of models. Water Resour. Res. 30, 523–533. doi:10.1029/93WR02885. Chen, Z., Govindaraju, R.S., Kavvas, M.L., 1994b. Spatial averaging of unsaturated flow equations under infiltration conditions over areally heterogeneous fields 2. Numerical simulations. Water Resour. Res. 30, 535–548. doi:10.1029/93WR02884. Childs, E.C., 1971. Drainage of groundwater resting on a sloping bed. Water Resour. Res. 7(5), 1256– 1263. doi:10.1029/WR007i005p01256.

CR IP T

Das, S., Behera, S.C., Kar, A., Narendra, P., Guha, S., 1997. Hydrogeomorphological mapping in ground water exploration using remotely sensed data - A case study in Keonjhar district, Orissa. J. Indian Soc. Remote Sens. 25(4), 247–259. doi:10.1007/BF03019366. Dhiman, S. C., 2012. Aquifer systems of India. Central Ground Water Board, Ministry of Water Resources, River Development and Ganga Rejuvenation, Government of India. Dickinson, R.E., 1984. Modeling evapotranspiration for three-dimensional global climate models. Clim. Process. Clim. Sensit. 29, 58–72. doi:10.1029/GM029p0058.

AN US

Doorenbos, J., Pruitt, W.O., 1977. Guidelines for predicting crop water requirements, Irrig. and Drain. Pap. 24. Food and Agric. Organ., Rome. Duffy, C.J., 1996. A two-state integral-balance model for soil moisture and groundwater dynamics in complex terrain. Water Resour. Res. 32(8), 2421–2434. doi:10.1029/96WR01049.

M

Dupuit, J., 1863. Etudes théoriques et practiques sur le mouvement des eaux dans les canaux découverts et á travers les terrains perméables, 2nd ed., Dunod, Paris.

ED

Eagleson, P.S., 1978. Climate, soil, and vegetation: 3. A simplified model of soil moisture movement in the liquid phase. Water Resour. Res. 14(5), 722–730. doi:10.1029/WR014i005p00722.

PT

Ebel, B.A., Loague, K., Montgomery, D.R., Dietrich, W.E., 2008. Physics-based continuous simulation of long-term near-surface hydrologic response for the Coos Bay experimental catchment. Water Resour. Res. 44(7). doi:10.1029/2007WR006442. Fan, Y., Bras, R.L., 1998. Analytical solutions to hillslope subsurface storm flow and saturation overland flow. Water Resour. Res. 34(4), 921–927. doi:10.1029/97WR03516.

CE

Gardner, W.R., 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85(4), 228–232. doi:10.1097/00010694-195804000-00006.

AC

Graham, C.B., Woods, R.A., McDonnell, J.J., 2010. Hillslope threshold response to rainfall: (1) A field based forensic approach. J. Hydrol. 393(1–2), 65–76. doi:10.1016/j.jhydrol.2009.12.015. Grayson, R., Blöschl, G., 2000. Spatial Patterns in catchment hydrology: observations and modelling, Cambridge Univ. Press, New York. Haitjema, H.M., 1995. Analytic element modeling of groundwater flow. Academic Press, San Diego, California, USA. doi:10.1016/B978-012316550-3/50005-2. Hazenberg, P., Broxton, P., Gochis, D., Niu, G.-Y., Pangle, L.A., Pelletier, J.D., Troch, P.A., Zeng, X., 2016. Testing the hybrid-3-D hillslope hydrological model in a controlled environment. Water Resour. Res. 52, 1089–1107. doi:10.1002/2015WR018106. 41

ACCEPTED MANUSCRIPT Hazenberg, P., Fang, Y., Broxton, P., Gochis, D., Niu, G.-Y., Pelletier, J.D., Troch, P.A., Zeng, X., 2015. A hybrid-3D hillslope hydrologicalmodel for use in Earth system models. Water Resour. Res. 51(10), 8218–8239. doi:10.1002/2014WR016842. He, B., Wang, Y., Takase, K., Mouri, G., Razafindrabe, B.H.N., 2009. Estimating land use impacts on regional scale urban water balance and groundwater recharge. Water Resour. Manag. 23(9), 1863–1873. doi:10.1007/s11269-008-9357-2. Heath, R. C., 1983, Basic ground-water hydrology, U.S. Geol. Surv. Water Supply Pap., 2220.

CR IP T

Hilberts, A.G.J., Troch, P.A., Paniconi, C., 2005. Storage-dependent drainable porosity for complex hillslopes. Water Resour. Res. 41(6), 1–13. doi:10.1029/2004WR003725. Hilberts, A.G.J., Troch, P.A., Paniconi, C., Boll, J., 2007. Low-dimensional modeling of hillslope subsurface flow: Relationship between rainfall, recharge, and unsaturated storage dynamics. Water Resour. Res. 43(3). doi:10.1029/2006WR004964. Hilberts, A.G.J., Van Loon, E.E., Troch, P.A., Paniconi, C., 2004. The hillslope-storage Boussinesq model for non-constant bedrock slope. J. Hydrol. 291, 160–173. doi:10.1016/j.jhydrol.2003.12.043.

AN US

Hrachowitz, M., Savenije, H.H.G., Blöschl, G., McDonnell, J.J., Sivapalan, M., Pomeroy, J.W., Arheimer, B., Blume, T., Clark, M.P., Ehret, U., Fenicia, F., Freer, J.E., Gelfan, A., Gupta, H.V., Hughes, D.A., Hut, R.W., Montanari, A., Pande, S., Tetzlaff, D., Troch, P.A., Uhlenbrook, S., Wagener, T., Winsemius, H.C., Woods, R.A., Zehe, E., Cudennec, C., 2013. A decade of Predictions in Ungauged Basins (PUB)—A review. Hydrol. Sci. J. 58(6), 1198–1255. doi:10.1080/02626667.2013.803183.

M

Ines, A.V.M., Mohanty, B.P., 2008. Near-surface soil moisture assimilation for quantifying effective soil hydraulic properties using genetic algorithm: 1. Conceptual modeling. Water Resour. Res. 44(6). doi:10.1029/2007WR005990.

ED

Jha, R.K., Sahoo, B., Panda, R.K., 2017. Modeling the water and Nitrogen transports in a soilpaddy-atmosphere system using HYDRUS-1D and lysimeter experiment. Paddy Water Environ. 15(4). 831–846. https://doi.org/10.1007/s10333-017-0596-9.

PT

Kale, R.V., Sahoo, B., 2011. Green-Ampt infiltration models for varied field conditions: A revisit. Water Resour. Manage. 25(14). 3505-3536. doi: 10.1007/s11269-011-9868-0.

CE

Khepar, S.D., Yadav, A.K., Sondhi, S.K., Siag, M., 2000. Water balance model for paddy fields under intermittent irrigation practices. Irrig. Sci. 19, 199–208. doi:10.1007/PL00006713.

AC

Kong, J., Shen, C., Luo, Z., Hua, G., Zhao, H., 2016. Improvement of the hillslope-storage Boussinesq model by considering lateral flow in the unsaturated zone. Water Resour. Res. 52(4), 2965–2984. doi:10.1002/2015WR018054. Kurylyk, B.L., MacQuarrie, K.T.B., Voss, C.I., 2014. Climate change impacts on the temperature and magnitude of groundwater discharge from shallow, unconfined aquifers. Water Resour. Res. 50(4), 3253–3274. doi:10.1002/2013WR014588. Larocque, M., Fortin, V., Pharand, M.C., Rivard, C., 2010. Groundwater contribution to river flows – using hydrograph separation, hydrological and hydrogeological models in a southern Quebec aquifer. Hydrol. Earth Syst. Sci. Discuss. 7(5), 7809–7838. doi:10.5194/hessd-7-7809-2010. List, F., Radu, F.A., 2016. A study on iterative methods for solving Richards’ equation. Comput. Geosci. 20, 341–353. doi:10.1007/s10596-016-9566-3. 42

ACCEPTED MANUSCRIPT Liu, J., Han, X., Chen, X., Lin, H., Wang, A., 2016. How well can the subsurface storagedischarge relation be interpreted and predicted using the geometric factors in headwater areas? Hydrol. Process. doi:10.1002/hyp.10958. Lu, G.Y., Wong, D.W., 2008. An adaptive inverse-distance weighting spatial interpolation technique. Comput. Geosci. 34, 1044–1055. doi:10.1016/j.cageo.2007.07.010. Matonse, A.H., Kroll, C., 2009. Simulating low streamflows with hillslope storage models. Water Resour. Res. 45(1). doi:10.1029/2007WR006529.

CR IP T

Matonse, A.H., Kroll, C.N., 2013. Applying hillslope-storage models to improve low flow estimates with limited streamflow data at a watershed scale. J. Hydrol. 494, 20–31. doi:10.1016/j.jhydrol.2013.04.032. Meerveld, I.T. van, Weiler, M., 2008. Hillslope dynamics modeled with increasing complexity. J. Hydrol. 361(1–2), 24–40. doi:10.1016/j.jhydrol.2008.07.019. Merill, L., Tonjes, D.J., 2014. A review of the hyporheic zone, stream restoration, and means to enhance denitrification. Crit. Rev. Environ. Sci. Technol. 44(21), 2337–2379. doi:10.1080/10643389.2013.829769.

AN US

Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Bingner, R.L., Harmel, R.D., and Veith, T.L., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. T ASABE. 50 (3), 885-900. doi: 10.13031/2013.23153. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522. doi:10.1029/WR012i003p00513.

M

NASA, 2015. Climatology Resource for Agroclimatology, , Accessed on 16 April 2015. Nash, J.E., Sutcliffe, J. V., 1970. River flow forecasting through conceptual models part I - A discussion of principles. J. Hydrol. 10(3), 282–290. doi:10.1016/0022-1694(70)90255-6.

ED

O’Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Comput. Vision, Graph. Image Process. 28(3), 323–344. doi:10.1016/S0734-189X(84)80011-0.

PT

Ouyang, Y., 2012. A potential approach for low flow selection in water resource supply and management. J. Hydrol. 454, 56–63. doi:10.1016/j.jhydrol.2012.05.062.

CE

Pan, Y., Weill, S., Ackerer, P., Delay, F., 2015. A coupled stream flow and depth-integrated subsurface flow model for catchment hydrology. J. Hydrol. 530, 66–78. doi:10.1016/J.JHYDROL.2015.09.044.

AC

Paniconi, C., Troch, P.A., Emiel van Loon, E., Hilberts, A.G.J., 2003. Hillslope-storage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 2. Intercomparison with a three-dimensional Richards equation model. Water Resour. Res. 39(11), 1317. doi:10.1029/2002wr001730. Parlange, J.-Y. and Haverkamp, R., 1989. Infiltration and ponding time. In Unsaturated Flow in Hydrologic Modeling, Theory and Practice. ed. H.J. Morel-Seytoux. 95-126. Kluwer Academic Publishers. Boston, MA. Philip, J.R., 1957. The theory of infiltration: 4. Sorptivity and algebraic infiltration equations. Soil Science. 84, 257-264. Sahoo, S., Sahoo, B., 2018. Modelling the variability of hillslope drainage using grid-based hillslope width function estimation algorithm. ISH J. Hydraul. Eng. 1–8. 43

ACCEPTED MANUSCRIPT doi:10.1080/09715010.2018.1441750. Sahoo, B., Sahoo, S., Swain, R., 2017. Modeling the river-aquifer flow-interaction using a coupled hsB-VPMM approach, In: ASCE-EWRI World Environmental and Water Resources Congress 2017: Groundwater, Sustainability, and Hydro-Climate/Climate Change, Sacramento, CA, USA, May 21–25, 2017, 172-182. doi:10.1061/9780784480618.0017.

CR IP T

Sivapalan, M., Takeuchi, K., Franks, S.W., Gupta, V.K., Karambiri, H., Lakshmi, V., Liang, X., McDonnell, J.J., Mendiondo, E.M., O’Connell, P.E., Oki, T., Pomeroy, J.W., Schertzer, D., Uhlenbrook, S., Zehe, E., 2003. IAHS Decade on Predictions in Ungauged Basins (PUB), 2003–2012: Shaping an exciting future for the hydrological sciences. Hydrol. Sci. J. 48(6), 857–880. doi:10.1623/hysj.48.6.857.51421. SMOS-Barcelona Expert Center (SMOS-BEC), 2012-2017, Land variables, Level 2 Soil moisture user data product, , updated regularly, Accessed on 08 March 2016.

AN US

Swain, R., Sahoo, B., 2015. Variable parameter McCarthy–Muskingum flow transport model for compound channels accounting for distributed non-uniform lateral flow. J. Hydrol. 530, 698–715. doi:10.1016/j.jhydrol.2015.10.030. Tarboton, D., 2014, TauDEM documentation, , Accessed on 08 January 2015. Troch, P.A., Carrillo, G., Sivapalan, M., Wagener, T., Sawicz, K., 2013. Climate-vegetationsoil interactions and long-term hydrologic partitioning: signatures of catchment coevolution. Hydrol. Earth Syst. Sci. 17(6), 2209–2217. doi:10.5194/hess-17-2209-2013.

ED

M

Troch, P.A., Paniconi, C., Emiel van Loon, E., 2003. Hillslope-storage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 1. Formulation and characteristic response. Water Resour. Res. 39(11), 1316. doi:10.1029/2002wr001728. Troch, P.A., Mancini, M., Paniconi, C., Wood, E.F., 1993. Evaluation of a distributed catchment scale water balance model. Water Resour. Res. 29, 1805–1817. doi:10.1029/93WR00398.

CE

PT

Troch, P.A., Smith, J.A., Wood, E.F., de Troch, F.P., 1994. Hydrologic controls of large floods in a small basin: central Appalachian case study. J. Hydrol. 156(1–4), 285–309. doi:10.1016/0022-1694(94)90082-5.

AC

Troch, P., Van Loon, E., Hilberts, A., 2002. Analytical solutions to a hillslope-storage kinematic wave equation for subsurface flow. Adv. Water Resour. 25(6), 637–649. doi:10.1016/S0309-1708(02)00017-9. Tromp-van Meerveld, H.J., Peters, N.E., McDonnell, J.J., 2007. Effect of bedrock permeability on subsurface stormflow and the water balance of a trenched hillslope at the Panola Mountain Research Watershed, Georgia, USA. Hydrol. Process. 21(6), 750–769. doi:10.1002/hyp.6265. van Genuchten, M.T., 1980. A Closed-form Equation for Prediccting Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J. 44, 892–898. doi: 10.2136/sssaj1980.03615995004400050002x. van Genuchten, M.T., Leij, F.J., Yates, S.R., 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. U. S. Environ. Prot. Agency, Ada, Oklahoma. doi:10.1002/9781118616871. 44

ACCEPTED MANUSCRIPT Weerts, A.H., Winsemius, H.C., Verkade, J.S., 2011. Estimation of predictive hydrological uncertainty using quantile regression: examples from the National Flood Forecasting System (England and Wales). Hydrol. Earth Syst. Sci. 15, 255–265. doi:10.5194/hess-15255-2011.

AC

CE

PT

ED

M

AN US

CR IP T

Woods, R., Sivapalan, M., 1999. A synthesis of space-time variability in storm response: Rainfall, runoff generation, and routing. Water Resour. Res. 35(8), 2469–2485. doi:10.1029/1999WR900014.

45