Hysteresis and scale in catchment storage, flow and transport

Hysteresis and scale in catchment storage, flow and transport

Jessica A. C. Davies ¹ * and Keith Beven ^1,2

1

Lancaster Environment Centre, Lancaster University, UK

2

Department of Earth Sciences, Uppsala University, Sweden

Abstract:

The closure problem of representing hydrological boundary fluxes given the state of the system has been described as the scienti fic ‘Holy Grail’ of hydrology. This relationship between storage state and flux should be hysteretic and scale dependent because of the differences between velocities and celerities in a hydrological system —effectively velocities are storage controlled, and celerities are controlled by storage de ficits. To improve our understanding of the nature of these relationships a new hydrology model is used (the Multiple Interacting Pathways or MIPs model) to explore the in fluence of catchment scale on storage –flow–transport relationships, and their non-linearities. The MIPs model has been shown to produce acceptable simulations of both flow and tracer, i.e. of both celerities and velocities, at the Gårdsjön catchment in Sweden. In this study the model is used to simulate scaled versions of the Gårdsjön catchment to allow us for the first time to investigate the influence of scale on the non-linearities in storage –flow–transport relationships, and help us steer the quest for the scientific hydrological

‘Holy Grail’. The simulations reveal the influence of scale on flow response in the nature of storage–discharge hysteresis and its links with antecedent storage; fractal-like systematic change of mean output travel times with scale; the effect of scale on input, output and storage residence time distributions; hysteric relations between storage and output travel times and links between storage and water table level hysteresis. © 2015 The Authors. Hydrological Processes published by John Wiley & Sons Ltd.

KEY WORDS catchment scale; residence times; storage –flux hysteresis; MIPs model; Representative Elementary Watershed Received 27 February 2014; Accepted 7 April 2015

INTRODUCTION

Arguably, the greatest challenge of hydrological model- ling remains getting the ‘right results for the right reasons’

(Kirchner, 2006). For too long, hydrological models that attempt to reproduce the processes of catchment hydro- logical response have been based on small-scale theory simply assumed to apply at larger scales —an approach that has been criticised at least since Beven (1989a).

Alternatives exist, such as the Representative Elemen- tary Watershed (REW) approach of Reggiani et al. (1999, 2000; Reggiani and Schellekens, 2003; Reggiani and Rientjes, 2005), which offers a scale-independent frame- work. The REW framework discretises the catchment into spatial elements termed REWs (which may be as large as the whole catchment or as small as the smallest de finable area), for which a set of generalised conservation equations of mass, momentum, energy and entropy are de fined. However, the REW approach raises the ‘closure problem ’ of formulating the boundary fluxes between calculation elements given the state of the system (Reggiani and Schellekens, 2003). To date, applications

of the approach have not addressed the ‘closure problem’, but rather simply assume that small-scale theory applies at the REW scale.

Consideration of the physics of the processes involved suggests that any such process representation should be hysteretic and scale dependent as a result of the differences between the celerities and velocities. This still does not seem to be fully understood by many hydrologists. Celerities or wave speeds represent the speed with which a perturbation to the flow propagates through the catchment to the outlet. Where that perturbation is a rainfall event then the distribution of celerities will control the shape of the hydrograph. The distribution of pore water velocities, on the other hand, will directly control the residence time and travel time characteristics of a conservative tracer for that event.

While celerities and velocities are clearly related it is important to understand that they are controlled by different mechanisms. Velocities are controlled by the pore size distribution and structure of the water filled pore space. Celerities are controlled by the local changes in storage associated with the effects of the perturbation.

Where the water table is responding, then the control is the storage de ficit that needs to be filled for the water table to rise, or the de ficit that is left when it falls.

Because this change of storage is often much less than the

*Correspondence to: J. A. C. Davies, Lancaster Environment Centre, Lancaster University, UK.

E-mail: j.davies4@lancaster.ac.uk

Published online 18 June 2015 in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1002/hyp.10511

© 2015 The Authors. Hydrological Processes published by John Wiley & Sons Ltd.

porosity of the soil then celerities are often much faster than the mean pore water velocity of the flow. For a more detailed discussion, including the simplest kinematic analysis, see Beven (1989b, 2006, 2012a Ch. 9) or McDonnell and Beven (2014).

Hysteresis is where an output depends not only on the current state of a system, but also on the input history and storage trajectory in time. The boundary fluxes required to close the mass/transport/energy balances will be reliant not only on the current storage state, but the history of wetting and drying, and as such can be described as hysteretic. The degree of hysteresis (i.e. the degree to which the storage varies between increasing and decreasing flux trajectories) should be scale dependent as the lags between pressure-wave celerities and transport velocities become more apparent as the scale increases. This scale dependence is important in controlling the difference between hydrograph and transport responses.

Nevertheless, there are very few modelling approaches that take any explicit account of the differences between celerities and velocities. Indeed, there are very few modelling approaches that have been tested against both flow and tracer responses, and those that have do not always show the right types of responses for both without modi fications (e.g. Christophersen et al., 1985; Lindström and Rodhe, 1986; de Grosbois et al., 1988; at Birkenes and Page et al., 2007 at Plynlimon). Thus, one of the ways of addressing the challenge of hydrological modelling at the catchment scale is to test models as hypotheses for reproducing both flow and tracer obser- vations (McDonnell and Beven, 2014).

