Velocities, celerities and the basin of attraction in catchment response

Velocities, celerities and the basin of attraction in catchment response

Keith Beven ^1,2,3 * and Jess Davies ¹

1

Lancaster Environment Centre, Lancaster University, Lancaster LA1 4YQ, UK

2

Department of Earth Sciences, Uppsala University, Uppsala, Sweden

3

IDYST, University of Lausanne, Lausanne, Switzerland

Abstract:

Catchment systems are interestingly nonlinear, but their dynamics are constrained from being unduly chaotic by mass and energy balance requirements. There have been no attempts in hydrology that we know of that have tried to map both the flow and transport dynamics of a catchment in any form of phase space. In part, this is because of the high dimensionality of the space – time patterns of response; in part because there is suf ficient uncertainty about the input and output fluxes estimated by measurement that this might be expected to obscure any attractor-like behaviour. In this study we explore the basin of the catchment attractor for the Multiple Interacting Pathway (MIPs) model that in previous papers has been shown to give good results for the small Gårdsjön catchment in Sweden. MIPs is based on particle tracking techniques and gives results for both the flow responses and for the travel and residence time responses of water in the catchment. Here it is used to provide consistent values of fluxes, total storage, travel time distributions and residence time distributions for a long simulation period. The nature of those responses in storage and input dimensions is then investigated. The results suggest that the range of behaviours is hysteretic in interesting ways and constrained by the forcing inputs, with space filling of trajectories in the basin of attraction as should be expected of a forced dissipative system. The range of behaviours exhibited de fines a space that the responses of any simpler emulator model will need to span. © 2015 The Authors Hydrological Processes published by John Wiley & Sons Ltd.

KEY WORDS nonlinearity; basin of attraction; travel times; residence times; hysteresis; storage –discharge Received 17 January 2015; Accepted 24 September 2015

INTRODUCTION: NONLINEAR DYNAMICS IN CATCHMENT SCIENCE

We know that catchment responses to rainfall and evapotranspiration forcing are nonlinear in their nature, making the problem of predicting storm discharges challenging. Discharges in any event are expected to vary with the spatial and temporal patterns of rainfall intensity in that event, the pattern of antecedent conditions prior to that event and various thresholds in the system response, including connectivity thresholds and the fill-spill hypoth- esis (e.g. Spence and Woo, 2003; Tromp-van Meerveld and McDonnell, 2006; Hopp and McDonnell, 2009;

Jencso et al., 2010; Tetzlaff et al., 2014), the generation of preferential flows (e.g. Beven and Germann, 1982, 2013; Beven, 2010; Chappell, 2010) and initiation of surface runoff (e.g. Kidron et al., 1999; Gomi et al., 2008;

Marshall et al., 2014). Identi fication of the controls on the catchment nonlinearity is also made more dif ficult by

uncertainties in observations of both the forcing data and the catchment discharge and evapotranspiration output, and the lack of observation techniques for the internal states of the system at hillslope and catchment scales.

Acknowledging the nonlinearity, however, brings us into the realm of developments in theoretical nonlinear dynamics and consequent inferences about chaotic behaviour, self-organization and the concept of the attractor. There have been many examples of the application of the principles of nonlinear dynamics to environmental systems; indeed one of the first demon- strations of chaos was in a simple three equation model of the atmosphere by Lorenz (1963, 1969). Since then, there have been many expositions of the principles of nonlinear dynamics in both the technical and popular science literature (e.g. Gleick, 1988; Lorenz, 1993; Thompson and Stewart, 2002; Smith, 2007; Strogatz, 2014).

Applications in hydrology have been reviewed by Sivakumar (2000, 2008, 2009, 2012) and have included the application of nonlinear dynamics to turbulence (Ruelle and Takens, 1971; Bai-Lin, 1983), to the fractal properties of rainfall fields (Schertzer and Lovejoy, 1996), concepts of self organization in drainage networks and

*Correspondence to: Keith Beven, Lancaster Environment Centre, Lancaster University, Lancaster LA1 4YQ, UK.

E-mail: k.beven@lancaster.ac.uk

Published online 21 November 2015 in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1002/hyp.10699

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

vegetation stands (e.g. Maritan et al., 1996; Rinaldo et al., 1996; Sherratt and Lord, 2007; Kleidon and Schymanski, 2008; Peterson et al., 2009; Schymanski et al., 2009;

Sheffer et al., 2013); dominant modes of response in the relationship between rainfall and runoff (Porporato and Ridol fi, 1996, 1997; Sivakumar, 2008; Sivakumar and Singh, 2012) and concepts of hydrological resilience (Peterson et al., 2012).

In this paper we want to look at interpreting the response of a catchment to time series of rainfalls and evapotrans- piration forcing as an example of the trajectory of a nonlinear system in phase space. In particular, we intend to examine the tendency for the catchment behaviour to have a constrained basin of attraction, and to explore the organization of the system trajectories in representations of the phase space. In doing so we will use the outputs of a recent model of the Gårdsjön catchment in Sweden, which has been shown to provide good fits to both discharge and tracer data (Davies et al., 2013) and can therefore be used to examine storage/residence time trajectories as well as storage/ flow trajectories (McDonnell and Beven, 2014).

