Catchment-scale non-linear groundwater-surface water interactions in densely drained lowland catchments

Freely discharging lowland catchments are characterized by a strongly seasonal contracting and expanding system of discharging streams and ditches. Due to this rapidly changing active channel network, discharge and solute transport cannot be modeled by a single characteristic travel path, travel time distribution, unit hydrograph, or linear reservoir. We propose a systematic spatial averaging approach to derive catchment-scale storage and discharge from point-scale water balances. The effects of spatial heterogeneity in soil properties, vegetation, and drainage network are lumped and described by a relation between groundwater storage and the spatial probability distribution of groundwater depths with measurable parameters. The model describes how, in lowland catchments, the catchment-scale flux from groundwater to surface water via various flow routes is affected by a changing active channel network, the unsaturated zone and surface ponding. We used observations of groundwater levels and catchment discharge of a 6.6 km 2 Dutch watershed in combination with a high-resolution spatially distributed hydrological model to test the model approach. Good results were obtained when modeling hourly discharges for a period of eight years. The validity of the underlying assumptions still needs to be tested under different conditions and for catchments of various sizes. Nevertheless, at this stage the model can already improve monitoring efficiency of groundwater-surface water interactions. Correspondence to: Y. van der Velde (ype.vandervelde@wur.nl)


Introduction
Catchments without real hillslopes, with an unconsolidated soil, a dense artificial drainage system, and with high inputs of nutrients due to intensive agriculture can be found in lowland landscapes all over the world.Polluted surface waters are an important environmental issue in all these catchments, with nutrient loads far exceeding loads in most mountainous catchments.Recent research on catchment scale discharge and transport modeling, however, was mainly oriented towards sloped catchments, creating concepts and models that are inappropriate for lowland catchments e.g.TOPMODEL by Beven and Kirkby (1979); ARNO by Todini (1996); Representative Elementary Watershed (REW) approach as implemented by Zhang et al. (2006); HBV by Lindström et al. (1997).
Typically, lowland catchments have a soil with sand, clay, and peat layers, sometimes interspersed with gravelly layers, with a shallow groundwater table.The absence of significant slopes makes groundwater the dominant contributor to stream discharge, either via direct inflow through the stream bed or through man-made drainage systems (De Vries, 1994;Wriedt et al., 2007;Tiemeyer et al., 2007).This groundwater flux is driven by continuously changing groundwater level gradients towards draining ditches and streams rather than by a fixed regional bedrock or surface elevation slope as is a common assumption for sloped catchments.Direct runoff occurs only when the infiltration capacity of the soil is exceeded by heavy rainfall or when the phreatic level rises to the soil surface.Freely discharging lowland catchments are characterized by a strongly seasonal contracting and expanding system of discharging ditches and streams (Ernst, Y. van der Velde et al.: Catchment-scale non-linear groundwater-surface water interactions 1978;De Vries, 1995).In hillslope hydrology this changing active channel network is reflected in the hydrological connectivity (Ocampo et al., 2006;Molenat et al., 2008) between the riparian and upland zones.Due to this rapidly changing active channel network, discharge and solute transport cannot be modeled by a single characteristic travel path, travel time distribution, unit hydrograph, or linear reservoir.This highly non-linear, transient behavior is well recognized (Beven and Kirkby, 1979;Van de Griend et al., 2002).Many approaches incorporated a variable contributing area concept for the description of stream discharge, but most of them focused on hillslopes (TOPMODEL based on a kinematic wave approach, Beven and Kirkby, 1979), direct runoff (PDM rainfall-runoff model based on spatial distribution of soil moisture, Moore, 1985), or characteristic soil-segments (Lazzarotto et al., 2006), making them over-parameterized and needlessly complicated for applications to large catchments, or even irrelevant for relatively flat lowland areas.Moore (1985Moore ( , 2007) ) proposed a probability distribution for soil moisture storage to include the spatial variability of discharge generation, but did not relate this to a distribution in groundwater levels.Discharge generation in lowland catchments, however, is driven to a far greater extent by the distribution of groundwater levels than it is by the soil moisture content of the top layer.Seibert et al. (2003) explored this interaction between groundwater level and unsaturated soil moisture and concluded that runoff models for catchments with shallow groundwater levels should explicitly include unsaturated zone storage coupled to groundwater levels.Wriedt et al. (2007), Ocampo et al. (2006) and Molenat et al. (2008) showed that hydrological connectivity through channel activity or high groundwater tables can be one of the major controls of nitrate transport within a catchment.A spatially distributed hydrological model can in principle calculate these spatial and temporal groundwater dynamics but has a huge data demand and to model correct contributions of specific flow routes (overland flow, tube drain flow, or groundwater flow) to the total discharge, very small spatial and temporal resolutions would be needed.This causes long building and calculation times and makes such models tedious to operate and calibrate.Rainfall-runoff models with variable source area concepts, on the other hand, can effectively calculate fluxes of individual flow routes when measurements are available, but their storage volumes are often inaccurate.Both aspects, an accurate separation in flow route contributions and accurate storage volumes, are essential for catchment-scale solute transport modeling.Molenat et al. (2007), Ocampo et al. (2006) and McDonnell (2003) reached a similar conclusion and suggested that for an accurate description of nitrate transport a classic "variable source area" model is not the way forward.
The objectives of this paper are to formulate expressions for catchment-scale water fluxes from the unsaturated zone to the groundwater and from groundwater to the stream network.The expressions need to incorporate spatial and tem-poral groundwater variations and should calculate realistic storage changes within the catchment.We apply these equations to a lowland agricultural catchment in The Netherlands (Hupsel Brook catchment, 6.6 km 2 ) and evaluate their performance.

