The Wageningen Lowland Runoff Simulator ( WALRUS ) : a lumped rainfall – runoff model for catchments with shallow groundwater

We present the Wageningen Lowland Runoff Simulator (WALRUS), a novel rainfall–runoff model to fill the gap between complex, spatially distributed models which are often used in lowland catchments and simple, parametric (conceptual) models which have mostly been developed for sloping catchments. WALRUS explicitly accounts for processes that are important in lowland areas, notably (1) groundwater–unsaturated zone coupling, (2) wetnessdependent flow routes, (3) groundwater–surface water feedbacks and (4) seepage and surface water supply. WALRUS consists of a coupled groundwater–vadose zone reservoir, a quickflow reservoir and a surface water reservoir. WALRUS is suitable for operational use because it is computationally efficient and numerically stable (achieved with a flexible time step approach). In the open source model code default relations have been implemented, leaving only four parameters which require calibration. For research purposes, these defaults can easily be changed. Numerical experiments show that the implemented feedbacks have the desired effect on the system variables.


Introduction
Lowlands, especially those in river deltas, are often densely populated and centers of agricultural production, economic activity and transportation. Therefore, socioeconomic consequences of natural hazards are specially large in these areas. In addition, the lack of topography increases their vulnerability to flooding (coastal, fluvial 20 and pluvial), climate change, and deterioration of water quality. To mitigate natural and human disasters, hydrological models can be used by water managers as a tool for risk assessment and infrastructure design. Lowlands, defined here as areas with shallow groundwater tables, exist all over the world: 13 % of world's land surface has groundwater tables shallower than 2 m and 22 % shallower than 4 m ( Fig. 1 Water managers in lowland areas often use complex hydrological models. MIKE-SHE (Refsgaard and Storm, 1995), HEC-RAS (Brunner, 1995) and SOBEK (Deltares, 2013) have detailed schematisations of surface water networks to simulate the complex flow routing in intensively drained areas. HYDRUS (Simůnek et al., 2008) and SWAP (Van Dam et al., 2008) have detailed vertical schematisations to simulate unsaturated-5 saturated zone coupling. Regional groundwater models, such as MODFLOW (McDonald and Harbaugh, 1984), account for seepage and lateral groundwater flow. Combinations of several of these models can be used to account for groundwater-surface water feedbacks, such as SHE (Abbott et al., 1986), HydroGeoSphere (Therrien et al., 2006) SIMGRO (Querner, 1988;Van Walsum and Veldhuizen, 2011) or NHI (Prinsen and 10 Becker, 2011). However, complex models have important disadvantages and simple models important advantages: 1. Overparameterisation -model parameters account for differences in response times or recession shapes between catchments with the same dominant processes (represented by the model structure). With too many parameters, an in- 15 appropriate model structure can be compensated for by mathematically fitting the model to the calibration data (Kirchner, 2006). An overparameterised model may perform well during calibration, but unsatisfactorily during validation (Perrin et al., 2001) and in different (future) climate regimes (e.g. Seibert, 1999).
2. Parameter identification -the risk of parameter dependence and equifinality 20 (where different combinations of parameter values to lead to similar results, Beven and Binley, 1992;Uhlenbrook et al., 1999) increases with the number of parameters. With one objective function, only typically 3 to 5 parameters can be identified (Jakeman and Hornberger, 1993;Beven, 1989). Multi-objective calibration allows more parameters to be calibrated (e.g. Gupta et al., 1998 which can be used both for multi-year water balance studies and for single rainfallrunoff events. The model was designed to have an understandable model structure that incorporates the most important processes and feedbacks, with fewer than 6 parameters of which the values do not change with the temporal resolution at which the model is run.

5
In this paper we present WALRUS. Section 2 contains descriptions of two contrasting lowland field sites which were used in the model development. In Sect. 3, we describe several lowland-specific hydrological processes and their representation in WALRUS. In Sect. 4 we explain the model structure in detail. Section 5 contains the implementation of the model in code and Sect. 6 the conclusions. A detailed model evaluation 10 is discussed in an accompanying manuscript available at Hydrology and Earth System Sciences Discussions (Brauer et al., 2014).

Two contrasting lowland catchments
Field experience and data from two contrasting field sites in the Netherlands have been used to develop the model structure. The Hupsel Brook catchment in the east 15 has slightly sloping surfaces and drainage is driven by gravity (freely draining). The Cabauw polder in the west is flat and its water levels are controlled with weirs and supply of surface water. The Hupsel Brook catchment is 6.5 km 2 , its soils consist of 0.2-11 m of loamy sand on an impermeable clay layer and land cover is mostly grass (59 %) and some maize (33 %). The Cabauw polder is 0.5 km 2 , its soils consist of about 20 70 cm heavy clay on peat and land cover is 80 % grass, 15 % maize and 5 % surface water. A more detailed description of both catchments and observations is presented in an accompanying paper (Brauer et al., 2014). From the Hupsel Brook catchment, we used combined observations of groundwater and soil moisture (from a neutron probe at 12 depths, ranging from 0.15 to 2.05 m) at six 25 locations, which represent the spatial variability in the catchment well. Potential evapotranspiration was estimated with the method of Thom and Oliver (1977). During the GMDD growing seasons  of 1976 through 1982 daily sums of actual evapotranspiration (ET act ) have been computed with the energy budget method: net radiation was measured and wind and temperature profile measurements were used to estimate sensible and ground heat flux. Evapotranspiration was then estimated as residual of the energy budget (for more information see Stricker and Brutsaert, 1978). 5 In addition, we used data of 1993: discharge data measured with a type of H-flume, groundwater depths measured at the meteorological station and total phosphorus, nitrate and chloride concentrations measured at the catchment outlet.
From the Cabauw polder, we used daily soil moisture data from four arrays of six TDR sensors between 5 and 73 cm deep from the period 2003-2010. Groundwater 10 depth was measured at the same location. Potential evapotranspiration ET pot was estimated with the method of Makkink (1957) and ET act is determined by measuring net radiation, ground heat flux and Bowen ratio (with an eddy covariance set-up) and closing the energy balance (Beljaars and Bosveld, 1997). Because ET act estimated with this method was on average 4 % higher than ET pot during well-watered conditions, we 15 divided ET act by 1.04.

Characteristics of lowland catchments
In this section we discuss some characteristics which affect hydrological processes in lowland catchments. We discuss how they are represented in some widely used rainfall-runoff models and how they are accounted for explicitly in WALRUS.

Groundwater-unsaturated zone coupling
Whereas in most models percolation is assumed to be driven by downward gravitational forces only, the vertical profile of moisture content in lowland soils is influenced by capillary forces associated with the presence of a shallow groundwater table. Percolation is slower and evapotranspiration remains high in dry periods, because storage Introduction  Hopmans and van Immerzeel, 1988;Stenitzer et al., 2007). Therefore, the vadose zone and the groundwater zone form a tightly coupled system and feedbacks should be included in models for lowland catchments (Chen and Hu, 2004). In addition, when groundwater rises to the soil surface, the unsaturated zone shrinks and its storage capacity decreases. It is therefore important 5 to include a dynamic unsaturated zone in the model, which is influenced by the surface fluxes precipitation and evapotranspiration as well as by the (dynamic) groundwater table below. Many conceptual rainfall-runoff models, e.g. HBV, the Sacramento model and the Wageningen Model, contain separate reservoirs for soil moisture and groundwater, al-10 lowing only downward movement of groundwater without considering feedbacks. One version of PDM does reduce recharge when the soil ceases to be freely draining (Moore, 2007). Catchments can also be simplified to a single nonlinear reservoir, without discriminating between the saturated and unsaturated zone (Kirchner, 2009), which yielded limited success in the lowland Hupsel Brook catchment (Brauer et al., 2013). 15 Quasi-steady state approaches have also been developed for implementation in distributed models e.g. by Koster et al. (2000); Bogaart et al. (2008) and Van Walsum and Groenendijk (2008).
WALRUS contains one soil reservoir, which can be divided effectively by the (dynamic) groundwater table into a groundwater zone and a vadose zone. The condition 20 of this soil reservoir is described by two strongly dependent variables: the groundwater depth and the storage deficit (the effective thickness of empty pores). The water balance in the whole soil reservoir is maintained through the storage deficit, while the groundwater depth is only used as pressure head to compute the groundwater drainage flux. The groundwater table reacts to changes in storage deficit (after rain or evapotran-25 spiration) by moving towards an equilibrium between storage deficit and groundwater depth. Although the soil moisture profile is not simulated explicitly, this implementation enables upward movement of groundwater when the top soil has dried through evapotranspiration. Introduction