While this is a virtual reality, it has the advantage that everything is known about the simulation (including the tracking of all particles of water through the catchment) such that inferences about the behaviour of the catchment are not in fluenced by either aleatory or epistemic uncertainties in either input or output variables (see Beven, 2012, 2015b for a discussion).

The MIPs model is, however, computationally intensive to run, even for this small catchment. One aim of this work is, therefore, to de fine the range of behaviours in both discharge and residence time responses that would need to be reproduced by any simple emulator model or closure scheme as discussed in Beven (2006).

CATCHMENT RESPONSES AND THE CONCEPT OF AN ATTRACTOR

In the theory of nonlinear dynamic systems an attractor is a set of numerical properties towards which the system evolves from a variety of initial conditions. For a system with finite dimensions, the trajectories of the system can be represented as an evolving N dimensional vector, with N as the number of dimensions de fining the phase space of the system. The attractor can be a point, curve or manifold in the space. It can also be fractal, in which case it is generally referred to as a strange attractor, a term coined by Ruelle and Takens (1971). Chaotic solutions of the Lorenz system are examples of a strange attractor with the well-known ‘butterfly wing’ form of trajectories.

Closed chaotic systems characteristically are associated with strange attractors, exhibiting bifurcations that are

sensitive to the initial starting point, but non-chaotic strange attractors also exist (Grebogi and Pelikan, 1984).

Where an attractor exists, then it is associated with a basin of attraction which is the area of the phase space within which points with arbitrary initial conditions will evolve towards the attractor. Where there are multiple attractors, then the trajectory will depend on which basin of attraction the initial conditions lie in.

Much of the theory of nonlinear dynamic systems is concerned with systems for which the trajectories evolve without external forcing. Catchments are open nonlinear systems for which the external forcings of rainfall (or other water inputs) and energy (driving evapotranspira- tion) are all important, and consistently force the system away from the trivial point attractor of zero active storage and zero fluxes, the point of maximum entropy for a storage system of this type. Catchments are also dissipative in the sense that rates of change of storage or energy in the system are less than the supply rate.

Dissipative systems can also exhibit chaotic divergences, multiple modes of organization and long-range correla- tions, as has been shown, for example, for turbulence.

Clearly all possible initial conditions and potential storage conditions in the catchment are within the basin of attraction of the zero flux end point. The issue then is whether the trajectories are structured in ways that might be useful and how the system moves towards the point attractor, especially in situations of frequent forcing as commonly the case in humid catchments. Concepts such as the principle of maximum entropy production, with its foundation in thermodynamics, have been invoked to understand the organization and behaviour of forced systems (e.g. Kleidon et al., 2010; Zehe et al., 2010). As an additional complication, a number of authors have suggested that the forcing itself may have chaotic properties (e.g. Rodriguez-Iturbe et al., 1989; Sivakumar, 1999; Sivakumar et al., 1999), which could result in catchment discharges also having chaotic properties (Jayawardena and Lai, 1994; Porporato and Ridol fi, 1996, 1997; Islam and Sivakumar, 2002).

Methods of studying nonlinear dynamic systems have

also been developed to assess the dimensionality of the

dynamic response (for examples in hydrology see

Porporato and Ridol fi, 1996, 1997; Sivakumar et al.,

2002, 2007; Sivakumar, 2005). These methods are mostly

based on correlation analysis and the derivation of

Lyaponov exponents from time series data of the system

response variables to determine the effective dimensional-

ity of these very high dimensional systems. Such methods

have a strong similarity with methods in stochastic system,

analysis and the interpretation of hydrological time series

as the result of deterministic chaotic dynamics is not

without its criticisms. In particular there has been some

contention about whether the data are adequate to identify

low dimensional chaos (e.g. Schertzer et al., 2002;

Sivakumar et al., 2002; Sivakumar, 2005) or whether a stochastic description is more appropriate than inferring low dimensional chaos (Koutsoyiannis and Pachakis, 1996; Koutsoyiannis, 2006, 2011). Here we explore the trajectory of a small catchment subject to external forcing, in terms of both flows and residence time distributions, to try and learn more about the structure of the response that might guide the formulation of new types of hillslope and catchment process representations (Beven, 2006).

THE CASE STUDY: MIPS APPLIED TO THE GÅRDSJÖN CATCHMENT

The Multiple Interacting Pathways (MIPs) model is a numerical hydrological model based on random particle tracking techniques that provides coherent simulation of plot to small catchment scale water flow and transport (Davies et al., 2012, 2013, 2015). This allows the simulation of discharges, storages, water age including input, output and storage residence times, and water origin

—making MIPs a good tool for exploring non-linearities and non-stationarities in catchment dynamics for both flow and transport. The MIPs model also provides a solution which directly represents a continuum of flow pathways in the subsurface, in contrast to many solutions which use assumptions of homogeneity or dual-domain representa- tions (see Davies et al., 2013 for a fuller discussion).

This model combines three concepts to characterize water flow and transport in a manner that implicitly represents heterogeneities in the subsurface. These are: random particle tracking for both flow and transport, pathway velocity distributions and pathway exchanges. Water in the catchment is represented as discrete packets, or ‘particles’

of water which each represents volumes of water moving at similar velocities. They enter the catchment via rainfall, can be lost from storage by evapotranspiration or can flow

towards the outlet according to mechanistic equations that are consistent with Darcy ’s law in the saturated zone but which allow preferential flows and bypassing.