The Multiple Interacting Pathways (MIPs) model is an exception as it represents both water flow and transport in one coherent solution and has been shown to reproduce both flow and transport for the small Gårdsjön catchment in Sweden within a hypothesis testing framework, using data collected by Allan Rodhe and others (Davies et al., 2011, 2013). We should, of course, expect a rather complex scale dependence of both runoff and residence times in hillslopes and catchments because of the nonlinear nature of the hydrological responses over complex wetting and drying cycles. Nonetheless, in the past, it has been suggested that catchments tend to more linear responses at larger scales in humid catchments, though less so in semi-arid or arid catchments (e.g.

Pilgrim, 1976; Dooge, 1986; Wood et al., 1990; Beven, 1991; Goodrich et al., 1997, Sivapalan et al., 2002).

Many catchment models take advantage of the expecta- tion that the routing of any effective volume of hillslope runoff in a storm will be more linear than the relationship between the total inputs to a catchment and the runoff (as in the models that incorporate unit hydrograph or linear transfer function concepts).

Here, we concentrate on the non-linearities and hysteresis inherent in the response of headwater areas, as if they were REWs in the sense of Reggiani and others. In this way we hope to suggest what might be needed in moving towards better closure schemes for the REW approach.

To do so, we will take a ‘virtual experiment’ approach (of the kind taken by Weiler & McDonnell 2004) by using the MIPs model to investigate the in fluence of catchment scale on the hysteretic relationships between water storage transport and flow. These relationships would be dif ficult to study based on observational data alone. It would be necessary to have good closure of the water balance for nested catchments at different scales, but which then would unavoidably each have their own unique characteristics of topography, soils, geology and land use (Beven, 2000), and uncertainties in each of the water balance terms (Beven, 2010, 2012b). However, creating and studying catchments within models negate these problems, provided that the model is scale independent and maintains mass balance —which MIPs do.

In what follows, the MIPs modelling concept is described, and its use in the creation of a suite of simulations of physically scaled catchments based on the Gårdsjön G1 catchment in Sweden is detailed. Results from these scaled simulations are then discussed, and the nature of scale-related hysteresis is explored.

The MIPs model

The Multiple Interacting Pathways (MIPs) model simulates both water flow and transport in one coherent solution, with direct representation of a continuum of preferential flows. This is an important departure from current methods of hydrological flow and transport representation where flow and transport are treated separately (usually through bolt-on transport solutions to continuum differential equation or storage accounting models) leading to a loss of coherence in representation of storage, flows and residence time. The MIPs model is also an attempt to overcome the rudimentary way in which preferential flows have been treated within dual-domain solutions (see Davies et al., 2013 for a full discussion).

The MIPs model achieves simultaneous heterogeneous flow and transport simulation through the combination of three concepts: random particle tracking, pathway velocity distributions and pathway exchange probabilities.

Water within the simulation is discretised into packages

or ‘particles’ whose movement is defined mechanistically

and whose position can be tracked throughout their

lifetime within the flow domain. Other particle properties

can be tracked alongside particle position, such as time of

entry, water age and source, time of exit and chemical

attributes. In this way, saturation and storage are

determined by the volume of particles within any de fined

area of the simulation domain, flow is determined by the number of particles flowing through a surface within the domain and transport is simulated through the analysis of particle attributes. For example, output travel times can be derived by examining the input time attributes of each particle on exit, and water origin can be determined by attributing a source label to each particle.

A continuous spectrum of preferential flow pathways can be simulated by the MIPs model using velocity distributions. At a given depth in the simulated soil pro file, the local mean velocity is approximately Darcian in that it is determined by the hydraulic conductivity, the water filled porosity and a constant local hydraulic gradient. However, the particle velocity is sampled from a continuous distribution of velocities around this mean, resulting in a solution where the flux per unit area integrates to the Darcian velocity, whilst representing transport at a range of velocities analogous to preferential flow features. Each particle then represents a volume of water flowing at that velocity in that part of the flow domain regardless of which pore spaces it is contained in.

In doing so, capillarity is neglected as only being important to slowly flowing matrix water, and gravity is assumed to drive the flow both vertically in the profile and laterally downslope. In de fining the velocity distribution, we are attempting to de fine the range and frequency of velocities permitted by preferential flow structures in the soil. There is no restriction on the velocity distribution other than it conforms to the mean velocity determined by hydraulic conductivity, porosity and gradient. In particular non-symmetrical or multi- modal distributions are possible including the simple (one parameter) exponential distributions used by Davies et al.

(2011, 2013) in previous applications to Gårdsjön.

Finally, pathway exchange probabilities determine how the pathways interact. Probabilities determining the likelihood of exchange between pathways can be de fined symmetrically (all pathways are equally likely to interact) or non-symmetrically (this can be a way of introducing capillarity effects, see Ewen, 1996) and can be based purely on frequency or on other system states such as soil moisture. If pathway exchange probabilities are set to zero, the velocities may be thought of as a Lagrangian representation of the particle velocity over its lifetime in the domain, but this would then be dif ficult to relate to any local velocity characteristics.