Shallow groundwater and plant water stress
Vegetation in lowland catchments is hardly affected by water stress, which is one of the drivers for high agricultural production. Water is not only made available through physical processes (capillary rise), but also through physiological ones: when plants have exhausted the readily available moisture in the top soil, deeper roots are used 5 (Zencich et al., 2002), and vertical roots grow deeper (Canadell et al., 1996;Weir and Barraclough, 1986;Teuling et al., 2006) and more quickly (Zeng et al., 2013). Because plants adapt to spatial variability in moisture content, water uptake and its vertical distribution depend primarily on the availability of moisture in the whole root zone (Jarvis, 1989). As roots in lowlands often extend to close to the groundwater table, plants can 10 adapt fully to dry periods and evapotranspiration reduction hardly occurs (Schenk and Jackson, 2002). This dynamic system of different plant species with varying stages of root development and spatially and temporally varying groundwater depths is complex, but not all complexity may be necessary to include in a model for runoff simulations (Van der Ploeg et al., 2012). 15 In some rainfall-runoff models for areas with deep water tables, a separate root zone is included, e.g. in SWAT and TOPMODEL, which exhibits a different behaviour than the unsaturated zone below. We assume that in lowlands, this distinction cannot be made because the whole unsaturated zone can be used by plant roots. The variation of plant species within a catchment is sometimes represented by running a model for 20 different vegetation types separately and multiplying the resulting discharge output with the fraction of that vegetation type (Van Dam et al., 2008).
We assume that in lowlands the whole unsaturated zone can be used by plant roots. In WALRUS, spatial variation in vegetation cover is not modelled explicitly to reduce the risk of overparameterisation (data on the detailed functioning of the system are scarce) 25 and because the entire system of feedbacks between plants and water is complex on small scales, but likely less complex on larger scales. The effect of vegetation diversity on potential evapotranspiration can be accounted for by preprocessing. Introduction

Wetness-dependent flowroutes
When the soil wetness increases, different flowpaths are activated: from groundwater flow (Hall, 1968), to natural macropores (Mosley, 1979;Beven andGermann, 1982, 2013;McDonnell, 2003) and drainpipes (Tiemeyer et al., 2007;Rozemeijer et al., 2010a;Van der Velde et al., 2010b) and eventually to surface runoff (Dunne and Black, 5 1970; Brauer et al., 2011;Appels et al., 2011). Figure 2 provides examples of discharge mechanisms in lowland catchments at different scales. Stream water chemistry is increasingly being used to detect hydrological flow paths (e.g. Soulsby et al., 2004;Tetzlaff et al., 2007;Delsman et al., 2013). Records of phosphorus, nitrate and chloride concentrations measured at the outlet of the Dutch Hupsel 10 Brook catchment confirm the activation of different flow routes at different stages of catchment wetness (Fig. 3). The activation of drainpipes in September is indicated by increasing nitrate concentrations and overland flow during peaks by decreasing chloride and nitrate concentrations and increasing phosphorus concentrations.
The contribution of preferential flow and macropore flow can be considerable and 15 needs to be accounted for in the model structure (Beven and Germann, 1982;Weiler and McDonnell, 2004;Hansen et al., 2013). Drainpipes can be viewed as man-made macropores (Herrmann and Duncker, 2008) and account for a large fraction (up to 80 %) of drainage in lowlands ( Van der Velde et al., 2011;Turunen et al., 2013). When local storage thresholds are exceeded and quick flowpaths are activated, a sudden 20 increase in local discharge occurs (the fill and spill hypothesis by Tromp-van Meerveld and McDonnell, 2006), but at the catchment scale, sudden changes in discharge are hardly ever observed, because spatial variability in groundwater depth, drainpipe depth and microtopography cause these thresholds to be reached at different moments at different locations (Appels, 2013).