However, to represent water flow within a subsurface that includes a multitude of pathways and velocities rather than assuming a mean velocity through an homogenous medium for which Darcy ’s equations apply, particle velocity distributions are used which integrate to provide a mean velocity. This distribution then is a representation of all the possible pathways in which water could travel. There is much flexibility in defining this velocity distribution (e.g. it could be asymmetrical and multi-modal); however, it must always be consistent with the de finition of the hydraulic conductivity and saturated porosity across the soil pro file in order to satisfy the constraint of integration to the Darcian mean velocity under saturated conditions.

The third element of the MIPs concept is pathway exchanges. Transition probabilities, which may be condi- tioned on soil moisture or other system states, are used to represent the exchange of particles between pathways, i.e.

between different velocities within the velocity distribution.

In this paper the MIPs model has been applied to the small G1 Gårdsjön experimental catchment located on the West coast of Sweden. The G1 catchment has an area of 6300 m ² , and it varies in altitude between 123 and 143 m.

a.s.l. Its vegetation is dominated by Norway Spruce with some Scotch Pine, and it has a podzolic soil underlain by gneissic granodiorite. The long-term mean annual pre- cipitation (between 1971 and 2001) is 1065 mm, and mean annual temperature is 6 °C.

Previously, the MIPs model has been successfully applied to field observations from the G1 catchment, reproducing a plot scale tracer experiment (Nyberg et al., 1999; Davies and Beven, 2012) and catchment scale isotope experiment data (Bishop et al., 2011; Davies et al., 2013). The methodology for applying MIPs to a catchment is described in detail in Davies et al. (2013), and the same final model is applied here. The most salient feature of modelling the G1

Figure 1. Exponential function used to describe saturated conductivity (K

) decline with depth which is characteristic of the soils within the G1 catchment,

and the cumulative velocity distribution of particles at depths 0.5 m and 0 m. Particle velocities at any given depth are distributed exponentially around the

mean velocity (as determined by K

and the porosity) so that the majority of particles are travelling slowly, with exponentially fewer faster particles

catchment is characterizing a consistent hydraulic conductivity –porosity–pathway velocity distribution relationship.. The saturated hydraulic conductivity at this site is characterized by a high near-surface transmis- sivity (as described by Nyberg, 1995; and Bishop et al., 2011), which rapidly declines with pro file depth, as macroporosity and the matrix porosity of the soil decreases.

As such, a non-linear exponential function is used to characterize this transmissivity pro file within the model as shown in Figure 1, while at each level particle velocities are chosen from an exponential distribution. An assumption that the faster velocities nearer to the surface are a result of increased porosity allows a simple relationships between porosity, transmissivity and velocity distributions to be developed (see Davies et al., 2011, 2013).

In the MIPs model, only the particle velocities are speci fied explicitly. Unlike in differential equation models, celerities controlling the flow response cannot be derived directly (see McDonnell and Beven, 2014, for an extended discussion of the velocity/celerity issue).

Celerity effects occur in MIPs, however, instigated by the

movement of particles in relation to local storage and the rise and fall of the water table on different parts of the slope. In particular as the catchment becomes wetter and water tables rise towards the soil surface, more particles are incorporated into the downslope fluxes causing the hydrograph to rise. The model can also allow for return flow and saturation excess runoff to occur once the water table reaches the surface. Davies et al. (2013) show how both celerity-controlled hydrographs and velocity- controlled residence time responses can be reproduced reasonably well at the Gårdsjön catchment.

Here, 30 years of rain gauge data (Temnerud et al., 2014) are used to drive the Gårdsjön MIPs model in order to provide a long-term output for phase analysis. Rainfall inputs are known on a daily basis, are assumed to be spread evenly over the catchment and are spread evenly over the day in the absence of further information on timings.

Evapotranspiration was assumed to take on a seasonal sinusoidal form with maximum of 3 mm/day and winter minimum of zero. The catchment is initialized from dry, and consequently the bulk of the analysis, where storage

Figure 2. a: Discharge phase plots with varying time step for a 1-year period. dt is 0.02 days (50 time steps per day). Plots give phase results with Q summed over 1, 2, 5 10, 25 and 50 time steps (0.02 to 1 day) each with a time lag of n*dt. b: Poincaré section plots in {v, w} dimensions of trajectories in

the Q

, Q

, Q

space for the same integral time steps

residence times are concerned, is concentrated on the latter 20 years. The 10-year run-in period gives suf ficient time for higher values of the storage residence times to develop.

Maximum storage times witnessed within the simulation are approximately 8 years, with 95% of the stored water having storage times less than 4 years.

The particle tracking approach allows the residence time of particles entering the system to be analysed. Thus residence times in the catchment can be determined for the water in each increment of input (input residence time), for all the storage in the catchment at any given time (storage residence time) and for any increment of output from the catchment in a time step (output residence times). These residence times are, except under special conditions, different from each other (e.g. Rinaldo et al., 2011) and will vary with the sequence of wetting and drying of the catchment.

ANALYSIS AND RESULTS

There are a number of different ways of representing the dynamic response of the catchment in terms of phase plots with dimensions representing the controls on the system

response. One type of phase plot to represent the response is to plot the output from the catchment at time t + 1 as a function of the out flow at the previous time (Takens, 1981;

Porporato and Ridol fi, 1996, 1997; Sivakumar, 2008).