Applying MIPs at different spatial scales

Within this paper, 3D MIPs representations of scaled hypothetical versions of the Gårdsjön G1 catchment are used as ‘virtual experiments’ (in the parlance of Weiler and McDonnell, 2004 who call for such exploratory work in the improvement of hydrological understanding) to

explore the nature of scale-dependent hysteresis on storage –flux and storage–residence-time relationships in a controlled and directly comparable manner.

The MIPs concept allows such an exercise, as the concept is readily scalable. The particle tracking meth- odology used within MIPs is largely grid-scale indepen- dent: the grid is only used after moving the particles in a grid-free manner at the end of each time step to determine the local water table levels, and because the movement is approximated by a kinematic flux relationship the local water table position does not in fluence the local hydraulic gradient.

The G1 catchment is a small (6300 m ² ) catchment, forested with Norway Spruce and underlain by gneissic rock, in the Gårdsjön area near the west coast of Sweden.

The catchment has a long central valley bordered by steep slopes. The catchment has predominantly podzolic soils, which are generally thinner on the steep slopes and deeper in the valley regions. The G1 catchment has been monitored since 1988 and has been the site of numerous manipulation experiments —many involving Allan Rodhe (Rodhe, 1985; Nyberg et al., 1993, 1999; Rodhe et al., 1996; Bishop et al., 2011; Seibert et al., 2011) — including the ‘roof’ experiment where a plastic roof was installed below the canopy to intercept natural rainfall, below which the catchment was sprinkled with water pumped from the nearby Gårdsjön lake. The wealth of driving data (sprinkled fluxes) and the small nature of the catchment makes it a good candidate for this scaling experiment albeit that the available input data have a limited (daily) time resolution.

The G1 catchment is virtually scaled downwards and upwards to provide three simulation catchments of varying scale, as shown in Figure 1. Initial representation of the G1 catchment (G1M), in terms of both topography and hydrological representation, is made on a 5 × 5 m grid. The catchment is scaled downwards to an area of 1575 m ² (a quarter of the original size) using the same topography, but changing the grid size to 2.5 × 2.5 m (denoted as G1S). Likewise the catchment is scaled upwards to an area of 25 200 m ² (four times the original size) by doubling the initial grid size to 10 × 10 m (denoted as G1L). For this initial investigation, the elevations at each grid point are not changed, such that the larger and smaller catchments also have shallower and steeper slopes. Keeping the number of scales examined small also keeps the computational requirements at a manageable level, as the MIPs model is currently computationally expensive. Each scale is run with the same sequence of input rates.

The 3D MIPs model for the initial catchment is

de fined as described in Davies et al. (2013), where the

G1 catchment-scale isotope experiment was simulated,

and as such the description of the methodology is not

repeated here. The parameters of this model are based on physically measured properties, and as such were not calibrated.

The only changes required to make the scaled versions of the catchment are in the de finition of the grid size (as explained above) and to the discretisation parameters. The particle volume is scaled from 2 l in the initial catchment to 0.5 and 8 l for the smaller and larger scale catchments, G1S and G1M, respectively. This is done in order to keep the number of particles entering each grid cell per mm of rainfall the same, to ensure that effects seen in the results are not attributable to discretisation scale change and only attributable to catchment scale change.

It is important to note that in the simulation experiments that follow the particles move continuously through space with the velocity and direction assigned at the current time step, i.e. the grid size enters into the calculation only in the determination of a local water table according to the total volume of particles moving laterally within that grid square. Mass balance is maintained by de finition because all particles are accounted for. Because of the assumption of gravity drainage according to the local slope, there is no need to calculate a hydraulic gradient between grid squares —the movement of each particle is kinematic in that respect (Davies and Beven, 2012). This minimises any numerical effects of a changing grid size.

Scaled catchment results

The three scaled models were initialised using a spin- up period of 3 years in order to wet up the catchment.

Driving data from the G1 experimental catchment monitoring were used in order to provide typical patterns of wetting and drying. Following the spin-up period, the catchments were driven with data from the G1 catchment spanning from February 1990 to December 1993.

In the results that follow, discharges, travel and residence times and storages derived from these virtual experiments are used to explore five assertions regarding scale-dependent catchment response, namely that:

i. Increasing scale damps the flow response,

ii. Increasing scale intensi fies storage discharge hystere- sis,

iii. Mean output travel time scales with catchment area systematically,

iv. The nature of input, output and storage travel and residence time distributions is in fluenced by catchment scale,

v. And water table levels exhibit hysteretic behaviour which is in fluenced by catchment scale.

These are discussed in turn in the subsections that follow.

i. Increasing scale damps flow response