The basics
We seek to develop a water flux model for densely drained lowland catchments without snow cover.The model should be able to describe the dynamic saturated groundwatersurface water contact interface.The interaction between the saturated and unsaturated zone is expected to help generate peak discharges during wet periods by amplifying the precipitation signal toward the saturated zone (Seibert et al., 2003).Surface ponding and water storage by filling dry ditches and stream branches, on the other hand, is expected to dampen peak discharges during wet periods.Both types of interactions are included in the model description.
In lowland catchments groundwater discharge from the saturated zone to the surface water system is the most important discharge generating process.It can occur as flow into tube drains (q dr [LT −1 ]), flow into ditch and stream beds, and as overland flow from groundwater seepage when the phreatic level is above the soil surface.Both overland flow and groundwater seepage into ditches and streams occur because groundwater levels rise above the level of the water layer on the soil surface (which may also be the stream/ditch bed) and therefore they both received the notation: q ex [LT −1 ].Discharge generation is generally described by a linear reservoir with a threshold, driven by groundwater heads, H (x,y,t): A location x, y [L] at time t [T] starts to generate discharge, q i (x,y,t) [LT −1 ], when the groundwater head, H (x, y, t) [L], is larger than a threshold groundwater head, H thres,i (x,y,t) [L].The resistance that this water flux has to overcome is denoted by r i (x,y,t) [T].The subscript i denotes the type of flux (groundwater flow towards surface water and surface ponds: i=ex and tube drain flow: i=dr).Since we limit ourselves to groundwater discharging directly into the surface water, the discharge flux can only be non-zero along the wet perimeters of stream beds, along the tube drains below groundwater level and at the soil surface when the phreatic level reaches the surface and overland flow occurs.
For stream/ditch and overland flow (q ex ), H thres,ex is the surface water level, or the soil surface elevation when there is Hydrol.Earth Syst.Sci., 13, 1867-1885, 2009 www.hydrol-earth-syst-sci.net/13/1867/2009/ no water storage on the soil surface.For tube drain discharge (q dr ), H thres,dr , is the elevation of the drain tube.A catchment can be viewed as a population of such point-scale linear reservoirs with individual values for H , H thres,i and r i .The draining area A q,i [L 2 ], i.e. the area of the catchment where groundwater and surface water are in direct contact, is then defined by: with 1 {var} an indicator function that is 1 when variable var is true and 0 when var is false and A [L 2 ] the catchment area.The values of H (x,y,t), the groundwater level, and H thres,i (x,y,t), the surface water level, are strongly time dependent and may cause the drainage area A q,i to vary strongly in time.In relatively flat lowland catchments with dynamic and shallow groundwater levels, A q,i has been observed to change considerably over time (Ernst, 1978;Wriedt et al., 2007;De Vries, 1995).This is a combined effect of groundwater tables that lose contact with surface water or tube drains during dry periods (compare the wet and dry state in Fig. 1) and of high surface water levels during wet periods, raising the threshold groundwater head.Consequently, models that use one linear reservoir to calculate groundwater flow towards the surface water network, which rely on the assumption that A q,i is constant with time, fail to describe groundwater discharge in lowland catchments or need multiple reservoirs to model discharge.Often, a fast-and slowresponse reservoir arranged in parallel are used.Although conceptually straightforward, this modeling strategy does not fully recognize the system dynamics and its parameters cannot be directly linked to observable catchment properties.
In lowland catchments a huge simplification can be made in upscaling Eqs. ( 1) and (2) when the change in saturated groundwater storage related to a change in groundwater level is expressed by a change in the thickness of the unsaturated zone, u[L].It is important to realize that from hereon, we will use the change in unsaturated zone thickness to express the change in saturated storage.In lowland catchments with a shallow phreatic groundwater system, the spatial variation of u(x,y,t) is heavily affected by the distance between draining ditches or tube drains (see how the groundwater level during a wet day is affected by ditches and drains in Fig. 1).This yields a spatial distribution between draining ditches or drains of point-scale values of u(x,y,t) that is mainly influenced by soil type, drainage depth and distance, and recharge flux.A lowland catchment typically has a dense network of ditches and drains with many different drainage depths and distances between ditches and drains.Thus the spatial distribution of u(x,y,t) at any given time over the entire catchment is the sum of the spatial distributions of u(x,y,t) at that time between actively draining ditches and drains.According to the central limit theorem, summing n distributions of weakly correlated random variables with finite means and variances, will yield a normal overall distribution for sufficiently large n (Feller, 1971).The key characteristic of our model is that the distribution of point-scale u(x,y,t) for the entire catchment is described by a Normal distribution function, f u (u(t), u(t) , σ u (t)) with mean unsaturated zone thickness, u [L], and standard deviation, σ u [L].From hereon the Normal distribution will be denoted to by f u (t), reflecting that each time has a unique spatial distribution of unsaturated zone thicknesses.The validity of this Normality assumption will be assessed in the Results and Discussion section.The locations with negative values for u(x,y,t) described by f u (t) indicate locations with a seepage face (i.e.groundwater is higher than the soil surface).This negative fraction of the distribution will be used to calculate the exfiltration fluxes of groundwater to the surface water (q ex ).The spatial structure of u within the catchment is lost, but the mean and variance of the values of u are preserved.Hence, no information on the location of a flux is available and consequently water cannot be routed downstream within the catchment.The model requires that the catchment characteristics are statistically homogeneous so that all local distributions of u have a mean and variance within the same order of magnitude and that the local distributions are to some degree independent.Therefore, it is not possible to choose a catchment size larger than typical rainfall and potential evaporation patterns, or to have significant trends or discontinuities in stream network densities or soil properties within the catchment.However it is Y. van der Velde et al.: Catchment-scale non-linear groundwater-surface water interactions possible to couple multiple models to account for these spatial discontinuities.These are the preliminaries from which the model is developed below.

Mass balance equation
The basis of the model is the mass balance equation for the saturated zone, the unsaturated zone and surface storage for each vertical column in the landscape (no changes in water density are assumed): ∂s surf (x,y,t) ∂t + ∂s unsat (x,y,t) ∂t + ∂s sat (x,y,t) ∂t = p(x, y, t) − e act (x, y, t) − l sat (x, y, t) With s[L] reflecting storage within a vertical column located at horizontal coordinates x, y, at time t.The subscripts surf, unsat and sat refer to storage of surface water/ponds, unsaturated soil water and saturated groundwater respectively.Rainfall is denoted by p[LT −1 ] and evapotranspiration by  4) requires the fluxes between these compartments.The fluxes between the unsaturated and the saturated soil are denoted by j [LT −1 ] while q[LT −1 ] denotes the fluxes from soil to the surface storage and vice versa: ∂s surf (x,y,t) ∂t = 1 {ssurf (x,y,t)>0} (p(x, y, t) − e act (x, y, t)) −q inf (x, y, t) + q ex (x, y, t) + q dr (x, y, t) −l surf (x, y, t) (5) ∂s unsat (x, y, t) ∂t =1 {s unsat (x,y,t)>0} (p(x, y, t)−e act (x, y, t)) ∂s sat (x,y,t) ∂t = j rch (x, y, t) − j cap (x, y, t) + q inf (x, y, t) −q ex (x, y, t) − q dr (x, y, t) −l sat (x, y, t) − o(x, y, t) Subscripts of q denote the infiltration from surface storage into the unsaturated zone, inf, exfiltration of groundwater to the surface water and surface ponds, ex, and groundwater flow towards tube drains, dr.Subscripts of j denote capillary up rise of groundwater to the unsaturated zone, cap, and the recharge of the saturated zone by unsaturated soil water, rch.Note that the flux into the drains appears in the surface water budget (Eq.5).Although counterintuitive, it signals that tube drain discharge no longer flows through the porous medium.Similarly, l surf comprises lateral fluxes of water both over the land surface, and through drain tubes.Both q dr and the tube drain contribution to l surf can only be nonzero for (x,y) located directly above a drain tube.Note that we assume that perched water tables do not occur.Therefore, one of the storages s surf or s unsat is necessarily zero and consequently the atmospherical forcings, p and e act , act on the active reservoir (s surf >0 or s unsat >0).All subsurface flows towards drains and surface water bodies are incorporated in l sat and all overland flows towards the surface water and flow from adjacent streams, ditches, and drains are incorporated in l surf .Figure 2 summarizes all fluxes that are described by this model.Equation ( 4) represents a point-scale mass balance.By integrating over the catchment area A [L 2 ], a catchment-scale mass balance can be obtained.In doing so, lateral flow components within A cancel out, and only the lateral flow over the boundary of A affects the mass balance.Thus we obtain: .We dropped the reference to the spatial and temporal coordinates for clarity.Of particular interest is the last term of Eq. ( 8) because this term represents the total catchment discharge by surface water at any given time.