Given the nonlinearity of the response, the relationship will vary depending on the length of the time step (see Figure 2).

The form of these plots depends heavily on the forcing from time step to time step. Where there is no forcing then there will be a tendency towards Q(t + 1) being a fraction of Q(t) (the classical recession constant concept in hydrological analysis). The fraction will be smaller for longer time steps, so that the slope of the relationships should get lower. Despite the nonlinearity of the soil transmissivity function in the Gårdsjön case (see Figure 1), the relationship stays rather linear in these plots, albeit being obscured by the effects of greater volumes of forcing at the longer time steps. This causes the hysteresis in the trajectories to increase with longer time steps. Such hysteresis is expected from consideration of the velocities and celerities in the catchment (Beven, 1989; McDonnell and Beven, 2014; Davies and Beven, 2015).

These modelled relationships will however be affected by the way the daily inputs have been distributed evenly

Figure 3. Hydrographs and storage –discharge/storage–mean output residence time relationships for two similar event periods with differing antecedent conditions. Outputs are given on a 5*dt basis (0.1 day timestep). The left column (plots in black) gives events between September and October 1980 where the catchment is relatively dry beforehand, and the right (plots in grey) show events between January and February 1982, where the catchment is relatively wet. It can be seen from the hydrographs that the catchment is more responsive under wet antecedent conditions. The storage –discharge relationships (Q vs S) show that there is more hysteresis under dry conditions, as opposed to wet. The storage –mean output residence time relationships

also display hysteresis

Figure 4. Hydrographs and storage –discharge/storage–mean output residence time relationships for an event period during October–November 2008.

Outputs are given on a 5*dt basis (0.1 day timestep). This sequence demonstrates how the hysteresis in the storage–discharge (Q vs S) and storage–

mean output residence time ( τ

vs S) relationships evolve under a sequence of events that takes the catchment from relatively dry conditions to relatively wet

Figure 5. Outputs for years 1989 to 2009. Top row: flow difference versus storage versus rainfall for days where rainfall > 0. Bottom row: flow

difference versus storage versus ET for days where rainfall = 0

throughout wet days, lacking further information about the sub-daily timing. There is an expectation that the hysteresis will be greater under dry antecedent conditions than for a similar event under wet antecedent conditions.

This is shown in Figure 3 which contrasts two similar events with quite different antecedent conditions and Figure 4 which shows a sequence of similar events during a wetting up sequence following a dry period.

There is hysteresis in both the storage –discharge relationship and in the storage –output–mean residence time relationship. Some small fluctuations resulting from the stochastic nature of the MIPs model can be seen in both the hydrographs, and more importantly in the fluctuations in mean residence times of the outflow discharge from time-step to time-step. This variability could be reduced by the use of smaller particles or averaging over multiple runs of the model at the expense of signi ficantly more computing time.

An alternative representation of the same data can be achieved by plotting two-dimensional Poincaré sections of trajectory intersections in three dimensions of fluxes or residence times at the current time, at time t
1 and time t
2. The relevant transformations to the Poincaré section

coordinates {v, w} are given, for the case of discharge Q, by (see Porporato and Ridol fi, 2003).

v ¼ ffiffiffi 2 p

2 ð
Q _t
2 þ Q _t
1 Þ (1a)

w ¼ ffiffiffi 6 p

6 ð
Q _t
2 þ 2Q t
Q _t
1 Þ (1b) Poincaré sections for equivalent time steps taken from the simulated flow data for Gårdsjön are shown in Figure 2B. The third quadrant of these plots (where Q t > Q t
1 > Q t
2 represent recession periods and for time steps above about 10 h) shows rather linear behaviour. It is generally accepted that there is a tendency towards more linear recession characteristics in catchment responses as catchment scale increases, but here the responses are for a small catchment with a strongly nonlinear soil transmissivity pro file. The

Figure 6. Outputs for years 1989 to 2009. Flow difference versus mean storage residence time versus storage for days where rainfall > 0

greater scatter in the first quadrant (generally lower flows) and at shorter time steps suggests higher order dynamics and stochastic noise in the response. Various other methods to explore the order of the nonlinearity and relative importance of stochasticity are reviewed by Porporato and Ridol fi (2003) and Sivakumar et al.

(2007).

Here we have chosen to use a phase space using summary hydrological variables. There is an expectation that for any given event or time step, the hydrological response will be a function of the input to the catchment and the antecedent wetness of the catchment. In particular, we are interested in the differing responses of the flow and transport in these dimensions, and the way in which the difference in celerities and velocities might lead to hysteresis in the response (Beven, 2006; McDonnell and Beven, 2014; Figures 3, 4). A representation of catchment storage –discharge–transport non-linearities over a much longer simulation period can be seen in the plots of Figures 5 and 6, where the catchment trajectories are shown in terms of storage, forcing inputs and the change in flow for time steps of 1 day. Normally, of course, it is dif ficult to make an accurate assessment of

storage in the catchment system because storage cannot be measured directly and inferred values from tracking the water balance are subject to important uncertainties in rainfalls, evapotranspiration and discharge observations.

In the MIPs model it is a simple matter of keeping a count of the particles in the system at any one time. Because each particle is labelled with the time at which it entered the catchment, this also gives a current estimate of the residence time distribution of water in the catchment.