25
Parametric models often divide water between fast and slow routes. In the GR4J model (Perrin et al., 2003) this division is fixed, the PDM model (Moore, 1985) uses a wetness-dependent probability distribution to express the spatial variability in quickflow contribution, and in the Wageningen Model (Stricker and Warmerdam, 1982) the division depends on groundwater storage. In WALRUS, the storage deficit determines the division of rain between a soil reservoir (slow routes: infiltration, percolation and groundwater flow) and a quickflow reservoir (quick routes: drainpipe, macropore and overland flow).

Groundwater-surface water feedbacks
Surface water is an important feature in lowland landscapes (Fig. 2). The aim of manmade drainage networks in controlled catchments is to optimize groundwater depths by adjusting surface water levels (Krause et al., 2007). During discharge peaks, backwater feedbacks can occur and high surface water levels reduce groundwater drainage or 10 may even cause infiltration (Brauer et al., 2011).
Most parametric rainfall-runoff models do not simulate surface water levels, and therefore parametric models for vertical flow in the unsaturated zone are often coupled to a distributed groundwater model for studies on groundwater-surface water interactions (Krause and Bronstert, 2007;Sophocleous and Perkins, 2000;Lasserre et al., 15 1999; Van der Velde et al., 2009).
In WALRUS, surface water forms an integral part of the model structure. Drainage depends on the difference in water level between the surface water and groundwater reservoirs (rather than groundwater levels alone), allowing for feedbacks and infiltration of surface water into the soil.

Seepage and surface water supply
Regional groundwater flow is common in lowland areas and upward or downward seepage can be a large term in the water budget. Surface water is often supplied to raise groundwater levels for optimal crop growth, to avoid algal blooms (by maintaining flow velocity), to reduce brackish seepage in coastal areas below sea level, or to prevent Introduction peat oxidation. In addition, the water can be removed from the catchment by pumping (Van den Eertwegh et al., 2006;Te Brake et al., 2013;Delsman et al., 2013). Usually, distributed models are used for regional groundwater flow (MODFLOW), surface water supply and extraction (MIKE-SHE, SOBEK) and control operations (Van Andel et al., 2010) and the effect of changing surface water levels on runoff generation 5 is not taken into account.
In WALRUS, seepage and surface water supply or extraction are added to or subtracted from the soil or surface water reservoir. These external fluxes affect the whole system through the groundwater-surface water feedbacks and saturated-unsaturated zone coupling described in the previous sections.

Model description
In this section we provide a detailed description of all model components: reservoirs, states, fluxes and feedback mechanisms. The model contains several relations between model variables which can be specified by the user. We implemented defaults for these relations, such that WALRUS can be used directly by practitioners, while re- 15 taining the option to change them for research purposes.

General overview
WALRUS is a water balance model with three reservoirs and fluxes between the reservoirs. The model can be split into 5 compartments (Fig. 4, for abbreviations of variables, see Tab. 1): 20 1. Land surface -at the land surface, water is added to the different reservoirs by precipitation (P ). A fixed fraction is led to the surface water reservoir (P S ). The soil wetness index (W ) determines which fraction of the remaining precipitation percolates slowly through the soil matrix (P V ) and which fraction flows towards the surface water via quick flow routes (P Q ). Water is removed by evapotranspiration from the vadose zone (ET V ) and surface water reservoir (ET S ).
2. Vadose zone within the soil reservoir -the vadose zone is the upper part of the soil reservoir and extends from the soil surface to the dynamic groundwater table (d G ), including the capillary fringe. The dryness of the vadose zone is charac-5 terised by a single state: the storage deficit (d V ), which represents the effective volume of empty pores per unit area. It controls the evapotranspiration reduction (β) and the wetness index (W ).
3. Groundwater zone within the soil reservoir -the phreatic groundwater extends from the groundwater depth (d G ) downwards, thereby assuming that there is no 10 shallow impermeable soil layer and allowing groundwater to drop below the depth of the drainage channels (c D ) in dry periods. The groundwater table responds to changes in the unsaturated zone storage and determines together with the surface water level groundwater drainage or infiltration of surface water (f GS ).

4.
Quickflow reservoir -all water that does not flow through the soil matrix, passes 15 through the quickflow reservoir to the surface water (f QS ). This represents macropore flow through drainpipes, animal burrows and soil cracks, but also local ponding and overland flow.

5.
Surface water reservoir -the surface water reservoir has a lower boundary (the channel bottom c D ), but no upper boundary. Discharge (Q) is computed from the 20 surface water level (h S ).
6. External fluxes -water can be added to or removed from the soil reservoir by seepage (f XG ) and to/from the surface water reservoir by surface water supply or extraction (f XS ).
The area of the surface water reservoir a S is the fraction of the catchment covered 25 by ditches and channels, which is supplied by the user and can generally be derived 1369 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | from maps. The area of the soil reservoir a G is the remainder (1 − a S ). The area of the quickflow reservoir is taken equal to a G , but this is arbitrary since the outflow depends on the volume of water in the reservoir and a parameter (see Sect. 4.8). In the following sections the processes occurring within and between each compartment are discussed.