Figure 2 gives the hydrographs for the scaled simula- tions between January 1991 and 1993. Discharge is derived in the MIPs model by counting the number of particles exiting the catchment over a time step. It can be seen from these hydrographs that changes in scale alone in fluence the hydrological flow response, producing a more damped response as catchment scale increases. Delay to peak flow following a rainfall event is increased with scale, because of greater travel distances and lower slopes such that the magnitude of the peak flow per unit area is reduced. Likewise, increases in travel distances and lower slopes for propagation of saturated conditions result in longer recessions after peak flow. There is also a minor effect on the total speci fic discharge, which decreases with catchment scale because of the longer travel times of the water allowing further evapotranspiration. This effect will

Figure 1. Scaled versions of the G1 catchment at Gårdsjön used within the scaling simulations

be enhanced by the assumption in this study that the grid elevations are kept constant such that average slope decreases with scale.

ii. Increasing scale intensi fies storage–discharge hysteresis It is well known that storage regulates flux in a non- linear fashion creating hysteresis in the storage–discharge relationship (Beven, 2006, Zehe et al., 2007; Detty and

McGuire, 2010). The MIPs model produces hysteretic behaviour as demonstrated in Figure 3 which gives the storage –discharge relations for each catchment. To normalise the plots between scales, fractional discharge is plotted against fractional storage, where fractional discharge Q(t) is the discharge volume Q v (t) expressed as a fraction of the potential storage volume S vtotal of the catchment:

Figure 2. Hydrographs for the scaled G1 catchments between 1991 and 1993

Figure 3. Fractional storage –discharge relationships for the three scaled catchments between January 1991 and January 1993. Fractional storage is the current storage as a fraction of the total storage. Fractional discharge is the current discharge as a fraction of the total storage. The plots show that the relationship between catchment storage and discharge is hysteretic, i.e. that the discharge at any point in time not only depends on the current storage (if

this was the case these relations would fall on one trajectory), but also on the history of wetting and drying

Q t ð Þ ¼ Q v ð Þ=S t vtotal (1) And similarly the fractional storage is the storage volume at time t S v (t) expressed as a fraction of the potential storage of the catchment:

S t ð Þ ¼ S v ð Þ=S t vtotal (2) Hysteresis is evident in these plots as the fractional storage increases on one trajectory as fractional discharge increases following a rainfall event, and falls on a separate trajectory as the discharge subsequently de- creases. However, it is easier to compare hysteresis across the three simulations by looking at smaller time slices.

Hence, to look at part of this data more closely, Figure 4 gives hydrographs and storage –discharge relations for two sets of events in winter/spring 1991 and 1992. Events during winter and spring are focused upon as the in fluence of evapotranspiration during these periods are low, allowing clearer interpretation of storage discharge relationships. It can be seen most clearly in the storage – discharge plot for the 1992 event sequence (bottom right pane, Figure 4) that antecedent storage conditions control the degree of hysteresis to some extent as the decline in initial storage between similar events at the start of March 1992 and April 1992 results in greater hysteresis during the April event when compared to the March event.

There are many reasons for storage –discharge hysteresis occurring in catchments through the factors controlling celerities and velocities, e.g. the effects of antecedent conditions on the storage de ficit profile for the propagation

of event perturbations, the role of horizontal preferential flow pathways and vertical fingering during wetting in bypassing matrix storage, changes in contributing areas between wetting and drying phases, changing flow path connectivi- ties, and fill and spill mechanisms. In the case of the MIPs model here, hysteresis is introduced via the representation of preferential pathways —the fastest pathways drain first in the drying phase, leaving water travelling in slower pathways. It will also be a function of changes in contributing areas during wetting and drying, i.e. the catchment flanks are more active during wetting than in drying. Antecedent storage will in fluence both of these effects as well as the deficit profile controlling local celerities in the saturated zone.

Scale also acts as a control on storage –discharge hysteresis. It can be seen in both of the event periods highlighted by Figure 4 (although perhaps most clearly in the 1991 event period) that increasing catchment scale increases the degree of hysteresis (there is more diversion between wetting and drying limbs in the storage-discharge relationships). This is attributable to the longer time scales involved in wetting and drying phases as spatial scale is increased and, in this case, averages slopes decrease.

iii. Mean output travel time scales with catchment area systematically

The output travel time distribution is the age distribution of the water making up the output at any time step. With increasing catchment size and travel pathways, it follows that the mean output travel time τ out should increase with catchment area.

Figure 4. Hydrographs and storage –discharge relationships (as in Figure 3) across scales for two event periods in winter–spring 1991 (top panels) and

winter –spring 1992 (bottom panels). Arrows are given in the storage–discharge plots to denote direction

Output travel time can be dynamically calculated using the MIPs model by examining the input times of particles exiting the catchment within any time period. Figure 5 gives a time series of mean output travel time calculated over a daily time period for the scaled G1 catchments.

Figure 5 shows that in the MIPs representation mean output travel time is increased with scale, and that the mean output travel time is non-stationary because of changes in storage and flow conditions.

The inset in Figure 5 gives the overall mean output travel time over 1991 and 1992 for each catchment. It can be seen that this overall mean scales approximately as area to the power of the third, i.e. τ out / A ¹⁼³ suggesting that mean travel times scale in a systematic manner, displaying fractal-like behaviour.