Travel times can also be inferred for both increments of the rainfall inputs, and increments of the predicted discharges from the catchment.

In Figure 5 the top figures are for rainfall forcing, and the bottom for evapotranspiration forcing for days when zero rainfall was recorded. In both cases, flow differences get smaller as storage decreases, as expected. The wide range of responses for a given rainfall, conditional on the antecedent storage, can be clearly seen, as well as the tendency for the range of responses to increase as the total storage in the catchment increases. Interestingly, even some days of high rainfall can show negative changes in discharge (perhaps because the hydrograph has already peaked). In Figure 6, which relates the discharge response

Figure 7. Outputs for years 1989 to 2009. Flow difference versus storage residence time quantiles (5

in blue, 50

in cyan and 95

in red) versus storage

for days where rainfall > 0

to the changes in storage and mean storage residence time in the catchment, it is also clear that mean storage residence time is a hysteretic function of storage, showing a range of behaviours within the bounds of wetting and drying envelope curves. Given the nature of the forcing at this site, these cannot really be considered to be the primary wetting and primary drying curves that are commonly assumed in soil physics. Given more extreme wet and dry conditions than included within the 20 years of data at Gårdsjön it is still possible that the envelope would be extended. However, there is no evidence in these plots that there are multiple modes of behaviour in the dynamics.

Figure 6 presents only the mean residence times for storage, but the mean hides a broad distribution of residence times for individual particles moving at different Lagrangian velocities through the system. Figure 7 breaks down the same information in the form of the percentiles of the storage residence time distribution for different values of total storage in the system. This figure clearly shows the non-stationarity of storage residence time distributions. It can be seen that rapid transit of water through the system (represented by the 5 ^th percentile) can occur under both wet and dry conditions, while long residence times can be much much longer. Here those

longer residence times are the product of the assumptions made about velocity characteristics within the soil. These were shown to reproduce the tracer behaviour at Gårdsjön reasonably well but only over a period of months (Davies et al., 2013). Thus, it has to be expected that the longer residence times in such a system might be poorly constrained by short periods of calibration data (see also Beven, 2015a). The time series of these quantiles show some interesting behaviours analogous to a form of fill and spill hypothesis (e.g. Tromp-van Meerveld and McDonnell, 2006). The 95% quantile in particular shows a linear increase with time over relatively dry periods, until the system is wetted suf ficiently to displace some of the older particles and cause a rapid drop in the value of this quantile.

As soon as it is allowed that residence time distributions are varying through time with catchment wetting and drying then it is necessary to distinguish between residence times for the storage, and the residence times or travel times for an increment of input or an increment of output. These are related, but different (e.g. Rinaldo et al., 2011). This can be seen in Figure 8 which shows the equivalent mean residence time information for different discharge outputs at a given total storage. The shapes of the envelope curves are quite different, and, in general, the mean residence

Figure 8. Outputs for years 1989 to 2009. Flow difference versus mean output residence time versus storage for days where Q > 0

times are shorter for the discharges than for the storage.

Mean residence times under dry conditions can be small (although given the generally low discharges, this will represent only a small number of particles). Figure 9 shows the equivalent quantiles for discharge output residence times. These are much more constrained than for the storage in the system and show a clearer response to increase in storage.

DISCUSSION

Catchment systems are forced dissipative dynamic systems subject to distributions of input events that will differ in their responses depending on the antecedent pattern of storage in the catchment. The events vary in their pattern of intensity in time, duration and total volume. In general we expect the distribution of rainfall volumes to be skewed towards smaller events, with the largest events having low probability. It is often assumed that a larger event might always happen so that the range of catchment responses might not be encapsulated by

even a long historical record. Thus in considering the results from the 20-year run of data for the Gårdsjön catchment, it should be remembered that a larger forcing might occur in future, resulting in an excursion from the behaviours plotted in the previous section.

The different controls on velocities and celerities complicate the relationships between storage, discharge and residence or travel times in the catchment. The effects of a forcing event, distributed over the catchment area, take time to propagate down the hillslopes and affect the outputs. Discharges will be primarily controlled by the celerities propagating pressures waves in the saturated zone and preferential flow pathways. Travel times will be controlled by the distribution of velocities in both unsaturated and saturated zones, including the generation of surface runoff (in the MIPs model of Gårdsjön this occurs only as a saturation excess process because surface soil conductivities are high).

It is clear from these long simulations however that while the basin of attraction of the response is largely constrained by the magnitude and sequences of the events in the rainfall records, within that basin there is no clear

Figure 9. Outputs for years 1989 to 2009. Flow difference versus output residence time quantiles (5th in blue, 50th in cyan and 95th in red) versus

storage for days where Q > 0

structure, just the potential for a space filling of points under different conditions of rainfall, evapotranspiration and storage. The only de finite attractor will be the point [0,0] in the storage- flow space where given the conditions of the simulation there may be storage above the undulating bedrock surface that may be depleted only after the discharge reaches zero. The behaviour of the residence times is interesting in that for storage residence times tend to increase with decreasing storage with the approach to the point attractor (Figure 6), but output residence times can show a wide range of values for a given storage at lower discharges with low values when storage is low (Figure 8). This is perhaps a result of the lack of capillarity effects in the MIPs model, allowing some faster flow pathways to remain active even under dry conditions, albeit that such very dry conditions are rare in the Gårdsjön climate.