5
Because the soil reservoir has no lower boundary and the surface water reservoir no upper boundary, the groundwater depth d G is measured with respect to the soil surface and the surface water level h S with respect to the channel bottom. The channel bottom c D , with respect to the soil surface, is used to compute the difference in level, which is necessary for the computation of groundwater drainage. The quickflow reservoir level 10 h Q is measured with respect to the bottom of that reservoir. The storage deficit d V is an effective thickness, instead of a level or depth.

Precipitation and wetness index
Precipitation (P ) is divided between the 3 reservoirs: a fixed fraction a S falls directly onto the surface water (P S ) and the remainder is divided between the vadose zone (P V ) 15 and the quickflow reservoir (P Q ). The wetness index (W ) gives the fraction of the rainfall that is led to the quickflow reservoir and ranges from 0 (dry -all water is led to the soil reservoir) to 1 (wet -all water is led to the quickflow reservoir). The wetness index is a function of storage deficit (d V , Sect. 4.4). This relation can be supplied by the user, but as default a cosine function has been implemented, which starts at 1 when the soil 20 is completely saturated (d V = 0) and drops to zero when d V is equal to the wetness parameter c W [mm], which has to be calibrated: The effect of this variable division between quick and slow flow paths is investigated by running WALRUS twice for an artificial example: with and without the variable W . Six rainfall events with a duration of one day and an intensity of 2 mm h −1 , separated by four dry days yield the same quickflow f QS and discharge Q response when the divider is not depending on soil moisture storage, but f QS and Q increase in case of 5 a wetness-dependent divider (Fig. 5). The storage deficit d V decreases quickly during rainfall events and increases slowly in dry intervals. The variable wetness index W follows d V without delay and the groundwater depth d G responds with a delay caused by the unsaturated zone (represented by its relaxation time parameter c V , see Sect. 4.6).
With a variable W , the groundwater level rises quickly at first, but more slowly at the 10 end, because less water is led to the soil reservoir when it is already wet. This numerical experiment shows that the variable wetness index ensures that WALRUS can simulate feedbacks between groundwater, vadose zone and quickflow and that variables at the soil surface do not only influence variables in the ground (as in most models), but also the other way around.