iv. The nature of input, output and storage travel and residence time distributions are in fluenced by catchment scale

Output travel times refer to the time an exiting volume of water has spent within the flow domain. In contrast, the input travel time refers to the time an entering volume of water will spend travelling in the flow domain. Input travel time is derived in the MIPs model by examining the output time of particles that entered the catchment during a de fined input time period.

As highlighted in Rinaldo et al 2011, the input and output travel time distributions should not be expected to be the same. In addition, there is storage residence time, which is the age of the water in the flow domain at any point in time. This can be derived in MIPs by examining

the input times of water in the simulated storage volume at any point in time.

It is to be expected that flow domain scale will in fluence each of these travel and residence time distributions. Figure 6 gives the cumulative distributions of the input, output and storage travel and residence times, calculated for each day in 1991 and 1992 within the scaled MIPs simulations. As scale increases the cumulative distributions are pushed towards the bottom left, indicating that:

• input travel times increase with scale because of increased pathway lengths and lower slopes;

• discharge in smaller catchments is made up of relatively more new water and discharge in larger catchments contains relatively more older water;

• the age distribution of the storage volume skews towards older waters in larger catchments.

There is a large variation in the travel and residence time distributions, particularly the output travel time distribution. This variation is related to antecedent storage and flow in a non-linear manner. Figure 7 illustrates how mean output travel time varies with storage and discharge over the first time period that was examined in relation to storage –discharge hysteresis in Figure 4 (February to April 1991). This time series has been split into three event slices: a), b) and c), the travel time –discharge and travel time –storage relations for which are given for each event. These relationships between travel time and storage/discharge also display hysteretic behaviour. The

Figure 5. Daily mean output travel time for each catchment scale, with inset figure showing overall mean output residence time plotted against catchment

area. The A

line is plotted alongside the overall mean output residence time

Figure 6. Cumulative distributions of input (top row), output (middle row) and storage (bottom row) travel and residence times for G1S, G1M and G1L for each day between January 1991 and January 1993. Mean travel and residence time distributions are shown in grey solid lines, and the 10th and 90th

percentile in dashed grey

Figure 7. The top panes focus in on the events from 20 February to 20 April 1991, giving precipitation and discharge. The event series is split into three slices: a), b) and c). For three time slices the storage-to-mean-output-travel-time (S –τ

) (middle panes) and the (lower panes) discharge-to-mean-output- travel-time (Q –τ

) relations are given. Clear hysteresis and catchment scale effects are seen in the Q –τ

relations, with hysteresis increasing with scale.

The relation between S and τ

is less clear. Hysteresis is displayed during events a) and c); however, there is no fixed trend with scale. In event a) the

smallest catchment produces the most hysteresis, and in event c) the largest catchment produces the most hysteresis

travel time –discharge relationships distinctly show hys- teresis loops that consistently increase in magnitude with scale. The travel time –storage relations are less clear.

There is very little hysteresis between mean residence time and storage during event b), although the largest catchment displays some looping. Event c) displays more branching between wetting and drying limbs, with greater hysteresis at larger scale again. Within event a), however, the trend is reversed, with the G1S scale presenting the most hysteresis.

v. Water table levels exhibit hysteretic behaviour which is in fluenced by catchment scale

Hysteretic behaviour is often observed in water table levels: water tables are observed to rise rapidly during rainfall and recede much more slowly after cessation of rain (e.g. Gillham, 1984; McDonnell, 1990; Vidon, 2012). Indeed, this has been observed to occur at the G1 catchment (Seibert et al 2011).

This behaviour is also witnessed in the MIPs simulations. Figure 8 gives the water table levels (given in height as a fraction of total depth to bedrock) for locations corresponding to locations of groundwater monitoring at the G1 catchment (see Nyberg 1995).

Rapid rise and delayed recession in the water table are particularly noticeable in GV8, GV141 and GV21 where the response to rainfall is more ‘flashy’. Because of the damping of response with scale, increased scale results in

more widespread levels of increased saturation. This raising of the water table leads to more pronounced responses in the more moderately saturated slopes (e.g.

GV10 and GV132), in comparison to those in the smaller catchments.

In order to examine this behaviour across space and time more thoroughly the difference between the mean wetting and mean drying water table velocity is examined. The mean water table velocity is calculated for each simulated grid cell water table level during water table rises v rise and falling phases v fall . These are then summed (v rise þ v fall ) to provide a measure of the difference in rise rate to fall rate. A positive value in v rise þ v fall denotes that, on average, the grid cell ’s water table rises more quickly than it falls. The cumulative density function (cdf) of v rise þ v fall is plotted in Figure 9 alongside the spatial variation in v rise þ v fall . The cdf shows that all of the catchment simulations produce positive values in each grid cell, indicating that everywhere in the catchment water tables, on average, rise faster than they fall. However, the degree to which water tables rise faster than they recede varies across each catchment and between catchment scales. Because of the increased area of initial saturation at larger scales, the water table level response to rainfall is more rapid in the slopes, and recessions longer lived in the valley regions, accounting for the movement of the cdf curve to the left and a general convergence towards similar values of v rise þ v fall across the catchment. In the smallest