This is a reminder that all the results presented here are model responses, not the real catchment. This means that there are no concerns about uncertainties in the input data, observed outputs or catchment characteristics, but also means that the representation of flows and travel times are necessarily approximate. As noted earlier, the MIPs model has previously given good simulations of both flow and residence times for the Gårdsjön catchment, but clearly there is still scope for the model to differ from the real response so that the inferences must be considered conditional on both the model structure and the treatment of the inputs and on the speci fic (small) scale of the catchment. Different topographies, slope lengths, eleva- tion differences, and convergence/divergence, different climates and different soil/bedrock characteristics would change the detailed nature of the hysteresis and characteristics of the basin of attraction, but perhaps not the overall conclusions of the study.

CONCLUSIONS

One of the expectations in setting out to produce this paper was that there would be limited structure in the basin of attraction of a forced catchment hydrology simulator (see also the discussion in Beven, 2015a). This we have shown.

As noted above, the catchment is a dissipative system that might require a strong feedback mechanism to show multiple basins of attraction (as, for example, in the studies of Dent et al., 2002; Ridol fi et al., 2006; Peterson et al., 2009; or the potential for feedbacks between soil moisture and rainfalls at much larger scales, e.g. Entekhabi et al., 1992). No such feedback mechanism is incorporated into the MIPs catchment representation used here.

There are many other forms of analysis that could have been carried out on this data set to explore the complexity of the dynamics (e.g. Poporato and Ridol fi, 2003; and

Sivakumar et al., 2007), but here the main interest has been to explore the possibility of any structural forms that might need to be represented in an emulator model that could re flect the scale-dependent effects of celerities and velocities in producing the hysteresis of discharge fluxes and non-stationary residence times (Beven, 2006). The basin of attraction behaviour presented in this paper will provide a useful test of any such emulator (see also the emulation of hysteresis of a hydraulic flood model reported in Beven et al., 2009).

There may be other ways of using summary variables in the representation of the complex catchment dynamics.

In refereeing this paper, Erwin Zehe (pers. comm.) suggested that knowledge of the particle positions and velocities might allow both kinetic and potential energy in the system to be integrated over all the pathways for particles in the catchment at any time. This would also implicitly account for the energy losses associated with individual pathways. Rainfall forcing could be treated as an energy gain, and evapotranspiration could be treated in terms of an energy loss (capillarity is neglected in this version of the model). This would be an interesting perspective on the dynamics that might be worth exploring in future.

ACKNOWLEDGEMENTS

This work was funded by the UK ’s Natural Environment Research Council (NERC) under grant reference NE/G017123/1. We are grateful to Kevin Bishop, Allan Rodhe and Lars Nyberg who were involved in the original experimental work at Gårdsjön and worked with us on the application of MIPs there; and to Erwin Zehe and other anonymous referees whose comments have significantly improved the presentation of this work.

REFERENCES

Bai-Lin HAO. 1983. Bifurcation, chaos, strange attractor, turbulence and all that —on intrinsic stochasticity in deterministic systems. Progress in Physics 3: 329–416.

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

191 –219.

Beven KJ. 2006. The Holy Grail of Scienti fic Hydrology: 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 Processes 24(12): 1537–1547.

Beven KJ. 2012. 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. 2015a. What we see now: event-persistence in predicting the responses of hydro-eco-geomorphological systems? Ecological Model- ling 298: 4–15. DOI:10.1016/j.ecolmodel.2014.07.019.

Beven KJ. 2015b. EGU Leonardo Lecture: Facets of Hydrology —epistemic

error, non-stationarity, likelihood, hypothesis testing, and communica-

tion. Hydrol. Sci. J. DOI:10.1080/02626667.2015.1031761.

Beven KJ, Germann PF. 1982. Macropores and water flow in soils. Water Resources Research 18(5): 1311–1325.

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

wrcr.20156.

Beven KJ, Young PC, Leedal DT, Romanowicz R. 2009. Computation- ally ef ficient flood water level prediction (with uncertainty). In Flood Risk Management: Research and Practice, Samuels P, Huntingdon S, Allsop W, Harrop J (eds). CRC Press: Balkema. ISBN:978-0-415- 48507-4.

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. Hydrol. Process. 25: 3950–3959. DOI:10.1002/

hyp.8355.

Chappell NA. 2010. Soil pipe distribution and hydrological functioning within the humid tropics: a synthesis. Hydrological processes 24(12):

1567 –1581.

Davies J, Beven KJ. 2012. Comparison of a multiple interacting pathways model with a classical kinematic wave subsurface flow solution.

Hydrol. Sci. J. 57: 203–216. DOI:10.1080/02626667.2011.645476.

Davies J, Beven KJ. 2015. Hysteresis and scale in catchment storage, flow, and transport. Hydrol. Process. DOI:10.1002/hyp.10511.

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, Rodhe A, Nyberg L, Bishop K. 2013. Integrated modelling of flow and residence times at the catchment scale with multiple interacting pathways. Water Resour. Res. 49(8): 4738–4750.

DOI:10.1002/wrcr.20377.

Dent CL, Cumming GS, Carpenter SR. 2002. Multiple states in river and lake ecosystems. Philosophical Transactions of the Royal Society of London B: Biological Sciences 357(1421): 635–645.