Evapotranspiration
Evapotranspiration (ET) takes place from the surface water reservoir (ET S ) and the vadose zone (ET V ). The actual evapotranspiration from the vadose zone depends on the potential evapotranspiration rate and the storage deficit (Fig. 6). The relation between the evapotranspiration reduction factor β and the storage deficit can be supplied by the 20 user. As a default, a two-parameter function has been implemented: The evapotranspiration reduction factor approaches one (no reduction) when the soil is saturated and decreases with storage deficit: first slowly, then more quickly and then more slowly again (although this end of the curve is never reached in practice). Eq. (2 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | d V the reduction factor is 0.5 (the inflection point). Note that Eq. (2) does not account for the effects of waterlogging on transpiration, although the net effect on ET is likely limited because of the compensating effect of soil evaporation. In addition, under extremely dry conditions Eq.
(2) will overestimate the soil moisture stress, but such conditions approach the limits of the range for which the assumptions behind WALRUS are valid.

5
Data from the two catchments (Sect. 2) are used to estimate ζ 1 and ζ 2 (Fig. 6). The scatter in the observed evapotranspiration data is very large, but when data points are collected in 25 mm wide sliding bins and averaged, a decrease in β with d V can be observed (orange-red line). In the Cabauw polder, the storage deficit is never large and therefore hardly any evapotranspiration reduction occurs. In the Hupsel Brook catch-10 ment, reduction is around 10 % when d V exceeds 300 mm, which corresponds to a rare groundwater depth of about 2 m (about 14 % of the data in Fig. 6 was obtained during the extremely dry summer of 1976).
The open water evaporation is assumed to be equal to the potential evapotranspiration (ET pot ) of a well-watered soil. A Penman approximation would be more appro- 15 priate, but for most catchments only one estimate of evapotranspiration is available. In addition, the area fraction of open water and consequently the error is small. No evapotranspiration from the surface water occurs when the surface water reservoir is empty. Because the groundwater and surface water reservoirs together cover the entire catchment area, no evapotranspiration occurs from the quickflow reservoir.

Storage deficit
The dryness of the vadose zone is expressed by the storage deficit (d V ), representing the volume of empty soil pores per unit area, or in other words, the depth of water necessary to reach saturation. The vertical profile of soil moisture is not simulated explicitly and, as WALRUS is a lumped model, neither is its horizontal variability. The 25 storage deficit controls the precipitation division between groundwater and quickflow (W ), evapotranspiration reduction (β) and the change in groundwater depth (d G ) and is itself the result of all fluxes into or out of the soil reservoir, both the vadose zone and the groundwater zone.
In the field, time series of storage deficit (d V ) can be estimated from soil moisture (θ [-]) profile data. For each depth the soil moisture content at saturation (θ s [-]) has to be determined, which can often be done by taking the highest measured soil moisture 5 content at that depth. The difference between the profiles of θ and θ s gives the profile of the fraction of soil filled with air (and the remainder, 1 − θ s , gives the soil particle fraction). The storage deficit is obtained by integrating this air profile over depth d from the groundwater table d G to the soil surface: (3)

Equilibrium storage deficit
For every groundwater depth d G , an equilibrium soil moisture profile exists where at all depths gravity is balanced by capillary forces, and no flow occurs. From this profile the equilibrium storage deficit d V,eq can be derived in the same way as d V , namely by integrating the volume of empty soil pores over depth. The relation between d V,eq and 15 d G can be estimated from combined observations of groundwater and soil moisture. By assuming that on average d V,eq equals d V , the relation can be read from a (d G ,d V )-plot and supplied to WALRUS. Alternatively, one can assume a relation based on parametrisations of steady-state (i.e. no-flow) profiles reported by e.g. Brooks and Corey (1964) with b the pore size distribution parameter [-] and ψ ae the air entry pressure [mm].
The air entry pressure raises the power law distribution above the groundwater table to allow for the capillary fringe (the saturated area above the groundwater table). The 5 parameters b, ψ ae and θ s differ per soil type and selected results from laboratory experiments by Clapp and Hornberger (1978) are given in Table 2 (see Cosby et al., 1984, for interpolations between soil types). When the part of the profile between the capillary fringe and the soil surface from Eq. (4) is substituted in Eq. (3), the relation between equilibrium storage deficit and groundwater depth becomes Heterogeneities, such as soil layering or disruption by plant roots, macrofauna and human activity, cause differences between laboratory and field observations. In Fig. 7 d G is plotted as a function of d V for several sites in the Hupsel Brook catchment and 15 Cabauw polder area with corresponding theoretical curves. We computed the temporal maximum θ per depth (at the meteorological station in the Hupsel Brook catchment and the average of four profiles in the Cabauw polder) and averaged over the entire measured depth (205 cm in the Hupsel Brook catchment and 72 cm in the Cabauw polder) to obtain a single value of θ s . For the Hupsel Brook catchment, we fitted b 20 while retaining ψ ae , but for the Cabauw polder it was necessary to fit both b and ψ ae to obtain curves which describe the data points relatively well. The values obtained with these fits are listed in Table 2. Note that the data are actual storage deficits, which may not be in equilibrium with the groundwater depth measured at the same time. In addition, sites differ considerably and will deviate from the catchment average.

Percolation and capillary rise
In practice, the soil moisture profile and storage deficit are never perfectly in equilibrium with the groundwater depth. Addition (e.g. through precipitation) and removal (e.g. by drainage or evapotranspiration) of water cause an imbalance between gravity and capillary forces, leading to downward (percolation) or upward (capillary rise) flow to-5 wards a new equilibrium situation. Because the flow decreases with proximity to the equilibrium, this equilibrium will only be reached asymptotically.
The exact profile of relative saturation is not simulated explicitly in WALRUS, but the temporal dynamics of d V and d G caused by the interactions between groundwater and vadose zone are taken into account. The groundwater depth responds to changes 10 in storage deficit. The change in groundwater depth is parameterised as a function of the difference between the actual storage deficit (computed from the water budget in the soil reservoir) and the equilibrium storage deficit corresponding to the current groundwater level: 15 with c V the vadose zone relaxation time constant, which determines how quickly the system advances towards a new equilibrium. Four situations may occur (illustrated in Fig. 8).
(1) Water is added to the vadose zone through percolation. The actual storage deficit is smaller than the equilibrium for the current groundwater depth. Water will flow downward and the groundwater level 20 will rise gradually to the depth corresponding to the actual storage deficit. (2) Water is removed from the vadose zone through evapotranspiration. The actual storage deficit exceeds the equilibrium for the current groundwater depth. Water will flow upward to replenish the shortage in the top soil and the groundwater level will drop gradually.
(3) Water is removed from the soil reservoir though drainage, downward seepage or 25 groundwater extraction. Air is sucked into the soil and the actual storage deficit increases. This happens instantaneously, because water is incompressible. Water will 1375 Introduction percolate to reach an equilibrium profile again and the groundwater level will drop gradually. (4) Water is added to the soil reservoir through infiltration from surface water or upward seepage. The storage deficit decreases directly and the groundwater table rises gradually.

5
Drainage of groundwater towards the surface water reservoir or infiltration of surface water f GS is computed as with d G the depth of the groundwater table below the soil surface, c G a reservoir constant [mm h] and c D the average channel depth [mm] (see also Table 1 and Fig. 4). The 10 parameter c G represents the combined effect of all resistance and variability therein and depends on soil type (hydraulic conductivity) and drainage density. The first term of Eq. (7), c D − d G − h S , expresses the pressure difference driving the flow. The second term, max((c D − d G ), h S ), expresses the contact surface (parameterised as a depth) through which the flow takes place. These terms can be compared to the pressure 15 head difference and layer thickness commonly used in groundwater models. The contact surface-term accounts for decreasing drainage efficiency when groundwater and surface water levels drop and headwaters run dry. With this term, the variable source area concept (Beven and Kirkby, 1979) is implemented effectively and without additional parameters. 20 When groundwater drops below the surface water level, infiltration will be computed with the same relation, decreasing to zero when the surface water reservoir is empty (the second term max((c D − d G ), h S ) becomes zero). The same parameter c G is used for both groundwater drainage and surface water infiltration to limit the number of parameters, even though the resistance may be different in practice. The groundwater-surface water feedback is illustrated by a numerical experiment. We ran the model for an artificial 3 h rainfall event with an intensity of 10 mm h −1 with and without using h S in the drainage flux computation. Including h S leads to a decrease in drainage f GS and even infiltration (negative f GS ) during the peak (Fig. 9, left panels). This causes an attenuation of the discharge peak and higher groundwater levels after 5 the peak. This feedback is an important characteristic of WALRUS: in most parametric models, surface water levels are not modelled explicitly and this feedback cannot take place.

Quickflow
The quickflow reservoir simulates the combined effect of all water flowing through quick 10 flow paths towards the surface water: overland, macropore and drainpipe flow. This reservoir can therefore be seen as a collection of ponds, small drainage trenches or gulleys, soil cracks, animal burrows and drainpipes. Quickflow f QS depends linearly on the elevation of the water level in the quickflow reservoir h Q , with a time constant (reservoir constant) c Q : Water cannot flow from the surface water into the quickflow reservoir. Therefore, a sudden surface water level rise caused by an increase in surface water supply or weir elevation does not affect the quickflow reservoir directly. The water level in the quickflow reservoir cannot be coupled to measurable variables 20 directly -groundwater level measurements show the combined effect of the seasonal variation of the groundwater depth and the high resolution dynamics of the quickflow reservoir. Even though quickflow is parameterised as a single linear reservoir, it is essential to include this reservoir to mimic the large and variable contribution of these flowroutes (see Sect. 3.3).

Surface water
The surface water level h S represents the water level in the average channel with respect to the channel bottom. The distance between channel bottom and soil surface c D is calibrated or estimated from field observations. The stage-discharge relation Q = func(h S ) specifies the relation between surface water level and discharge at the 5 catchment outlet (in mm h −1 ). It is provided by the user as a function, e.g. the relation belonging to the weir at the catchment outlet, or as a lookup table. A threshold level h S,min can be included in the stage-discharge relation to account for a weir or other water management structures. If applicable, a value or time series of h S,min should be provided. When the surface water level drops below the crest of a weir, discharge will 10 be zero, but because there may still be drainage, infiltration and evaporation, it is important to include standing water. A default stage-discharge relation with the shape of a power law with a default exponent x S of 1.5 has been implemented: for h S ≤ c D . The default exponent value 1.5 for x S is inspired by equilibrium flow in 15 open channels (Manning, 1889). The parameter c S corresponds to the discharge at the catchment outlet (in mm h −1 ) when the surface water level reaches the soil surface, comparable to the bankfull discharge. It can be calibrated or provided based on field observations. 20 All fluxes across the catchment boundary, except for the discharge at the catchment outlet, are combined in the external groundwater flow term f XG (downward or upward seepage and lateral groundwater inflow or outflow) and the external surface water flow term f XS (supply or extraction Because these fluxes are added to the soil reservoir or surface water reservoir, they influence other variables through the different feedbacks implemented in the model. Most parametric rainfall-runoff models do not contain a surface water reservoir and therefore surface water supply can only be added to discharge afterwards and the impact of surface water increase on groundwater level and the groundwater drainage 5 flux is not considered.

Seepage and surface water supply
To investigate the effect of WALRUS' set-up considering surface water supply, we modelled an artificial event with two model set-ups: (1) f XS is added to the surface water reservoir and groundwater-surface water feedbacks are considered (as implemented in WALRUS) and (2) f XS is added to Q afterwards and h S is not used in the 10 groundwater drainage computation. Adding f XS to the surface water reservoir causes a gradual increase in h S and gradually rising Q (Fig. 9, right panels). When f XS is added to Q afterwards, h S is not affected by f XS and only increases after rainfall, and Q rises and falls instantly after changes in f XS . When a larger fraction of the catchment is covered by surface water (a S ), the increase in h S and Q becomes more gradual, because 15 the supplied surface water volume is spread out over a larger surface. Including the groundwater-surface water feedback leads to an attenuated discharge peak, caused by a decrease in drainage as a result of a decreasing difference between d G and h S . In dry periods, f XS may cause h S to rise above d G , leading to infiltration of surface water, which indicates that seepage and groundwater-surface water feedback should be 20 implemented together.

Large-scale ponding and flooding
The quickflow reservoir simulates the effect of local ponding and overland flow, but large-scale ponding may also occur. When the storage deficit becomes zero (i.e. all soil pores are filled with water), the groundwater level will rise directly to the surface 25 (as observed by Gillham, 1984;Brauer et al., 2011). Storage deficit and groundwater depth continue to drop (i.e. become more negative) together as there are no capillary forces any more and water level and pressure head coincide -negative d V and d G 1379 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | express ponding depths. Note that the levels rise less quickly above ground as the storativity becomes 1. Unfortunately, few quantitative, catchment-scale observations exist of different fluxes during floods. Because WALRUS has no spatial dimensions, the complex process of overland flow must be simplified. It is assumed that when the groundwater or surface 5 water level rises above the soil surface, the groundwater drainage/surface water infiltration flux f GS will include overland flow and is instantaneous, because overland flow is much faster than groundwater flow. When the surface water level exceeds the soil surface, discharge becomes less sensitive to changes in surface water level, represented by an abrupt change in the stage-discharge relation. However, as soon as the surface water level exceeds the soil surface, the excess water is led to the soil reservoir directly and therefore h S hardly rises above the soil surface. Therefore, we keep the same stage-discharge relation when h S > c D as a default. When the modelled groundwater table reaches the soil surface, an abrupt change in catchment discharge occurs. This is in contrast to the gradual activation of different flowpaths when the catchment 15 effective groundwater table is below surface (as represented by the wetness index).
We investigated the option of making the surface water area fraction a S a function of h S , representing gradual widening of brooks and inundation of areas close to the surface water network, and thereby smoothing the effect of flooding on discharge at the catchment outlet. Unfortunately, this approach made the model structure less intuitive 20 and introduced more degrees of freedom to define the shape of this function. Because flooding of the surface water reservoir only occurs during extremely wet situations, we chose to keep the model structure simple and leave a S fixed.

Outlook: possible model extensions
Some processes are not taken into account in the core model yet, but a user could 25 easily add preprocessing and postprocessing steps to adapt WALRUS to catchmentspecific situations. (1) The potential evapotranspiration estimated at a meteorological station may not be representative for the collection of vegetation types in the catchment.

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Therefore, one could use land cover distributions and crop factors to determine the catchment average potential evapotranspiration.
(2) Currently, WALRUS is set-up to receive liquid precipitation, but preprocessing steps to account for snow and/or interception can be added. For example, the delay in precipitation input caused by snow accumulation and melt can be simulated with methods based on the land surface en-5 ergy balance (Kustas et al., 1994) or a degree-day method (Seibert, 1997). (3) Interception can be parameterised with a threshold. Only the rainfall which exceeds the threshold is used as input for the model. The intercepted water evaporates directly and is not subtracted from ET pot (Teuling and Troch, 2005). (4) Paved surfaces have a low infiltration capacity, which limits groundwater recharge. This can be parametrized by 10 decreasing the groundwater reservoir area a G , introducing a paved surface area and leading the fraction of the rainfall belonging to this area directly to the surface water. (5) For large catchments, the discharge pulse from the model can be delayed and attenuated in the channels. It is possible to add a routing function to account for the delay and attenuation. 15 Another possibility is to couple WALRUS to other models. The outflow Q of one catchment can be used as surface water supply f XS for another WALRUS-unit downstream. With this technique, one could make a chain of WALRUS units to model subcatchments (with possibly different catchment characteristics and therefore parameter values) separately. Groundwater flow from one unit to the next can be computed from groundwater 20 levels in adjacent cells and Eq. (7). This groundwater flow is added to or subtracted from the seepage flux f XG for both units. Regional groundwater flow from a distributed groundwater model can be added to or subtracted from the soil reservoir through the seepage flux f XG . This can be specified with a time series or an external groundwater level. The outflow of the model can be used as input for a hydraulic model. Discharge

Model implementation
In this section we describe some key parts of the model implementation, which affect the model application and performance.

Code set-up
The model code is written in R, but can be easily translated into any vector-oriented

Initial conditions
The model can (as default) compute initial conditions for all states automatically, based on a stationary situation (thereby avoiding long burn-in periods). The quickflow reservoir is initially empty. The initial surface water level is derived from the first discharge observation and the stage-discharge relation. The initial groundwater depth is computed 20 with the assumption that initial groundwater drainage (f GS ) is equal to the initial discharge. It is also possible to supply the fraction of the initial discharge originating from drainage G frac and the model will solve for d G,0 with the quadratic formula and then use the remainder of the discharge to compute the initial quickflow reservoir level: Alternatively, the initial groundwater depth can be supplied (or calibrated) by the user and h Q,0 is computed such that Q 0 = f GS,0 + f QS,0 again. The initial storage deficit is 5 assumed to be the equilibrium value belonging to the initial groundwater depth.

Parameters
WALRUS has four parameters which require calibration: c W , c V , c G and c Q . These parameters have a physical meaning and can be explained qualitatively with catchment characteristics. The channel depth c D and surface water area fraction a S can be esti-10 mated from field observations. When the default stage-discharge relation is used, the bankfull discharge c S and (if applicable) the weir elevation h S,min need to be supplied (or calibrated) as well. Parameters are catchment-specific, but time-independent, to allow a calibrated model to be run for both long periods and events. We did not implement a specific calibration routine in the model, but used the HydroPSO package, 15 which is a particle swarm optimization technique (Zambrano-Bigarini and Rojas, 2013). The user can define the (multi-)objective function.

Forcing
Forcing data can be supplied as a time series or as a function (e.g. a sine function for ET pot or a Poisson rainfall generator). Observation times do not need to be equidistant, 20 which is especially useful for tipping-bucket rain gauges. Forcing time series are converted to functions (e.g. cumulative P as function of time), which allows other time steps than used for the original forcing.

If-statements
If-statements associated with thresholds cause nonlinearities in a model and their abrupt changes hamper calibration, in particular when using gradient-based methods. It is therefore important to know that there are four causes for abrupt changes in the model: (1) the stage-discharge relation (supplied by a user) may show abrupt changes 5 at the elevation of the crest of the weir or at the soil surface; (2) no evaporation occurs from empty channel beds; (3) if the storage deficit becomes negative or exceeds the groundwater depth, the groundwater depth becomes equal to the storage deficit; (4) if either groundwater or surface water level exceeds the soil surface, overland flow is instantaneous.

Integration scheme
The model is implemented as an explicit scheme, because nonlinearities caused by feedbacks and if-statements do not allow for the use of an implicit scheme. The states at the end of the previous time step are used to compute the fluxes during the current time step, which are then used to compute the states at the end of the current time 15 step (Fig. 10). The output data file lists the sums of the fluxes during and the states at the end of each time step.

Time step
The user can specify at which moments output should be generated, for example with a fixed interval (i.e. each hour or day), with increased frequency during certain events 20 or after each millimeter of rainfall. The output time steps can be both larger and smaller than those of the forcing. An important feature of the model code is the flexible computation time step. The model first attempts to run a whole output time step at once, but the time step is decreased when (1) the rainfall sum, discharge sum or change in discharge, surface water level or groundwater depth during the time step exceeds a certain threshold, or when (2) the surface water level is negative at the end of the time step. The first criterion prevents numerical instability caused by the explicit integration scheme and a delayed the response to rainfall as a result of the explicit model code (it takes one step to update the surface water level and another for the discharge). The second criterion is neces-5 sary because the total surface water outflow, computed from water levels at the start of the time step and the time step size, can exceed the available water. Because this means that non-existing water flowed out, there is a physical reason to avoid this. The procedure of decreasing time steps is illustrated in Fig. 10 (3rd step). First the original time step is halved and the model is run for this substep (of course with the 10 forcing corresponding to this substep). When the criteria are still not met, the step size will be halved again and again until the criteria are met. When one substep is completed, the fluxes are stored and the states at the end of the time step are used as initial values for the next substep. Then the model is run for the remainder of the original time step and, if necessary, the substep is halved until the criteria are met. This 15 will continue until the end of the intended output time step is reached. The sum of the fluxes of the substeps and the states of the last substep are stored in the output file.
The effect of the variable time step is illustrated in Fig. 11, in which the output of the model ran with a fixed time step and with variable time steps is shown. Note the erroneous time delay and magnitude of the discharge peak when no substeps are 20 used.

Water balance
WALRUS is a mass conserving model, and therefore the model water budget, computed as implies an increase in water in the reservoir. The groundwater level does not appear explicitly in the water balance, because it only plays a role as a pressure level driving groundwater drainage and surface water infiltration fluxes, while the storage deficit accounts for volume changes in the whole soil reservoir.

5
The Wageningen Lowland Runoff Simulator (WALRUS) is a new rainfall-runoff model, which is suitable for lowlands where shallow groundwater and surface water influence runoff generation. The model includes: 1. Groundwater-unsaturated zone coupling -WALRUS contains one soil reservoir, which is divided effectively by the (dynamic) groundwater table into a groundwater 10 zone and a vadose zone. The condition of this soil reservoir is described by two strongly dependent variables: the groundwater depth and the storage deficit (the effective thickness of empty pores). This implementation enables capillary rise when the top soil has dried through evapotranspiration.
2. Wetness-dependent flowroutes -the storage deficit determines the division of 15 rain water between the soil reservoir (slow routes: infiltration, percolation and groundwater flow) and a quickflow reservoir (quick routes: drainpipe, macropore and overland flow).
3. Groundwater-surface water feedbacks -surface water forms an explicit part of the model structure. Drainage depends on the difference between surface water 20 level and groundwater level (rather than groundwater level alone), allowing for feedbacks and infiltration of surface water into the soil.
4. Seepage and surface water supply -groundwater seepage and surface water supply or extraction (pumping) are added to or subtracted from the soil or surface water reservoir. These external fluxes affect the whole system through The open source model code is implemented in R and the model is set-up such that it can be used by both practitioners and researchers. For direct use by practitioners, defaults are implemented for relations between model variables and to compute initial 5 conditions, leaving only four parameters which require calibration. For research purposes, the defaults can easily be changed. WALRUS is computationally efficient, which allows operational forecasting and uncertainty estimation by creating ensembles. An approach for flexible time steps increases numerical stability and makes model parameter values independent of time step size, which facilitates use of the model with the 10 same parameter set for multi-year water balance studies as well as detailed analyses of individual flood peaks.
Numerical experiments shows that the implemented feedbacks have the desired effect on the system variables: (1) the wetness-dependent division between slow and quick flowroutes results in more quickflow, less recharge and higher discharge peaks 15 during wet periods; (2) the surface water level attenuates drainage during discharge peaks or when surface water is supplied upstream. An exhaustive test of WALRUS, with calibration, several validation studies, sensitivity analyses and uncertainty analyses in two catchments, the freely draining Hupsel Brook catchment and the controlled Cabauw polder, is presented in a subsequent paper (Brauer et al., 2014).

Code availability
The complete model code can be obtained upon request (by emailing the first author). In addition, the code will be made available shortly through the R CRAN website. WALRUS is licensed under the GPL v3 licence.

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version
States ] q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq qq q q q q q q q q qq qqq q q q q q q q q q q q q q q q q q q q q q q q qq q q q ] q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qqq qq qq q q q q q q q q q q q q q q q q q qq q q q q q q q q q q 30 50 e Chloride 10 Sep 1 Oct 20 Oct 1993   Fig. 4. Overview of the model structure with the five compartments: land surface (purple), vadose zone within the soil reservoir (yellow/red hatched), groundwater zone within the soil reservoir (orange), quickflow reservoir (green) and surface water reservoir (blue). Fluxes are black arrows, model parameters brown diamonds and states in the colour of the reservoir they belong to. For a complete description of all variables, see Table 1 and Sec. 4.1. The names of the fluxes are derived from the reservoirs (for example f XS : f stands for flow, the X for external and the S for surface water -water flowing from outside the catchment into the surface water network). 44 Fig. 4. Overview of the model structure with the five compartments: land surface (purple), vadose zone within the soil reservoir (yellow/red hatched), groundwater zone within the soil reservoir (orange), quickflow reservoir (green) and surface water reservoir (blue). Fluxes are black arrows, model parameters brown diamonds and states in the colour of the reservoir they belong to. For a complete description of all variables, see Table 1 and Sec. 4.1. The names of the fluxes are derived from the reservoirs (for example f XS : f stands for flow, the X for external and the S for surface water -water flowing from outside the catchment into the surface water network).    (Wand and Jones, 1995 ig. 10. Illustration of the variable time step procedure. (a) The non-equidistant output time teps (purple) are used as first attempts for computations of fluxes (blue/green) and states urple), but during time step number 3, the precipitation sum is too large (panel b) and the tep is divided into substeps: it is halved and then halved again until the criterion was reached. ote that even though the size of output time step 2 is larger, it is not divided into substeps, ecause all criteria are met. The non-equidistant output time steps (purple) are used as first attempts for computations of fluxes (blue/green) and states (purple), but during time step number 3, the precipitation sum is too large (panel b) and the step is divided into substeps: it is halved and then halved again until the criterion was reached. Note that even though the size of output time step 2 is larger, it is not divided into substeps, because all criteria are met. ig. 11. The effect of variable time steps on the model output. An artificial case with a rainfall ent of 30 mm in the first hour and no evapotranspiration. The lines connect the discharge odelled at the end of a time step (instantaneous value), and do not represent the sum over e time step (which is given in the output file). The same parameter values as in Fig. 5 were ed. Fig. 11. The effect of variable time steps on the model output. An artificial case with a rainfall event of 30 mm in the first hour and no evapotranspiration. The lines connect the discharge modelled at the end of a time step (instantaneous value), and do not represent the sum over the time step (which is given in the output file). The same parameter values as in Fig. 5 were used.