Dimension reduction of the catchment scale mass balance equation
The integral formulation of the mass balance, Eq. ( 8), has two spatial dimensions and one time dimension, and generally will be impossible to evaluate in a practical way.We therefore seek a dimensional reduction approach in which we lump spatially distributed processes where possible while maintaining the characteristic behavior of a typical lowland catchment with realistic water storage changes inside the catchment.The characteristic behavior we focus on is defined by: -A continuously changing active drainage system defined by the contact zone between saturated groundwater and surface water, due to varying groundwater and surface water levels (Fig. 3a, b, c, and d).
-The unsaturated zone as an amplifier of rainfall and evapotranspiration fluxes towards and from the groundwater.This figure illustrates three locations (x 1 ,y 1 ), (x 2 ,y 2 ) and (x 3 ,y 3 ) within a cross section of a typical lowland field.The groundwater level at location (x 1 ,y 1 ) is above soil surface, which leads to ponding.Note that when the groundwater level is above soil surface there is no unsaturated zone.Infiltrating water from the pond into the saturated zone is denoted q inf .Exfiltrating water from the saturated groundwater into the pond is denoted q ex .A sink is denoted o and the lateral overland flow l surf .Location (x 2 ,y 2 ) has an unsaturated zone and consequently no surface storage.The flux from the unsaturated zone to the saturated zone is denoted j rch and the capillary flux from saturated to unsaturated zone j cap .This location is also tube-drained with a tube drain flux q dr .Note that surface storage and tube drainage can occur at the same location.Point (x 3 ,y 3 ) is located at a stream.Above the stream bed surface storage occurs.The exfiltrating, infiltrating, and lateral surface fluxes are treated the same way for a ponded location (x 1 ,y 1 ) and a stream/ditch location(x 3 ,y 3 ).Rainfall, p, evapotranspiration, e act , and lateral saturated groundwater fluxes, l sat , occur in all three locations.
-Ponding of parts of the soil during prolonged periods of rain (Fig. 3e, f).
As a first step, we eliminate the spatial dimensions in Eq. ( 8) by spatial averaging.Spatial averaging is simply obtained by carrying out the integration over A for that variable and dividing by A. Thus we obtain: Where denotes the spatial averaging operation over any A.
Note that the dimensional reduction changed the dimensions of all terms from [L 3 T −1 ] in Eq. ( 8) to [LT −1 ] in Eq. ( 9).When we choose the catchment such that its boundaries are zero-flux boundaries for the shallow groundwater, the boundary integral of the saturated lateral flux can be neglected.Even in the case of a large-scale background flow of groundwater passing through the catchment the net flux over S will be close to zero if no significant groundwater exfiltration or recharge of the aquifer occurs.The boundary integral of lateral fluxes of surface storage on the other hand, represents the total stream discharge from the catchment.This is of course the key flux that can be compared with discharge measurements.
The storage and flux terms in Eq. ( 9), are functionally dependent on the thickness of the unsaturated zone: low phreatic levels lowers e act (x, y, t), q ex (x, y, t) and q dr (x, y, t).In soils with a high infiltration capacity, overland flow, l surf (x, y, t), will be zero if u(x, y, t) is significantly larger than zero.We formalize this by declaring all local flux densities dependent upon u(x, y, t): Were J [LT −1 ] denotes a flux density or change in storage in Eq. ( 9), and g J () denotes a non-linear functional dependence on the variables in parentheses.The spatial average of J is: Y. van der Velde et al.: Catchment-scale non-linear groundwater-surface water interactions Note that the spatial dependence is replaced by a dependence of J on u through the probability density function (PDF) of u at the time of interest, which describes the spatial variation of u.Equations ( 10) and ( 11) reduce the problem of the spatial variation of the many terms in Eq. ( 9) to that of the variation in u and identifying g J () for the various J 's.
If we assume f u (t) to be Normal as discussed above, f u (t) is completely characterized by its mean u(t) and standard deviation σ u (t).By noting that during dry periods e act (x,y,t) tends to be large for small u(x,y,t), we can deduce that σ u (t) is relatively small during prolonged dry periods: shallow groundwater levels are lowered more than deep groundwater tables, reducing the variation of u(x,y,t) for large u(t) .During and shortly after rainfall, with ditches and drains discharging, u varies strongly within fields, increasing σ u (see also the cross section of Fig. 1).For prolonged rainfall, u(t) will reduce further, and the occurrence of ponding creates negative values of u.Eventually, when nearly the entire catchment is flooded, the water level above the soil surface will run approximately parallel to the groundwater level under dry conditions.Consequently, it is expected that σ u will tend towards the same relatively low value under very wet and very dry conditions.Based on these arguments and in the spirit of dimension reduction we will consider σ u (t) to be a function of u(t) that peaks at an intermediate value and tails off at the extremes.The exact functional dependence is a characteristic of the catchment topography, soil, and climate.A simple empirical four-parameter expression to approximate this relation is given by: where σ max [L] is the maximum standard deviation of u, occurring at u(t) =u sdmax [L].The minimum standard deviation, σ min [L], occurs for large and very small (negative) u(t) values.The shape parameter b[L] determines the steepness of the curve.The ability of this empirical function to describe the complex shape of the catchment-scale groundwater table will be assessed in the Results and Discussion section.

Storage and flux expressions
In this section the terms of the water balance, Eq. ( 9), are one by one expressed as functions of u and f u .Section 2.5 gives the final water balance equation, which is used to calculate catchment-scale fluxes and storages.

Temporal variations of average saturated storage
The point-scale saturated storage, s sat , is defined as: θ s (x, y, z)dz for u(x, y, t)>0 (13) where z[L] is the vertical coordinate, z 0 [L] is the elevation of the impermeable base or another suitable lower boundary, z s [L] is the elevation of the soil surface, and θ s is the saturated volumetric water content.Since we are interested in storage of water at a given horizontal location, the exact vertical location is of limited value.By noting that z s (x,y)−z 0 (x,y) is the local thickness T [L] of the subsurface affecting the catchment hydrological behavior, Eqs. ( 13) and ( 14) can be simplified to: Where z * is a transformed coordinate defined as z * =z−z 0 (x,y).
The first term on the right-hand-side of Eq. ( 15) is a location-specific constant, if temporal variations in θ s caused by soil tillage, biological activity etc. are neglected.It reflects the total pore space [L] in the column at (x,y).Likewise, the second term represents the total pore space in the unsaturated zone at (x,y).Averaging Eq. ( 15) gives the average groundwater storage of the catchment as the difference between the total and the unsaturated volumes of pores in the catchment: where we dropped the references to the spatial coordinates for clarity.The change of the average saturated storage is: The time derivative of θ s dz * is determined by the depth interval in the soil between the maximum and the minimum value of u(x,y,t).If θ s varies little within that interval, Eq. ( 17) simplifies to: If θ s and ∂ ∂t u(t) are uncorrelated random variables distributed over A, the average of their product equals the product of their averages: Hydrol.Earth Syst.Sci., 13, 1867-1885, 2009 www.hydrol-earth-syst-sci.net/13/1867/2009/ where θ s is evaluated between the highest and the lowest groundwater level.Applying Eq. ( 11) for positive values of u yields:

Temporal variations of average unsaturated storage
The unsaturated zone is assumed to be in hydrostatic equilibrium with the groundwater table at all times, making the volume of stored water in the unsaturated zone a function of the soil type and the water table.This assumption is only valid for shallow ground water tables, but has proven to be very useful in estimating the total amount of water in the unsaturated zone and its effect on groundwater table fluctuations (Kim et al.,1996;Bierkens, 1998).The equilibrium assumption implies that any water added to the soil (e.g. by precipitation) is transferred immediately to the groundwater.Similarly, any water removed from the unsaturated zone (e.g. by evapotranspiration) is immediately withdrawn from the groundwater.
The assumption of instantaneous equilibrium throughout the unsaturated zone implies that the soils will always be on the wet end of the soil water characteristic.We therefore use van Genuchten's (1980) expression with the dry-end residual water content equal to zero: Where α[L −1 ] and n are location-specific shape parameters and h(t)=z−z s (x,y)+u(x,y,t) is the height above the phreatic water level [L].The point-scale unsaturated zone storage, s unsat , can be obtained by integrating Eq. ( 21) for z ranging from z s −u to z s .Similarly to Eq. ( 11), the spatial average can be obtained by integrating across all positive values of u, where θ s represents the spatial average of the local vertically integrated θ s of the unsaturated zone, already introduced in Eq. ( 20).At catchment-scale, however, we do not define spatial average Van Genuchten parameters, α, and n, but we view them as effective parameters describing the storage behavior of the unsaturated zone of the catchment incorporating the effects of unsaturated zone heterogeneities.
The temporal derivative follows directly.Note that the assumption of instantaneous hydrostatic equilibrium of the unsaturated zone implies that ∂s sat (t)   ∂t and ∂s unsat (t)   ∂t in Eq. ( 9) have opposite signs and the absolute value of ∂s sat (t)   ∂t is always the largest (if the average thickness of the unsaturated zone increases, the saturated storage decreases, and the storage of the unsaturated zone increases).Effectively, the unsaturated zone amplifies the effects of the atmospheric fluxes on the groundwater table.

Temporal variations of average surface storage
Storage on the soil surface is assumed to occur only when groundwater levels rise above the soil surface.Ponding due to high rainfall intensities is assumed not to occur, which is valid for permeable soils in climates without long highintensity rainfalls.A linear relation is assumed between the surface storage depth, s surf [L], and the height of the groundwater level above soil surface at location (x,y) (i.e.negative values of u): where m[−] is a location-specific empirical constant with a value between 0 and 1 that gives the fraction of the excess water stored on the soil.If m=1, the negative u is entirely accounted for by the depth of the water layer on the soil surface.Consequently, no water is removed from the location by overland flow.For m<1, a water layer of thickness −m•u is stored on the soil surface, and the pressure head difference (m−1)•u generates overland flow.For m=0, no ponding occurs and all excess water is discharged.This relation underestimates the complexity of the generation of overland flow and groundwater flow towards surface water at the point-scale but it is expected that the averaging operation over the catchment, with its wide range of negative u values, gives a reasonable approximation of increased surface storage with decreasing average unsaturated zone thickness.Assuming independence between the factor m and u and applying Eq. (11) gives: The temporal derivative follows directly.Note that the time derivative, , has the same sign as, and is always smaller than ∂s sat (t) ∂t in Eq. ( 9): when the thickness of the unsaturated zone decreases, the saturated storage and the surface storage increase (with a thinner average unsaturated zone, there will be more ponding and therefore a higher surface storage).This term dampens the fluctuations in groundwater levels needed to maintain the water balance Eq. ( 9) and consequently dampens peak discharges.
Each negative thickness of the unsaturated zone translates into a fixed volume of stored water on the surface.This assumption implies that lateral surface fluxes cannot be stored elsewhere in the catchment (all available surface storage is always occupied) and that consequently surface water The catchment scale discharge can be calculated from the mass balance equation of the surface storage reservoir, Eq. ( 5): From the assumptions of instantaneous equilibrium of the surface storage reservoir, we can define q grw (t) , the groundwater exfiltration additional to the water needed to fill the surface storage, as: Note that the average infiltration flux density, q inf (t) , is zero when surface storage increases, i.e.

Groundwater exfiltration
Exfiltration of groundwater, q grw (x,y,t) [LT −1 ], defined by Eq. ( 27), is assumed to occur only when a groundwater head is higher than the level of the water layer stored on the soil surface.We also assume that groundwater exfiltration is proportional to the magnitude of the difference between the groundwater level u(x,y,t) and the surface storage level, s surf (x,y,t), yielding: With r grw (x, y, t) [T] the resistance that the water flux from soil to surface water must overcome.Replacing r grw (x, y, t) by its catchment scale average r grw , invoking Eq. ( 11) and introducing Eq. ( 23) gives the catchment-scale average groundwater exfiltration rate:

Tube drain discharge
Tube drain discharge occurs when the drainage depth, d dr (x,y) [L], is larger than u(x,y,t): with the function g(x, y)[−] equal to one above a drain tube and equal to zero elsewhere, and r dr (x, y, t)[T] denoting the resistance that the water flux from soil to tube drain has to overcome.However, drainage fluxes derived only at the exact location of drain tubes are of little practical value.We therefore introduce q * dr [LT −1 ] as the rate at which saturated flow towards nearby drain tubes removes water from a location (x,y) at time t.Consequently, this fraction of the total flow should be subtracted from the value of l sat (x,y,t) to maintain mass conservation.We then have: Where g * (x,y) equals one whenever (x,y) is in a tube drained field and is zero elsewhere.The drainage depth d * dr [L] gives the average drainage depth of the field in which (x,y) is located.Similarly r * dr [T] denotes the resistance to the flow towards and into the drain tube.When g * =0, d * dr and r * dr are undefined.In order to express q * dr as a function of u(t), we assume u(t) and d * dr to be independent.For the drained area of the catchment we may then write: Where the averaging operations have been carried out over the region within A where g * =1.Some of the very wet locations within a catchment (small u) are likely not to be drained.For example there are no drains under ditches and streams which are obviously the wettest locations in the catchment.For an accurate contribution of tube drain discharge to the total discharge under dry conditions it is important to define this fraction of the catchment (wet and undrained).When we would ignore this and assume drainage to be more or less uniformly distributed over the full range of u, the model will generate substantial tube drain discharge even under dry conditions.We therefore assume that a fixed fraction of the catchment area (A nd,wet Hence: Note that of course many of the undrained fields simply are dry enough without drain tubes.Therefore A nd,wet is smaller than the total undrained area.Equation ( 35) constitutes an additional condition that must be satisfied for q * dr to be nonzero.Extending Eq. ( 33) accordingly yields: Again, we determine the catchment average drainage discharge flux density by applying Eq. ( 11), taking into account that only the drained area A g * dA generates discharge:

Rainfall and evapotranspiration
Rainfall does not depend on u.We assume the catchment small enough for the rainfall rate p(x,y,t) to be uniform: p(t).Thus, p(t) =p(t).
In soils with shallow groundwater and a humid climate, transpiration by far exceeds evaporation when the plant cover is complete.In autumn and winter, cropped soils are bare, but the evapotranspiration rate in this period is low.The transpiration is assumed equal to the potential evapotranspiration, e pot [LT −1 ], as long as u(x,y,t) is smaller than some threshold.When u(x,y,t) exceeds that threshold, e act (x,y,t) drops to zero.It is expected that the averaging operation over the catchment with its wide range of local values of u produces a smoothly decreasing e act (t) as the catchment becomes drier.For a threshold u et (x,y,t) [L] we have: e act (x, y, t) = 1 {u(x,y,t)<u et (x,y,t)} • e pot (x, y, t) Applying Eq. ( 11) with u et (x,y,t) constant in time and space gives the average transpiration rate over the catchment: A more elaborate function such as a linear or exponential decline between two groundwater depths or a linear decline with unsaturated water content will only improve the results when the standard deviation of groundwater depth is small (<0.2).The averaging effect of the catchment will then be less, and only then the effect of the extra parameters of a more elaborate function will not be overruled by the averaging effect.For the entire Hupsel brook catchment we have chosen the most basic formulation as presented above.

The water balance as function of groundwater table fluctuations
In the previous sections all terms of the water balance, Eq. ( 9), have been made solely dependent on u(t) and f u (t).We now take Eq. ( 9) and substitute Eqs. ( 20), ( 22), and (24) for the three storage terms, maintain the precipitation term p(t) , and set the source/sink term o(t) to zero, and assume the net subsurface flux l sat across S to be negligible.Finally we insert Eq. ( 39) for the evapotranspiration, and Eq. ( 25) for the net surface water flux across S to obtain the water balance of the catchment: with the total discharge from the catchment, l surf (t) , derived from Eq. ( 28) combined with expressions for the individual flux terms, Eqs. ( 30), ( 37) and ( 39): where we assume zero travel time in the surface water.Note that e act is equal to e pot for negative values of u.Hence, e act in Eq. ( 28) is replaced by e pot .When we combine these two equations with a relation between u(t) and σ u (t) as Y. van der Velde et al.: Catchment-scale non-linear groundwater-surface water interactions given by Eq. ( 12), the model is complete.The advantage of the presented probability distribution function approach is that all point-scale threshold values for which a flux generating process is (de)activated have been translated into gradual changes and smooth transitions between fluxes at the catchment-scale, without introducing many new parameters.Therefore this model is stable in backwards iterations and during automatic calibration.
In this model, changes in saturated storage drive all catchment fluxes.The saturated storage change is dictated by the relation between mean and standard deviation of a Normally distributed thickness of the unsaturated zone.However, this relation cannot be derived by measuring catchment discharge only.When we want to apply this model, we need to derive this relation separately.Fortunately, it is possible to measure the spatial distribution of groundwater depth (=thickness of unsaturated zone) by measuring many randomly located groundwater depths or to use a spatially distributed groundwater model to derive the spatial distribution of groundwater depths.The latter method is less accurate because errors in the groundwater model propagate to the water balance model.
Other models such at TOPMODEL (Beven and Kirkby, 1979), the soil routine in HBV (Lindström et al., 1997) and the PDM rainfall runoff model (Moore, 1985) also use spatial distributions.These models have chosen slope type, soil type or soil moisture storage, of which the distributions remain constant in time, as the primary source of spatial variation.Because we deal with lowland catchments, the spatial distribution of groundwater levels drives discharge generation.This spatial distribution of groundwater depth, however, is not a constant in time but a function of storage.We defined relations between the distribution parameters and the storage.This resulted in a much more dynamical model driven by continuously changing groundwater head gradients.

Catchment characteristics
The Hupsel Brook catchment is located in the eastern part of The Netherlands (Fig. 4).The size of the catchment is about 6.6 km 2 , with the surface elevation ranging from 22 to 30 m above sea level.The soil texture class is mostly loamy sand with occasional layers of clay, peat and gravel of which the spatial extension is only marginally known (Wösten et al., 1985).A Miocene clay layer (20-30 m thick, starting at 0.5 to 20 m below the soil surface) forms an impermeable boundary for the unconfined water flow.The surface of this clay layer is carved by Pleistocene glacier erosion.
The entire catchment is densely drained with 68 km of ditches and many tube drains.The main brook is canalized (Fig. 4).A natural or reference situation is impossible to identify, because this catchment has been under continuous antropogenic change (canalization, re-meandering, land use change) for the last hundred years.The land use during the last ten years is mainly agricultural (maize and grass), with isolated farms and a few patches of forest.
The Hupsel brook catchment has a semi-humid sea climate with an annual precipitation of 500 to 1100 mm and an annual estimated evaporation of 300 to 600 mm, leaving an estimated sum of runoff and recharge of 200 to 800 mm year −1 .

Measured data
For the period 1994 through 2001 hourly weather data are available from a measurement station within the catchment operated by the KNMI (Royal Dutch Meteorological Institute) (Fig. 4).For the same period, discharges of the Hupsel brook were measured at the catchment outlet by the local waterboard with a 15 min interval using a calibrated weir.Groundwater levels were also recorded every 20 min in a monitoring well located at the meteorological station.For the period May 2007 through October 2008 weekly groundwater levels at 31 locations at a tube drained field site of 0.9 ha, located next to the meteorological station, were manually collected (Fig. 5).Within the catchment more than a 100 drilling logs were available to estimate the depth of the impermeable clay layer and the transmissivity.A Digital Elevation Model (DEM) was developed from radar data with a 5 m resolution.An estimate of the surface water levels in ditches and tributary brooks was obtained from this detailed DEM.

Groundwater model
The catchment water balance model requires the distribution of u, which obviously depends on the phreatic surface and the topography within the catchment.Since the former is not well-known and certainly not available with a high temporal resolution we resorted to modeling the phreatic aquifer of the Hupsel Brook catchment.We used a spatially distributed groundwater model with a 5 m resolution to test two major assumptions in the Theory section: -The Normality of the distribution of the thickness of the unsaturated zone within the catchment.
-The validity of Eq. ( 12) to describe the relation between the standard deviation of the thickness of the unsaturated zone at any given time and the average thickness of the unsaturated zone at that time.
The goal of this groundwater model, therefore, is not to represent the Hupsel Brook discharges and groundwater heads as accurately as possible, but to capture the most important flow processes like the wetting and drying of ditches and streams, tube drain drainage, and the effect of spatially distributed evapotranspiration so that we can establish the relation between the standard deviation of the thickness of the unsaturated zone and the average thickness of the unsaturated zone for the catchment.We therefore refrained from a detailed calibration of the model, since this was not expected to significantly change the relationship sought.
The groundwater model Modflow (McDonald and Harbaugh, 1988) was used to calculate the Darcian groundwa-ter flow with daily time steps for the period of 1994 through 2001.The model consisted of one unconfined layer of 740 by 800 cells.Transmissivity values were corrected for incisions of the brook and ditches and for the groundwater head, only taking into account the thickness of the wet crosssection (an unconfined simulation).Surface water levels were fixed to their annual average, with no flow of water from surface water to the soil allowed.Potential evapotranspiration was determined using the Makkink relation (Makkink, 1957) with temperature and global radiation measurements of the Hupsel meteorological station.To determine the actual transpiration for each cell, a relation with u was adopted.For 0<u≤0.7 m, e act =e pot .For 0.7<u<1.5 m, e act =e pot •(1.5−u)/0.8.For u≥1.5 m, e act =0.The effect of the unsaturated zone is modeled with an effective storage expressing the water layer needed for one meter of groundwater level rise.The value depends on soil type and the average local u, and varies between 0.08 for wet clayey soils and 0.26 for dry sandy soils.Because the main goal of this groundwater model was to mimic and not to exactly reproduce the natural groundwater flow these value were indicative and were not experimentally based.

Calibration and validation of the storage and flux model
The model developed in the Theory section (Eqs.12, 40 and 41) was calibrated on hourly measured catchment discharges for the period of 1 January 1994 to 1 January 1996, hourly measured groundwater depths at the meteorological station for the same period and an estimated yearly 59% contribution of tube drains to the total catchment discharge (estimation originates from Van der Velde et al., 2009).We have chosen an hourly time step because the time to peak of the catchment discharge after rainfall typically is a few hours.We adopted the fitted parameter values for Eq. ( 12) that relate σ u to u(t) from the groundwater model results and added 5 cm to σ min and σ max to account for additional soil surface elevation variation within 5×5 m model cells (this is an intuitive value and has not been validated by measurements).
Table 1 shows which model parameters were kept constant during calibration at their estimated value and which parameters were calibrated.Validation of the model was performed on similar data for the period 1996 through 2001.Within this period we selected the periods February 1997 through February 2000 and April 2001 through December 2001 (32 570 h) for the validation, because the quality of the catchment discharge data was good for these periods.Note that, for calibration and validation purposes, we had groundwater levels available for only a single location during this period.We considered those observations suitable, since the monitoring well was in the middle of a tube drained pasture field, approximately 100 m away from the nearest ditch.Therefore, we were confident that the values of u observed there were within the 20 percent (U 20 ) and 80 percent quantile (U 80 ) of www.hydrol-earth-syst-sci.net/13/1867/2009/ Hydrol.Earth Syst.Sci., 13, 1867-1885, 2009 Catchment fraction of un-drained and wet ---0.008all u within the catchment at all times.Including measured groundwater heads in the calibration (even at a single point) reduces the problem of model equifinality.The parameter estimation code PEST (Doherty, 2002) was used to optimize the model parameters for the objective function: is the measured discharge and l surf (t) •A is the modeled discharge at time step t.Variable EQ tot represents the error between measured and modeled fluxes, and EMQ tot accounts for the error in the cumulative mass flux during the simulation period between measured and modeled fluxes.The variable, EMq dr , accounts for the deviation in tube drainage contribution to the total discharge from the estimated 59%, and EH assures that the optimal parameter set gives a solution for which the measured groundwater head lies within the 20 to 80 percentile of the modeled distribution of groundwater depths.The weighting factors, 1.0, 2000, 2.0 and 5.0, for the respective components of the objective function were determined iteratively by running several optimization runs.These values ensure that each of the errors, Eqs.(43, 44, 45 and 46), contributed in the same order of magnitude to the final objective function, Eq. ( 42

Groundwater modeling
Figures 6 and 7 show the results of the distributed groundwater model for simulating measured discharges and groundwater heads.High groundwater heads and low discharges are overestimated and high discharges are underestimated.A sensitivity analysis showed that the phreatic storage coefficient was the most sensitive parameter to improve high flow or low flow model results.However, because we used a single coefficient for both flow conditions no significant improvements could be made.This also underlines the importance of the unsaturated-saturated zone interaction as implemented in our model in the Theory section.
For each time step, u was calculated from the Modflow results at all discretisation nodes.From this, the extent of water-filled drains, ditches and soil surface could also be found, thus allowing us to establish the extent of the active drainage network with time.Figure 8 illustrates the analysis for a dry and a wet day.The discharge was peaked, reflecting the efficiency of the drainage system during wet periods.The total spatial extent of the active drainage network differed dramatically between the dry and the wet situation.The necessity of the variable contributing area concept for groundwater flow to surface water is evident.The average u is much smaller during the wet period, as expected.For both events, the Normal distribution provides a good fit of the spatial distribution of u except for the hump around u=0.For the dry period this hump is caused by a few very deep incisions of ditches.Because there are only few deep incisions in the catchment the central limit theorem is not valid to describe With N the number of measurements 0.042 m 3 s −1 0.053 m 3 s −1 Q meas (t)−Q meas 2 0.87 0.78 their effect on u.During the wet period this hump is caused because the groundwater model removes all water above the average soil surface elevation in a grid cell not taking into account the possibility of ponding.The points in Fig. 9 represent the relation between the standard deviation, σ u , and the average thickness of the unsaturated zone, u , for the fitted Normal distributions for every simulated time step.Figure 9 corroborates the relationship between u and σ u hypothesized in the Theory section.Only one part (the string of outliers for 0.9< u <1.5 m) did not match the general trend.Figure 9 shows that Eq. ( 12) fits the data generated with the groundwater model well.

Field site results
Figure 10 shows the depth of the groundwater levels relative to the local surface elevation observed in the 31 monitoring wells installed in the 0.9 ha field.This graph quantitatively confirms that the spatial variation is large during wet periods and small during dry periods.The measured groundwater levels are spatially interpolated to obtain a groundwater table for the entire field site.Figure 11 shows the relation between u(t) and σ u within this field, together with a fit of Eq. ( 12).Within the range of groundwater depths measured at the field site Eq.( 12) gave a good fit.

Calibration results
Table 1 gives the fifteen parameters of the model.During the automatic calibration the eight parameters in the last column were kept constant, while the remaining parameters were allowed to vary between the maximum and minimum values of Table 1.We found significant non-uniqueness of the optimal dataset, which originates from the large correlation between storage and fluxes.Therefore, storage and fluxes should be determined separately (storage should not be derived from fluxes or fluxes from storages) by independently determining the parameters of Eq. ( 12) (that relate u(t) to σ u ), the parameters describing the unsaturated zone storage Eq. ( 22), and the surface storage s surf (t) .Since measurements of u(x,y,t) are only available at the field site, we added prior information to the PEST optimization to ensure that the optimal solution has: -Values for θ s , α and n close to the ranges for Dutch sandy soils reported by Wösten et al. (2001).
We visualized the model by means of eight characteristic curves: one representing σ u as function of u (Eq.12), three curves representing saturated, unsaturated and surface storage as a function of u (Eqs.16,22,24), and four curves giving the fluxes as a function of u .Figure 12 presents all eight curves.Figure 12b shows the relations between storage and u .The solid line represents the total pore space in the unsaturated zone, denoted by s sat,max − s sat , with s sat,max the total soil pore space.The difference between the curves for s sat,max − s sat and s unsat gives the catchment average air-filled pore space.Figure 12c shows the delicate balance between tube drain discharge q dr and discharge by streams ditches and overland flow q grw given by Eqs. ( 31) and ( 38): for u <0.9 m q grw is larger than q dr , for u >0.9 m q dr is larger than q grw .Figure 12d gives the fraction of precipitation that reaches and the fraction of potential evaporation that stems from the unsaturated zone.For small u , relatively large areas have surface storage (i.e.no unsaturated zone, see also s surf as a function of u in Fig. 12b).On locations with surface storage, precipitation is converted to discharge and evapotranspiration is subtracted from discharge appears in the last term of Eq. ( 41).
For large u evapotranspiration is reduced Eq.( 40).Both effects create the shape of the curves of Fig. 12d.Most of the discharge peaks were slightly underestimated, except for the discharge peak just after the summer dry period of 1994 which was simulated too high (Fig. 13).This resulted in an underestimation of the mean discharge by 3% (Table 2).Overall the hourly discharge was reproduced well (R 2 =0.88;Nash-Sutcliff (NS) coefficient=0.87,Nash and Sutcliffe, 1970).In contrast the root mean squared error (RMSE) was high compared to the average discharge.However the RMSE was dominated by errors during peak flow events between 0.1 and 1.5 m 3 s −1 , i.e. up to an order of magnitude larger than the average flux.
The model performed not so well for discharge events during dry conditions (the smaller discharge events around July, 1994 and July, 1995 in Fig. 13).We attribute this to the fact that under dry conditions only a small portion of the catchment generates discharge.Consequently, the number of fields involved in the discharge-generating process is too limited for the central limit theorem to apply.The assumption of a Normal distribution of u therefore becomes untenable.Figure 8d shows the distribution of u during a dry period.Overall, the normality of the PDF is convincing, but the generation of discharge in this situation is dominated by the few fields close to the sparsely distributed active drainage channels (including tube drains, Fig. 8b).These locations are represented by the small hump for u≈0 of the distribution of u.This hump is not described by the overall Normal distribution.
The groundwater levels measured at a single point at the field site were assumed to be within the 20% and the 80% quantile envelope of the spatial distribution of u, during calibration.This is visualized by Fig. 14.The measured data points lay within the dark gray area (the 20% to 80% quantile), but it is clear that the measured groundwater depths  were smaller than the modeled average groundwater depth, u .The measured location should therefore be relatively wet.The fact that the measurement field is tube-drained is consistent with this.
The total contribution of tube drains was somewhat lower than the 59% estimated by Van der Velde et al. (2009).This estimation, however, was based on the winter 2007-2008.Rainfall differences between years are likely to cause differences in the tube drain contribution.The sharp drops in tube drain contribution to total discharge in Fig. 13c during low discharge periods indicate a shift from tube drain discharge dominated to groundwater discharge dominated surface water.Only with surface water concentration measurements and a clear contrast between concentrations of tube drain flux and groundwater flux it is possible to calibrate A nd,wet A and to align these shifts with measured shifts in surface water concentrations.

Validation results
Table 2 shows that the average measured and calculated discharge for the validation period were close to those of the calibration period.The RMSE, however, increases to 53 L s −1 but also the extreme discharges during the validation period are much higher than during the calibration period.The R 2 and the NS coefficients of the validation decrease slightly to 0.85 and 0.78, respectively.The model performed well for the validation period, even for the high flows that were a factor two higher than the high flows of the calibration period (Fig. 15).
The model regularly overestimated discharge during autumn after a dry summer period (October, November and December 1994, 1997, 1999and 2001) and underestimated discharge during spring after a wet winter (March through June 1994, 1995, 1997, 1999and 2001).Possible sources of these errors are: -A slightly different relation between u(t) and σ u (t) when the groundwater table evolves from relatively parallel to the soil surface (low σ u (t) during summer) to a groundwater table with many large curvatures (high σ u (t)) between draining elements during autumn and winter than vice versa (from winter and spring to summer).
-The equilibrium assumption for the unsaturated zone storage overestimates unsaturated storage during evapotranspiration periods and underestimates unsaturated zone storage during infiltration periods.This reduces the precipitation amplifying nature of the unsaturated zone.
-Vegetation growth inside ditches and streams during summer and early autumn increases surface storage.After ditch cleaning in late autumn water is discharged more effectively with consequently higher peak discharges.-Systematic measurement and up-scaling errors in precipitation and evapotranspiration also contribute to the calculated errors in discharge.
Figure 16 shows distinct underestimations of the low flows as was already observed during calibration.Another difficulty with low flows is that they are far less accurate to measure because of the large dimensions of the weir and the abundant vegetation growth in and around the weir during summer.Particularly the latter leads to measurement errors that overestimate the true discharge, which would exaggerate the deviation from the 1:1 line in Fig. 16.
Infiltration excess overland flow is not incorporated in the model.Therefore, high discharge events due to high rainfall intensities, which occur mainly in summer, cannot be simulated accurately with the current model.This is shown in Fig. 16, where six discharge events that were measured were not simulated (the horizontal strands of data points under the 1:1 line).The modeled values of u deviated from the single-location values of u at the field site during the validation period (Fig. 17).Still, the deviations were nearly all contained within the envelope defined by U 20 and U 80 .

Conclusions
In lowland catchments without significant hillslopes, the depth to groundwater (thickness of the unsaturated zone) governs the various storage and flux terms in the water balance.We developed a model in which catchment-scale terms of the water balance are all expressed in terms of the PDF of the unsaturated zone thickness.By assuming this PDF to be Normal, a considerable reduction in the model complexity could be achieved.We demonstrated the ability of this parsimonious and uncomplicated model in a full calibrationvalidation cycle.While the potential of this novel approach to catchment modeling has been demonstrated, it is still unclear over which range of spatial scales the distribution of the unsaturated zone thickness remains Normal; small areas in particularly are likely to have deviating distributions.Furthermore, the relation between the shape of the distribution and properties of the drainage network and the transmissivity of the phreatic aquifer are still poorly understood.
This model offers great opportunities to improve our understanding of the interactions between groundwater and surface water.The model integrates subsurface and surface processes giving catchment-scale information about the interaction between groundwater and surface water, between groundwater and evapotranspiration, and the importance of the unsaturated zone and surface ponding during highdischarge events.
The model relies heavily on the relation between the average thickness of the unsaturated zone and its standard deviation.This relation cannot be derived from discharge data only, and needs to be derived from other data sources such as groundwater head measurements or a spatially distributed groundwater model to prevent large equifinality problems common to models of groundwater-surface water interactions (Beven, 2001).At the moment, measuring the average and standard deviation of the unsaturated zone thickness appears to be quite feasible in many catchments, possibly helped by remote sensing to quantify wet and dry fractions of a catchment.A well chosen nested-scales setup of discharge and groundwater head measurements could reduce the number of measurement needed.Monitoring programs in catchments aimed at determining the interactions between the groundwater and the surface water (which is relevant if the quality of the discharged water is of interest) should therefore be focused on quantifying the distribution of groundwater depths throughout the catchment in time.Reggiani et al. (1998) published a unifying framework for watershed thermodynamics, with conservation equations for mass, momentum, energy, and entropy, but leaving hydrologists struggling with the search for appropriate expressions for the interfaces between the different reservoirs (saturated zone, unsaturated zone, surface water) for their specific problems.In this paper we developed expressions for the interfaces between saturated zone, unsaturated zone, and surface water for typical lowland catchments based on the Normality assumption of the thickness of the unsaturated zone.
So far, the model has only been applied to one catchment and applications to new catchments will have to reveal the general applicability of the model concepts.The model results can be further improved by adding measurements of concentrations of selected compounds in various locations in the surface water, the soil, and the groundwater.Equally helpful are measurements that help quantify the contributions of individual flow routes, possibly leading to a water balance model that can accurately estimate the average travel time within the various reservoirs comprising the catchment.

Fig. 1 .
Fig. 1.Vertical cross-section of the Hupsel brook catchment in The Netherlands (see the main text for details).The surface elevation and the elevation of the impermeable thick clay layer are indicated, as well as the water levels of the brook that drains the catchment and of the ditches that discharge into the streams.Many of the fields in the catchment have tube drains, which are also indicated.Calculated groundwater levels on a wet (5 February 2001) and a dry day (8 July 1994) are also given.
e act [LT −1 ].The net lateral outward flux density through the subsurface is denoted by l sat [LT −1 ], and the net lateral outward flux density over the soil surface by overland flow, stream flow, and tube drain discharge by l surf [LT −1 ].No lateral fluxes in the unsaturated zone are assumed.Any sources and sinks are reflected by o[LT −1 ].The water balance of each of the storage compartments of Eq. ( represents the boundary of A at soil surface, l sat is the vertically integrated lateral flux density vector of the saturated zone [L 2 T −1 ], l surf is the vertically integrated lateral flux density of surface storage [L 2 T −1 ] and n is the outward normal vector of unit length [−] of S. The integrations convert flux densities [LT −1 ] to fluxes [L 3 T −1 ]

Fig. 2 .
Fig. 2. The water balance model describes fluxes at the point-scale.This figure illustrates three locations (x 1 ,y 1 ), (x 2 ,y 2 ) and (x 3 ,y 3 ) within a cross section of a typical lowland field.The groundwater level at location (x 1 ,y 1 ) is above soil surface, which leads to ponding.Note that when the groundwater level is above soil surface there is no unsaturated zone.Infiltrating water from the pond into the saturated zone is denoted q inf .Exfiltrating water from the saturated groundwater into the pond is denoted q ex .A sink is denoted o and the lateral overland flow l surf .Location (x 2 ,y 2 ) has an unsaturated zone and consequently no surface storage.The flux from the unsaturated zone to the saturated zone is denoted j rch and the capillary flux from saturated to unsaturated zone j cap .This location is also tube-drained with a tube drain flux q dr .Note that surface storage and tube drainage can occur at the same location.Point (x 3 ,y 3 ) is located at a stream.Above the stream bed surface storage occurs.The exfiltrating, infiltrating, and lateral surface fluxes are treated the same way for a ponded location (x 1 ,y 1 ) and a stream/ditch location(x 3 ,y 3 ).Rainfall, p, evapotranspiration, e act , and lateral saturated groundwater fluxes, l sat , occur in all three locations.

Fig. 3 .
Fig. 3. Pictures of the Hupsel brook catchment.(a) and (b), and (c) and (d) show the typical change in surface water level during a dry and a wet period resulting in changes in unsaturated zone thickness variation.(e) shows the large scale ponding that occurs during wet periods, and the resulting overland flow is shown in (f).

Fig. 6 .
Fig. 6.Observed daily groundwater levels at the weather station against Modflow calculations for the same location (2922 days).
). Evaluation of the objective function starts at time T start [T] (40 days), allowing for uncertainty in the starting value of u(0) , and runs until the time, T end [T].

Fig. 8 .
Fig. 8. (a) Measured daily discharge, with a wet and a dry day highlighted.The active (water draining) portion of the drainage network for the indicated dates is shown in (b) and (c).Panels (d) and (e) give the corresponding simulated distributions over the catchment of the thickness of the unsaturated zone (dots), and the fitted normal distribution (solid line).

Fig. 9 .
Fig. 9. Relation between daily values of the average thickness of the unsaturated zone, u(t) , and the standard deviation σ u of the spatial distribution of u, derived from the results of the groundwater model.The line represents a fit of Eq. (12).Arrow A shows the decline in variation when the catchment becomes dryer (larger value of u(t) ).Arrow B shows the decline in variation when the catchment becomes wet.

Fig. 10 .
Fig.10.Depth of groundwater levels relative to the local surface elevation (thickness of the unsaturated zone) for 31 wells at the field site within the Hupsel catchment (Fig.4).The dashed line shows the field average thickness of the unsaturated zone.The gray area represents the ranges between the 0.05 and 0.95 quantiles of a normal distribution around the mean.

Fig. 11 .
Fig. 11.Relation between measured average and standard deviation of thickness of the unsaturated zone.The line is a fit of Eq. (12).

Fig. 12 .
Fig. 12. Characteristic curves for the catchment scale: variation of the unsaturated zone thickness (a); unsaturated zone pore space s sat,max − s sat , unsaturated zone storage and surface storage (b); stream, ditch and overland flow discharge( q grw ), and tube drain discharge (c); and the fraction of precipitation that reaches and evapotranspiration that stems from the unsaturated zone (d).Reduction of precipitation that reaches the unsaturated zone occurs because part of the rain falls on locations with surface storage.Reduction of evapotranspiration is partly caused by surface storage (small u(t) ) and partly by and deep groundwater tables (large u(t) )).

Fig. 13 .
Fig. 13.(a) shows modeled discharge for calibration period.(b) shows the daily average error between measured and modeled discharge (measured-modeled) and (c) shows the modeled discharge subdivided into the contribution of tube drains, the contribution of stream, ditches, and land surface, and the contribution of direct rainfall in increasingly dark tones.

Fig. 14 .
Fig.14.The modeled spatial distribution of unsaturated zone thicknesses.The black dots are the measured groundwater depths at the meteorological station.Whenever the light gray area is below zero, more than 1% of the catchment soil surface contributes actively to discharge.Whenever the dark gray area is below 0.8 m, more then 20% of the catchment area has a groundwater level above the tube drainage level.Of this area a fraction of 0.6 (i.e.g * ) is tube drained and generates tube drain discharge.

Fig. 15 .
Fig. 15.(a) shows modeled and measured discharge for validation period.Only for the periods February 1997-February 2000 and April 2001-December 2001 we had good quality discharge data.(b) shows the daily average model error (measured-modeled discharge).(c) shows the contribution of individual flow routes to the total discharge.

Fig. 16 .
Fig. 16.Modeled versus measured hourly discharges for the validation period.The horizontally oriented strands of data points are observed high-discharge events which were not modeled.

Fig. 17 .
Fig. 17.Modeled versus measured hourly groundwater heads for the validation period in the monitoring well at the meteorological station.

Table 1 .
Calibration ranges and calibrated values for the model parameters (symbols explained in the main text).

Table 2 .
(Nash and Sutcliffe, 1970) results of the catchment model.RMSE is the Root Mean Squared Error.NS is the Nash-Sutcliff coefficient(Nash and Sutcliffe, 1970).R 2 is the squared Pearson correlation coefficient.