Entekhabi D, Rodriguez-Iturbe I, Bras RL. 1992. Variability in large-scale water balance with land surface –atmosphere interaction. Journal of Climate 5(8): 798–813.

Gleick J. 1988. Chaos: Making a New Science. ISBN:0-14-009250-1.

Gomi T, Sidle RC, Ueno M, Miyata S, Kosugi KI. 2008. Characteristics of overland flow generation on steep forested hillslopes of central Japan.

Journal of Hydrology 361(3): 275–290.

Grebogi O, Pelikan Y. 1984. Strange attractors that are not chaotic.

Physica D 13: 261–268. DOI:10.1016/0167-2789(84)90282-3.

Hopp L, McDonnell JJ. 2009. Connectivity at the hillslope scale:

identifying interactions between storm size, bedrock permeability, slope angle and soil depth. Journal of Hydrology 376(3): 378–391.

Islam MN, Sivakumar B. 2002. Characterization and prediction of runoff dynamics: a nonlinear dynamical view. Advances in water resources 25(2): 179–190.

Jayawardena AW, Lai F. 1994. Analysis and prediction of chaos in rainfall and stream flow time series. J. Hydrol. 153: 23–52.

Jencso KG, McGlynn BL, Gooseff MN, Bencala KE, Wondzell SM.

2010. Hillslope hydrologic connectivity controls riparian groundwater turnover: implications of catchment structure for riparian buffering and stream water sources. Water Resources Research 46(10): W10424.

DOI:10.1029/2009WR008818.

Kidron GJ, Yaalon DH, Vonshak A. 1999. Two causes for runoff initiation on microbiotic crusts: hydrophobicity and pore clogging. Soil Science 164(1): 18–27.

Kleidon A, Schymanski S. 2008. Thermodynamics and optimality of the water budget on land: a review. Geophys. Res. Lett. 35, L20404.

DOI:10.1029/2008GL035393.

Kleidon A, Malhi Y, Cox PM. 2010. Maximum entropy production in environmental and ecological systems. Phil. Trans. R. Soc. B 365:

1297 –1302. DOI:10.1098/rstb.2010.0018.

Koutsoyiannis D. 2006. On the quest for chaotic attractors in hydrological processes. Hydrological Sciences Journal 51(6): 1065–1091.

Koutsoyiannis D. 2011. Hurst –Kolmogorov dynamics and uncertain- ty. Journal of the American Water Resources Association 47(3):

481 –495.

Koutsoyiannis D, Pachakis D. 1996. Deterministic chaos versus stochasticity in analysis and modeling of point rainfall series. J.

Geophys. Res. 101: 26,441–26,451.

Lorenz EN. 1963. Deterministic nonperiodic flow. Journal of the Atmospheric Sciences 20(2): 130–141. DOI:10.1175/1520-0469(1963) 020 <0130:DNF>2.0.CO;2.c

Lorenz EN. 1969. Three approaches to atmospheric predictability. Bulletin of the American Meteorological Society 50: 345–349.

Lorenz EN. 1993. The Essence of Chaos. University of Washington Press:

Seattle.

Maritan A, Rinaldo A, Rigon R, Giacometti A, Rodríguez-Iturbe I. 1996.

Scaling laws for river networks. Physical review E 53(2): 1510.

Marshall MR, Ballard CE, Frogbrook ZL, Solloway I, McIntyre N, Reynolds B, Wheater HS. 2014. The impact of rural land management changes on soil hydraulic properties and runoff processes: results from experimental plots in upland UK. Hydrolog- ical Processes 28(4): 2617–2629.

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 Resour. Res. 50: DOI:10.1002/

2013WR015141.

Nyberg L. 1995. Soil and groundwater distribution, flowpaths, and transit times in a small till catchment. PhD thesis, Department of Hydrology, Uppsala University, Sweden.

Nyberg L, Rodhe A, Bishop K. 1999. Water transit times and flow paths from two line injections of

H and

Cl in a microcatchment at Gårdsjön. Sweden. Hydrol. Process. 13: 1557–1575. DOI:10.1002/

(SICI)1099-1085(19990815)13:11 <1557::AID-HYP835>3.0.CO;2-S.

Peterson TJ, Argent RM, Western AW, Chiew FHS. 2009. Multiple stable states in hydrological models: an ecohydrological investigation. Water Resources Research 40(8): W03406. DOI:10.1029/2008WR006886.

Peterson TJ, Western AW, Argent RM. 2012. Analytical methods for ecosystem resilience: a hydrological investigation. Water Resources Research 48: W10531. DOI:10.1029/2012WR012150.

Porporato A, Ridol fi L. 1996. Clues to the existence of deterministic chaos in river flow. Int. J. Mod. Phys. B 10: 1821–1862.

Porporato A, Ridol fi L. 1997. Nonlinear analysis of river flow time sequences. Water Resour. Res. 33(6): 1353–1367.

Porporato A., Ridol fi L. 2003. Detecting determinism and nonlinearity in river- flow time series. Hydrological Sciences Journal 48(5), 763–780.

Ridol fi L, D’Odorico P, Laio F. 2006. Effect of vegetation-water table feedbacks on the stability and resilience of plant ecosystems. Water Resour. Res. 42: W01201. DOI:10.1029/2005WR004444.