Figure 8. Fractional water table heights (height as a fraction of the depth to bedrock) for scaled simulation cells corresponding to groundwater tubes as in

Nyberg (1995)

catchment the quick draining nature of the slopes means that water tables remain low, and the valley region is more flashy, leading to a greater spatial difference in v _rise þ v fall .

DISCUSSION

Hysteresis concepts are not new to hillslope hydrology but in the past have mostly been concerned with the small-scale hysteresis in the soil moisture characteristics within a Darcian modelling framework that includes capillarity effects (Mualem, 1974; Jaynes, 1984; Kool and Parker, 1987; Haverkamp et al., 2002; Huang et al., 2005). It is worth noting in passing that the bypassing effects associated with dynamic non-equilibrium fluxes in unsaturated soils can similarly result in local hysteresis in ways that cannot be directly represented by the Darcy – Richards equation. The original experiments of Richards were, indeed, designed to exclude the possibility of non- equilibrium flows, although he did report on differences in the flux relationships between wetting and drying (Beven and Germann, 2013; Beven, 2014).

The hysteresis we consider in this paper is at the small catchment scale, integrating over the heterogeneity in soil properties, slope properties, length scales, convergence and divergence. At larger scales, we may expect further hysteresis induced by further processes such as hillslope and stream flow connectivity (Smith et al. 2013). The hysteresis here at these small catchment scales is produced essentially by the difference between transport velocities and celerities on the slopes. It can be considered at the scale of a local calculation element, but the local responses are conditioned by the upscale responses which will also be a function of scale. It can also be considered at the full catchment scale as Representative Elementary Watersheds in the sense of Reggiani et al. (1999, 2000), and the results presented here demonstrate the importance of accounting for such hysteresis in REW closure schemes, as first suggested by Beven (2006).

That is, however, easier said than done, because the hysteresis demonstrated does not appear to be simply structured. There are theories of hysteretic relationships that might be invoked in representing storage –flux relationships, but these tend to be stationary and rate independent in nature (e.g. Krasnosel ’skii and Pokrovskii,1989; Mayergoyz, 1991; as adopted in hydrology by O ’Kane, 2005; O’Kane and Flynn, 2007; and Flynn et al. 2008).

The hysteresis exhibited in the MIPs simulations of the different scales of REW is, as might be expected, rate and antecedent condition dependent (as was also the case in the analysis of observed data presented by Beven, 2006).

At the scale of the small Gårdsjön catchment we have shown how this can be modelled directly, at least under the relatively wet gravity driven conditions of these MIPs simulations. We have done so in a way that provides information about the different (time variable) residence time distributions in the system as well as the discharge responses. That is a start but, at least for the present, would be prohibitively expensive computationally at larger catchment scales made up of many such REW elements. This suggests that a summary functional representation of the hysteretic storage –flux relationship is needed to act as a general closure scheme for larger scale models (see, for example Ewen and Birkenshaw, 2007). How to develop such a relationship from either modelled or field data remains the subject of future work.

CONCLUSIONS

This exercise has demonstrated the power of the MIPs model concept for exploring hysteretic relations at different catchment scales as it is able to simultaneously represent storage, flows and transport processes, includ- ing the effects of preferential flows. The results have demonstrated how the storage –discharge hysteresis is antecedent wetness, rate and scale dependent in rather complex ways at both the catchment scale and at individual local grid elements. It has also been shown

Figure 9. Cumulative distributions for the difference between mean water table rising velocity v

and falling velocity v

are shown alongside the

spatial pattern for v

þ v

for each catchment size

how the different forms of residence time distributions are scale dependent and non-stationary in time. The hypothetical simulations have been useful for exploring these issues, and testing whether the model outputs are consistent with our conceptual understanding. While there have been no conceptual surprises in the results it is clear that developing a functional representation of the form of hysteresis that could be used in parameterising REW closure schemes in larger catch- ments made up of many REW elements is likely to be challenging. The results suggest that further explorations over larger scales, and changing topographies and slopes would be useful.

ACKNOWLEDGEMENTS

This work was funded by the UK Natural Environment Research Council (NERC) under grant reference NE/G017123/1. Thanks are extended to Allan Rodhe and his colleagues who carried out the G1 experiments at Gårdsjön that inspired this work.

REFERENCES

Beven KJ. 1989a. Changing ideas in hydrology: the case of physically based models. Journal of Hydrology 105: 157 –172.

Beven KJ. 1989b. Inter flow. In Proc. NATO ARW on Unsaturated Flow in Hydrological Modelling, Morel-Seytoux HJ (ed.) Kluwer: Dordrecht;

191 –219.

Beven KJ. 1991. Scale considerations. In Proc. NATO ASI on Recent Advances in the Modelling of Hydrological Systems, Bowles DS, O ’Connell PE (eds). Springer: Netherlands; 357–371.

Beven KJ. 2000. Uniqueness of place and process representations in hydrological modelling. Hydrology and Earth System Sciences 4(2):

203 –213.

