Multiple hydrological stable states and the probability of climate variability causing a threshold crossing

Many physically based models of surface and groundwater hydrology are constructed without the possibility of multiple stable states for the same parameter set because they do not include positive feedbacks. For such a conceptualisation, at the cessation of a transient hydrological disturbance of any magnitude the model will return to the initial stable state, and thus show an infinite resilience. To highlight and challenge this assumption a numerical distributed eco-hydrological model (coupled hillslope Boussinesq-vertically lumped vadose zone) was developed in which qualitatively different steady state water table elevation and stream flow exist for the same parameter set. It is used herein to quantify catchment resilience to climatic disturbances. Resilience investigations have traditionally used equilibrium or limit cycle continuation analysis to quantify resilience by estimating the state space location and number of stable states (henceforth referred to as attractors) and thresholds (henceforth referred to as repellor). This method requires the use of constant or stable and smooth inter-annual cycles of climate forcing. As groundwater recharge is often dominated by episodic climate events, it is an open question as to whether the multiple hydrological attractors identified from continuation analysis still exist under stochastic climate forcing. Using the above model, this was investigated by undertaking simulations with 100 stochastic climate replicates of 118 years for a range of saturated hydraulic conductivity (k_smax) parameter values. It allowed assessment of which parameter values: i) lead to bimodal water table elevations emerging, indicating multiple attractors; and ii) quantification of the probability of a shift from a shallow to a deep water table attractor, and vice versa. To assess bimodal behaviour, for each k_smax value the mean depth to the water table over the final year of simulation was calculated and a histogram derived from the 100 simulations. Multiple water table depth attractors (see below) and modes did emerge but were dependent upon the spatially-averaged limiting infiltration rate parameter, I_o. For I_o of 20 mm/day (and k_smax from 0.05 to 0.7 m/day) the system shifts from the deep to shallow attractor at the first major rainfall event, thus resulting in only one clear mode. Reducing I_o by 50% and 75% resulted in increasingly clear bimodal water table elevations. The probability of an attractor shift was also calculated and found to be very dependent upon whether the initial state was the shallow or deep water table attractor, with a shift from the deep to the shallow attractor being more probable than the reverse. Overall, this work: i) strengthens the important theory of multiple hydrological attractors; ii) expands quantitative resilience concepts and methods to better incorporate disturbances; and iii) extends characterisation of which catchments are likely to have multiple attractors. (Graph Presented).