Rinaldo A, Maritan A, Colaiori F, Flammini A, Rigon R. 1996.

Thermodynamics of fractal networks. Phys. Rev. Lett. 76: 3364–3367.

DOI:10.1103/PhysRev-Lett.76.3364.

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.

Rodriguez-Iturbe I, De Power FB, Shari fi MB, Georgakakos KP. 1989.

Chaos in rainfall. Water Resour. Res. 25(7): 1667–1675.

Ruelle D, Takens F. 1971. On the nature of turbulence. Communications in Mathematical Physics 20(3): 167–192. DOI:10.1007/bf01646553.

Schertzer D, Lovejoy S. 1996. Causal space-time multifractal processes:

predictability and forecasting of rain fields. Journal of Geophysical Research: Atmospheres (1984 –2012) 101(D21): 26333–26346.

Schertzer D, Tchiguirinskaia I, Lovejoy S, Hubert P, Bendjoudi H, Larchevêque M. 2002. Discussion of “Evidence of chaos in the rainfall- runoff process ” Which chaos in the rainfall-runoff process? Hydrolog- ical Sciences Journal 47(1): 139–148.

Schymanski SJ, Sivapalan M, Roderick ML, Hutley LB, Beringer J. 2009.

An optimality-based model of the dynamic feedbacks between natural vegetation and the water balance. Water Resources Research 45(1): W01412. DOI:10.1029/2008WR006841.

Sheffer E, Hardenberg J, Yizhaq H, Shachak M, Meron E. 2013. Emerged or imposed: a theory on the role of physical templates and self- organisation for vegetation patchiness. Ecology letters 16(2): 127–139.

Sherratt JA, Lord GJ. 2007. Nonlinear dynamics and pattern bifurcations in a model for vegetation stripes in semi-arid environments. Theoretical Population Biology 71(1): 1–11.

Sivakumar B. 1999. Rainfall dynamics at different temporal scales: a chaotic perspective. Hydrology and Earth System Sciences 5(4):

645 –652.

Sivakumar B. 2000. Chaos theory in hydrology: important issues and interpretations. Journal of Hydrology 227(2000): 1–20.

Sivakumar B. 2005. Correlation dimension estimation of hydrological series and data size requirement: myth and reality. Hydrological Sciences Journal 50(4): 591–603. DOI:10.1623/hysj.2005.50.4.591.

Sivakumar B. 2008. Dominant processes concept, model simpli fication and classification framework in catchment hydrology. Stochastic Environmental Research and Risk Assessment 22(6): 737–748.

Sivakumar B. 2009. Nonlinear dynamics and chaos in hydrologic systems:

latest developments and a look forward. Stochastic Environmental Research and Risk Assessment 23(7): 1027–1036.

Sivakumar B, Liong S-Y, Liaw C-Y, Phoon K-K. 1999. Singapore rainfall behavior: chaotic? J. Hydrol. Engng, ASCE 4(1): 38–48.

Sivakumar B, Persson M, Berndtsson R, Uvo CB. 2002. Is correlation dimension a reliable indicator of low-dimensional chaos in short hydrological time series? Water Resources Research 38(2): 3–1.

Sivakumar B, Jayawardena AW, Li WK. 2007. Hydrologic complexity and classi fication: a simple data reconstruction approach. Hydrol.

Process. 21: 2713–2728.

Sivakumar B., Singh V. P. 2012. Hydrologic system complexity and nonlinear dynamic concepts for a catchment classi fication framework.

Hydrology and Earth System Sciences 16(11), 4119–4131.

Smith LA. 2007. Chaos: A Very Short Introduction (Very Short Introductions, 159). Oxford University Press: Oxford, UK. ISBN:0192853783.

Spence C, Woo MK. 2003. Hydrology of subarctic Canadian shield: soil- filled valleys. Journal of Hydrology 279(1): 151–166.

Strogatz SH. 2014. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Studies in Nonlinearity), 2nd edition. Westview Press: Boulder.

Takens F. 1981. Detecting strange attractors in turbulence. In Dynamical Systems and Turbulence, Rand DA, Young LS (eds) (eds)., Lecture Notes in Mathematics 898. Springer-Verlag: Berlin, Germany; 366 –381.

Temnerud J, Andersson L, Arheimer B, Bishop K, Bringmark L, Lundin L, Löfgren S, Moldan F, Nilsson T, 2014. The CLEO-database: nine forested headwater streams in Sweden, 1977 –2009.

Tetzlaff D, Birkel C, Dick J, Geris J, Soulsby C. 2014. Storage dynamics in hydropedological units control hillslope connectivity, runoff generation, and the evolution of catchment transit time distributions.

Water resources research 50(2): 969–985.

Thompson JMT, Stewart HB. 2002. Nonlinear Dynamics and Chaos. John Wiley & Sons: New York.

Tromp-van Meerveld HJ, McDonnell JJ. 2006. Threshold relations in subsurface storm flow: 2. The fill and spill hypothesis. Water Resources Research 42(2): W02410.

Zehe E, Blume T, Blöschl G. 2010. The principle of ‘maximum energy

dissipation ’: a novel thermodynamic perspective on rapid water flow in

connected soil structures. Philosophical Transactions of the Royal

Society B: Biological Sciences 365(1545): 1377–1386.