Beven KJ. 2006. The holy grail of scienti fic hydrology: Q

¼ H ←S←R ð ÞA as closure. Hydrology and Earth Systems Science 10: 609 –618.

Beven KJ. 2010. Preferential flows and travel time distributions: defining adequate hypothesis tests for hydrological process models. Hydrolog- ical Proceses 24: 1537 –1547.

Beven KJ. 2012a. Rainfall –Runoff Modelling—The Primer, 2

. edition.

Wiley-Blackwell: Chichester UK.

Beven KJ. 2012b. Causal models as multiple working hypotheses about environmental processes. Comptes Rendus Geoscience, Académie de Sciences, Paris 344: 77 –88. DOI:10.1016/j.crte.2012.01.005.

Beven KJ. 2014. Penman Lecture: “Here we have a system in which liquid water is moving; let ’s just get at the physics of it” (Penman 1965).

Hydrology Research 45(6): 727 –736. DOI:10.2166/nh.2014.130.

Beven KJ, Germann PF. 2013. Macropores and water flow in soils revisited. Water Resources Research 49(6): 3071 –3092. DOI:10.1002/

wrcr.20156.

Bishop K, Seibert J, Nyberg L, Rodhe A. 2011. Water storage in a till catchment. II: Implications of transmissivity feedback for flow paths and turnover times. Hydrological Processes 25(25): 3950 –3959.

Christophersen N, Kjærnsrød S, Rodhe A. 1985. Preliminary evaluation of flow patterns in the Birkenes catchment using natural

0 as a tracer, Hydrologic and Hydrogeochemical Mechanisms and Model Ap- proaches to the Acidi fication of Ecological Systems, Rep. 10, Nord.

Hydrol. Programme, Oslo, Norway; 29 –40.

Davies J, Beven KJ, Nyberg L, Rodhe A. 2011. A discrete particle representation of hillslope hydrology: hypothesis testing in reproducing a tracer experiment at Gårdsjön, Sweden. Hydrological Processes 25:

3602 –3612. DOI:10.1002/hyp.8085.

Davies J, Beven KJ. 2012. Comparison of a multiple interacting pathways model with a classical kinematic wave subsurface flow solution. Hydrolog- ical Sciences Journal 57: 203 –216. DOI:10.1080/02626667.2011.645476.

Davies J, Beven KJ, Rodhe A, Nyberg L, Bishop K. 2013. Integrated modelling of flow and residence times at the catchment scale with multiple interacting pathways. Water Resources Research. DOI:10.1002/wrcr.20377.

De Grosbois E, Hooper RP, Christophersen N. 1988. A multisignal automatic calibration methodology for hydrochemical models: a case study of the Birkenes model. Water Resources Research 24(8): 1299 –1307.

Detty JM, McGuire KJ. 2010. Topographic controls on groundwater dynamics:

implications for hydrologic connectivity between hillslopes and riparian zones in a till mantled catchment. Hydrological Processes 24: 2222 –2236.

Dooge JC. 1986. Looking for hydrologic laws. Water Resources Research 22(9S): 46S –58S.

Ewen J. 1996. ‘SAMP’ model for water and solute movement in unsaturated porous media involving thermodynamic subsystems and moving packets:

2. Design and application. Journal of Hydrology 182: 195 –207.

Ewen J, Birkenshaw SJ. 2007. Lumped hysteretic model for subsurface storm flow developed using downward approach. Hydrological Pro- cesses 21: 1496 –1505.

Flynn D, Zhezherun A, Pokrovskii A, O ’Kane JP. 2008. Modeling discontinuous flow through porous media using ODEs with Preisach operator. Physica B: Condensed Matter 403.2(2008): 440 –442.

DOI:10.1016/j.physb.2007.08.070

Gillham RW. 1984. The capillary fringe and its effect on water-table response. Journal of Hydrology 67(1 –4): 307–324, ISSN 0022-1694 DOI: 10.1016/0022-1694(84)90248-8

Goodrich DC, Lane LJ, Shillito RM, Miller SN, Syed KH, Woolhiser DA.

1997. Linearity of basin response as a function of scale in a semiarid watershed. Water Resources Research 33(12): 2951 –2965.

Haverkamp R, Reggiani P, Ross PJ, Parlange JY. 2002. Soil water hysteresis prediction model based on theory and geometric scaling. In Environmental Mechanics: Water, Mass and Energy Transfer in the Biosphere. The Philip Volume. Geophysical Monograph 129. The American Geophyscial Union. 213.

Huang HC, Tan Y-C, Liu CW, Chen CH. 2005. A novel hysteresis model in unsaturated soil. Hydrological Processes 19: 1653 –1665.

Jaynes DB. 1984. Comparison of soil-water hysteresis models. Journal of Hydrology 75: 287 –299.

Kirchner JW. 2006. Getting the right answers for the right reasons:

Linking measurements, analyses, and models to advance the science of hydrology. Water Resources Research 42: W03S04 DOI: 10.1029/

2005WR004362

Kool JB, Parker JC. 1987. Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties. Water Resources Research 23: 105 –114.

Krasnosel ’skii M, Pokrovskii A. 1989. Systems with Hysteresis. Springer- Verlag: New York.

Lindström G, Rodhe A. 1986. Modeling water exchange and transit times in till basins using oxygen-18. Nordic Hydrology 17: 325 –334.

Mayergoyz ID. 1991. Mathematical Models of Hysteresis. Springer-Verlag:

New York.

McDonnell JJ. 1990. A rationale for old water discharge through macropores in a steep, humid catchment. Water Resources Research 26(11): 2821 –2832. DOI:10.1029/WR026i011p02821.

McDonnell JJ, Beven KJ. 2014. Debates —the future of hydrological sciences:

a (common) path forward? A call to action aimed at understanding velocities, celerities, and residence time distributions of the headwater hydrograph.

Water Resources Research 50: . DOI:10.1002/2013WR015141.

Mualem Y. 1974. A conceptual model of hysteresis. Water Resources Research 10: 514 –520.

Nyberg L. 1995. Soil- and groundwater distribution, flowpaths and transit times in a small till catchment, Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 97, 35 pp.

Nyberg L, Bishop K, Rodhe A. 1993. Importance of hydrology in the reversal of acidi fication in till soils, Gårdsjön, Sweden. Applied Geochemistry 8: 61 –66.

Nyberg L, Rodhe A, Bishop K. 1999. Water transit times and flow paths from two line injections of 3H and 36Cl in a microcatchment at Gårdsjön, Sweden. Hydrological Processes 13(11): 1557 –1575.

O ’Kane JP. 2005. The FEST model—a test bed for hysteresis in hydrology

and soil physics. Journal of Physics: Conference Series 22: 148.

O ’Kane JP, Flynn D. 2007. Threshold switches, thresholds and hysteresis in hydrology from the pedon to the catchment scale: a nonlinear systems theory. Hydrology and Earth System Sciences 11: 443 –459.

Page T, Beven KJ, Freer J. 2007. Modelling the chloride signal at the Plynlimon catchments, Wales using a modi fied dynamic TOPMODEL.

Hydrological Processes 21: 292 –307.

Pilgrim DH. 1976. Travel times and nonlinearity of flood runoff from tracer measurements on a small watershed. Water Resources Research 12(3):

487 –496.

Reggiani P, Hassanizadeh SM, Sivapalan M, Gray WG. 1999. A unifying framework for watershed thermodynamics: constitutive relationships.

Advances in Water Resources 23: 15 –39.

Reggiani P, Rientjes THM. 2005. Flux parameterization in the representative elementary watershed approach: application to a natural basin. Water Resources Research 41: W04013. DOI:10.1029/

2004WR003693.

Reggiani P, Schellekens J. 2003. Modelling of hydrological responses: the representative elementary watershed approach as an alternative blueprint for watershed modelling. Hydrological Processes 17(18): 3785 –3789.

Reggiani P, Sivapalan M, Hassanizadeh SM. 2000. Conservation equations governing hillslope responses: exploring the physical basis of water balance. Water Resources Research 36(7): 1845 –1863.

Rinaldo A, Beven KJ, Bertuzzo E, Nicotina L, Davies J, Fiori A, Russo D, Botter G. 2011. Catchment travel time distributions and water flow in soils. Water Resources Research 47: W07537. DOI:10.1029/

2011WR010478.

Rodhe A. 1985. Groundwater contribution to stream flow in the Lake Gårdsjön area. Ecological Bulletins 37: 75 –85.

Rodhe A, Nyberg L, Bishop K. 1996. Transit times for water in a small till catchment from a step shift in the oxygen 18 content of the water input.

Water Resources Research 32(12): 3497 –3511.

Seibert J, Bishop K, Nyberg L, Rodhe A. 2011. Water storage in a till catchment. I: Distributed modelling and relationship to runoff.

Hydrological Processes 25(25): 3937 –3949.

Sivapalan M, Jothityangkoon C, Menabde M. 2002. Linearity and nonlinearity of basin response as a function of scale: discussion of alternative de finitions. Water Resources Research 38(2): 1012.

Smith T, Marshall L, McGlynn B, Jencso K. 2013. Using field data to inform and evaluate a new model of catchment hydrologic connectivity.

Water Resources Research 49: 6834 –6846. DOI:10.1002/wrcr.20546.

Vidon P. 2012. Towards a better understanding of riparian zone water table response to precipitation: surface water in filtration, hillslope contribution or pressure wave processes? Hydrological Processes 26:

3207 –3215. DOI:10.1002/hyp.8258.

Weiler M, McDonnell J. 2004. Virtual experiments: a new approach for improving process conceptualization in hillslope hydrology. Journal of Hydrology 285(1): 3 –18.

Wood EF, Sivapalan M, Beven KJ. 1990. Similarity and scale in catchment storm response. Reviews Geophysics 28: 1 –18.

Zehe E, Elsenbeer H, Lindenmaier F, Schulz K, Bloschl G. 2007. Patterns

of predictability in hydrological threshold systems. Water Resources

Research 43: W07434. DOI:10.1029/2006WR005589.