PERSiST : the precipitation , evapotranspiration and runo ff simulator for solute transport

Introduction Conclusions References Tables Figures


Introduction
Understanding the fate and transport of pollutants in surface waters is dependent on credible simulations of movement of water through the landscape.In principle, this means partitioning precipitation (P ) into runoff (R) and evapotranspiration (E ) so as to satisfy the relationship P = R − E .There are a wide range of approaches to simulating hillslope and catchment-scale water fluxes, depending on the purpose of the modeling exercise (Kampf and Burges, 2007).Numerous models have been developed for flood flow forecasting (Bergström and Singh, 1995;Vehviläuinen, 2007 inter alia) or as a tool to better understand the hydrology of well-instrumented catchments (Fenicia et al., Figures Back Close Full 2011; Kavetski and Fenicia, 2011;Hrachowitz et al., 2013).While these models have proven extremely useful for both operational and research hydrology, they have not necessarily been well suited to simulating solute transport.Hydrological modeling is often not an end unto itself, but is performed in order to gain understanding about pollutant fate and transport.A number of models have been developed for simulating pollutant transport and transformations in catchments and surface waters including HBV-NP (Andersson et al., 2005), NAM-Mike 11 (Andersen et al., 2006), HYPE (Lindström et al., 2012) and the INCA family of models (Whitehead et al., 1998;Wade et al., 2002).These are all semi-distributed bucket-type models in which a single perceptual model of the runoff generation process is used to represent the movement of water through a number of stores in the catchment.While these models work well in many cases, greater flexibility in the representation of the perceptual model of the runoff generation process would be desirable.For example, Bernal et al. (2004) observe that the perceptual model of runoff generation used in the INCA model, which is derived from a perceptual model based on observations of temperate catchments, is not ideal for simulating hydrochemistry in intermittent Mediterranean streams.
The approach taken to simulating water fluxes is strongly influenced by the perceptual model (Fenicia et al., 2011;Beven, 2012) of the system.That is to say, the modeler's preconceptions and prior identification of important processes and pathways influence the manner in which runoff generation process is represented in the model, thereby dictating model structure.In practice, this means that the use of models for simulating water fluxes is influenced by the subjective experiences of the model developer, including the types of system he or she is most familiar with, background knowledge and degree of mathematical sophistication.
A perceptual model of water fluxes in which the landscape is represented as one or more buckets, or generic containers that receive, store or transmit water, has a long history, starting with the Sacramento Soil Moisture Accounting model which was developed in the late 1960s (Burnash et al., 1973).While this approach has been criticized recently (Gupta et al., 2012) it is still widely used for both operation hydrological ap-Figures

Back Close
Full plications.Notably, the HBV (Bergström and Singh, 1995) and WSFS (Vehviläuinen, 2007) models are widely used for operational flood forecasting in Scandinavia and the IHACRES (Jakeman et al., 1990) model and its variants for operational purposes in Australia.Simple bucket-type models such as HYMOD (Wagener et al., 2001) are also a key component of hydrologic research programs developing new Bayesian techniques for model calibration (Vrugt et al., 2009) Because they have been developed for different regions, spatial scales and climatic conditions, most currently used rainfall-runoff models can be exceedingly difficult to apply successfully outside of the conditions for which they were designed.There is an implicit "one size fits all" approach in most rainfall runoff models where a single perceptual model of the runoff generation process is assumed to be applicable under all conditions.Clearly, this is not an ideal situation.One possible response to the "one size fits all" problem is the development of modular frameworks in which individual model components can be assembled that better represent the modeler's perceptual model of the runoff generation process.Two recent examples of such frameworks are the snowmelt-runoff modeling framework of Smith and Marshall (2010) and the SU-PERFLEX rainfall runoff modeling framework (Fenicia et al., 2011;Kavetski and Fenicia, 2011).The SUPERFLEX framework allows a modeler to develop multiple model structures by combining and connecting a variable number of stores so as to provide a lumped, catchment-scale representation of the runoff generation process.One clear advantage of these frameworks over traditional rainfall runoff models is that they provide a great deal of flexibility in model structure, thereby facilitating a more credible representation of underlying perceptual models.
Flexible modeling frameworks are needed in hydrology if progress is to be made in understanding and simulating the transport of solutes including nutrients and contaminants.Refsgaard (2007)  Unfortunately, the development of an entirely new model is not often an option, and the model user is forced to use someone else's perceptual model.This is not always ideal as it can introduce known but unavoidable uncertainties related to model structure.The approach presented here is similar to that advocated by the developers of the SUPERFLEX modeling framework (Kavetski and Fenicia, 2011;Fenicia et al., 2011) insofar as it helps to overcome this problem.One potential difference between the model presented here and the SUPERFLEX approach is that most published SUPERFLEX applications use a single, lumped representation of a catchment whereas the model framework presented here has been designed from the outset to be semi-distributed.
Here, we present a description and first application of PERSiST, the Precipitation, Evapotranspiration and Runoff Simulator for Solute Transport.PERSiST is a semidistributed bucket-type modeling framework which can be easily modified by model users to incorporate multiple perceptual models of the runoff generation process.Here, we describe the structure and assumptions of the PERSiST model, present a demonstration application to simulate flows at 8 sites along the main branch of the Thames River in the UK, provide a brief description of a simple Markov Chain Monte Carlo calibration strategy and give some information about model parameter sensitivity.

Model description
PERSiST is a watershed-scale hydrological model suitable for simulating runoff generation processes at a catchment to landscape scale.It is a conceptual, daily time- (v) abstraction and discharge from industrial sources including drinking water supply; (vi) improved simulation of biogeochemically important low-flow events, (vii) a full water balance and (viii) an ability to simulate age of water in different stores throughout the landscape.The model has been implemented as a series of first-difference equations and has been designed to be as simple as possible but no simpler.
By giving the model user an ability to represent different patterns of water storage and connectivity between stores in the landscape, PERSiST can be used to evaluate different perceptual models and explore the effects of model structural uncertainty on runoff prediction.The model has limited data requirements.It requires daily time series of air temperature and precipitation from one or more sites as driving data.Thus, the model can also be used for projecting possible future patterns of runoff by using downscaled temperature and precipitation time series from regional or global climate models.PERSiST is calibrated against stream flow measured at one or more points in a river.The necessary spatial data to run the model include descriptions of all relevant hydrological response unit types; subcatchments areas, the proportional coverage of different landscape types within each sub catchment; and reach (river or stream) information including length and average width.When available, additional data on abstraction and discharge volumes or stream flow velocity can aid in model calibration.Data on soil moisture or depth to groundwater will also aid in model calibration by providing opportunities for soft calibration sensu Seibert and McDonnell (2002).sub catchment is associated with a reach.Landscape units consist of one or more connected water stores.These stores can be conceptualized as buckets that receive inputs of water from the atmosphere, other buckets and potentially river water.Landscape units are analogous to hydrological response units (sensu Wade et al., 2001).

Land cover types and watersheds
There are some parameters in PERSiST which are specified by individual land cover types and are applicable across the entire watershed (Table 1).These include the landscape-scale snow threshold temperature (l 1 ) and the snowfall (l 2 ) and rainfall (l 3 ) multipliers.When the air temperature is below the snow threshold temperature, precipitation is assumed to fall as snow and accumulate in the snowpack.The depth of snowfall is calculated by multiplying the observed precipitation by the snowfall multiplier.When air temperature is above the snow threshold temperature, precipitation is assumed to fall as rain.The depth of rainfall is estimated by multiplying the observed depth of precipitation by the rainfall multiplier.The rate of snowmelt (mm d −1 ) is estimated as the lesser of (i) the difference between measured air temperature and snow threshold temperature ( • C) multiplied by the degree day melt factor (l 4 ; mm • C −1 d −1 ) and (ii) the depth of the snowpack.Actual evapotranspiration is calculated in a similar manner to that presented by Durand (2004).However, instead of using Penman potential evapotranspiration as the baseline, a degree day evapotranspiration parameter (l 5 ; mm • C −1 d −1 ) which defines the maximum rate at which evaporation (i.e.potential evapotranspiration) can occur when air temperatures are above the growing degree day threshold (l 6 ; • C) is used.
When air temperatures are below the growing degree day threshold, it is assumed that no evapotranspiration occurs.The actual rate of evapotranspiration can be limited by the moisture content of the soil, as described below in the section on bucket properties.
The land cover specific maximum possible rate of evapotranspiration (E (l ); mm d −1 ) is calculated using a temperature index based on the difference between observed air temperature and the growing degree day threshold ( ).
The actual rate of evapotranspiration can be less than the maximum potential rate, depending on the amount of moisture available in the bucket from which water is returning to the atmosphere.

Sub catchment and reach
A watershed is represented as one or more subcatchments.Within a watershed, the stream is divided into reaches and it is assumed that there is one sub catchment per reach (Fig. 1).So as to deal with differences in climate throughout the watershed, it is possible to associate a different temperature and precipitation time series with each sub catchment.PERSiST is able to represent both branched and main-stem stream topologies.The reach is assumed to represent the main stream channel within a sub catchment.
To run PERSiST, additional parameters are specified for subcatchments and reaches (Table 1).To deal with possible elevation gradients in precipitation dynamics, it is possible to specify snowfall (s 1 ) and rainfall (s 2 ) multipliers at the sub catchment level.The effective snowfall and rainfall multipliers are determined by multiplying the landscapescale and subcatchments-scale parameter values.The sub catchment area (s 3 ) and the proportional cover of each landscape unit type must also be specified.Reach parameters include length (r 1 ), width (r 2 ) and the parameters necessary to determine flow velocity (v) as a function of stream flow (Q). (2) Rates of water abstraction from and effluent input to individual reaches may be specified either as constant values (r 5 , r 6 ; m 3 s −1 ) or as time series of daily average flow values.Figures

Back Close
Full The effects of land use change on hydrology can be simulated in PERSiST by allowing the proportional cover of landscape unit types to vary over time.

Bucket
A landscape unit type is comprised of one or more stores of water.These stores are conceptualized as well mixed volumes of water held in buckets.Buckets can store water, return it to the atmosphere through evapotranspiration, transfer it to other buckets or to surface waters.Each bucket has the following properties (Fig. 3, Table 1).The depth of water in the bucket at time t is expressed as z t and has units of millimetres.The maximum depth of water (b 1 ) is the maximum depth of water, in millimetres, that can be held in a bucket.The retained water depth (b 2 ) is the depth of water below which water no longer freely drains.When water is below this depth, normal runoff controlled only be the depth of water in the bucket does not occur but evapotranspiration can continue and drought-related runoff (described below) can occur.Water drains from a bucket with a characteristic time constant (b 3 ) with units of days.The depth of water draining on day t is calculated as follows: Water can be returned to the atmosphere through evapotranspiration.While water returned to the atmosphere through evaporation is released from surfaces, water returned as transpiration may be derived from different depths, depending on the root structure of the vegetation community in the landscape unit.Thus, there is a need to simulate different rates of evapotranspiration from different buckets.Typically, the rate will be highest in a surface bucket as it will include evaporation, to unity within a landscape unit so the base degree day evapotranspiration (b 4 ) is consistent with the total evapotranspiration produced within a landscape unit.When the depth of water in a bucket is above the retained water depth, evapotranspiration occurs at the maximum rate.When the depth of water is below the retained water depth, the rate of evapotranspiration can be limited as follows: Changing values of b 5 adjust the rate at which evapotranspiration slows when the water store, represented as the depth of water in a bucket, is below the retained water depth.
A value of 0 indicates that the rate of evapotranspiration will be unchanged by soil moisture status while high values (10+) effectively stop evapotranspiration when there is no longer any freely draining water.
The amount of water that can be added to a bucket in any given time step is limited.It cannot exceed the infiltration rate (b 6 ) or the difference between the maximum and current depths of water (b 1 − z).Infiltration excess is that water which is prevented from percolating due to the infiltration capacity being exceeded and it is routed to the adjacent stream.Water that is prevented from percolating due to the maximum storage capacity of the receiving box being exceeded is referred to as saturation excess, and will be routed back to the quick box as described in the section on landscape units.
Low flow events can have a disproportionate importance for surface water solute chemistry.Small increases in flow, which have a negligible effect on the overall water budget, can transport high concentrations of solute to streams.This has been observed for the flushing of nitrate to the Thames (Jin et al., 2012) or organic carbon from boreal catchments (Ledesma et al., 2012) inter alia.In many hydrological models, rainfall during dry conditions is assumed to contribute only to recharging soil and groundwater.However, a small fraction of the rain may contribute to stream flow and solute transport, even when soils are very dry.This phenomenon is simulated in PERSiST using the drought runoff fraction (b 7 ).When the depth of water in a bucket is below the re-Figures

Back Close
Full tained water depth, the amount of water entering the bucket which contributes to runoff is estimated by multiplying the total input by the drought runoff fraction.Note that the default behavior is for all water inputs to contribute only to recharge when depth of water in the bucket is below the freely draining threshold (b 2 ).There are two special bucket types.Quick flow buckets simulate surface processes while bidirectional buckets can receive inputs of water from the river through infiltration or inundation.Quick flow buckets receive inputs of rainfall and snowmelt.Saturation excess flow generated by other buckets in a landscape type is routed through a quick bucket to the adjacent reach.Each landscape unit type must include quick flow buckets for which the relative area indices sum to unity.
The vast majority of water flows within a watershed are from land to surface waters.However, flows from rivers to the land can also occur.During flood conditions, a river can overflow its banks and inundate the surrounding land.It is also possible for water to infiltrate from a river to the surrounding land.Both of these phenomena can be biogeochemically important.Inundation can deposit large amounts of sediments and nutrients (Bayley, 1995) while infiltration can alter rates of nitrogen processing (Grischek et al., 1998).For either inundation or infiltration to occur, a bucket must be identified for bidirectional flows.It is only possible to have one inundation and one infiltration flow bucket in a landscape unit type.
Inundation can only be simulated for bidirectional quick flow buckets.Inundation is simulated when the depth of water in the reach (d ) exceeds the bucket specific inundation threshold (b 9 ).The volume of water inundating per unit time (V Inundate ; m 3 d −1 ) is estimated as the fraction of the total flow through the reach outflow which occurs at a depth exceeding the inundation threshold.A rectangular channel cross section and uniform flow velocity throughout the reach is assumed when calculating inundation.

Conclusions References
Tables Figures

Back Close
Full The depth of water inundating the bucket per unit time (z Inundate ; mm d −1 ) which corresponds to V Inundate is calculated by dividing the inundation volume by the relative area of the bucket in the relevant landscape cover unit all multiplied by the sub catchment area.
Inundating water is treated exactly the same way as other inputs to a quick flow bucket.
Infiltration from the stream to the land can be simulated in PERSiST by moving water from the reach to either a regular or quick flow bucket.The height of the water column in a bucket is determined by dividing the depth of water in a bucket (z) by its porosity (b 10 ).Infiltration occurs when the depth of water in the reach (d ) exceeds the height of the water column of the receiving bucket (d − r 7 > z/b 10 ), where r 7 is an offset to account for cases where the water depth for the stream and the riparian soil are measured against different reference levels.The volume of infiltrating water per unit time (V Infiltrate ; m 3 d −1 ) is calculated in a similar manner as for inundating water.After some algebraic rearrangement, V Infiltrate can be expressed as follows: The increase in water depth in the bucket per unit time due to infiltration (z infiltrate ; mm d −1 ) can be calculated as follows:

Landscape unit
A landscape (or hydrological response) unit type consists of one or more buckets linked together in a user-specified manner so as to represent the perceptual model of the 8646 Introduction

Conclusions References
Tables Figures

Back Close
Full runoff generation process.Water is routed from the landscape unit to the stream.Thus, there is no movement of water between landscape units but PERSiST is able to represent perceptual models, for example, of runoff generated from recharge and discharge areas by appropriate combinations of buckets.Each landscape unit must contain one or more quick buckets to receive inputs of precipitation.Flows of water between buckets are described using a square matrix (Table 2).Element (i , j ) of the square matrix represents the fraction of water leaving bucket i which is added to bucket j .Diagonal elements (i , i ) on the square matrix define the fraction of water leaving the bucket that is routed directly to the stream.Off-diagonal elements (i , j ; j > i ) in the upper quadrant define the fraction of water leaving the row bucket entering the column bucket.The values of cells in each row from the diagonal to the rightmost entry must sum to 1.
Elements of the square matrix below the diagonal are used to identify the quick bucket to which saturation excess flow can be routed.Figure 2 illustrates a simple landscape unit consisting of three buckets which are able to generate stream flow, the square matrix associated with this landscape unit is shown in Table 2.The direct runoff bucket functions as a quick bucket.It is able to receive inputs of precipitation from the atmosphere.Of the water leaving the direct runoff bucket, 10 % is routed directly to the reach and 90 % to the soil water bucket.Of the water leaving the soil water bucket, 40 % is routed to the reach and 60 % to the groundwater bucket.All water leaving the groundwater bucket is routed directly to the reach.The value of 1 in the square matrix below the diagonal in row 2 indicates that saturation excess runoff from the soil water bucket is routed to the direct runoff bucket.Note that saturation excess runoff can only be generated from buckets that are immediately adjacent to a quick bucket (e.g.saturation excess runoff from the groundwater bucket is not possible in the configuration presented here).Introduction

Conclusions References
Tables Figures

Back Close
Full The configuration shown in Fig. 2 is similar to the terrestrial hydrological representation used in the INCA model (Whitehead et al., 1998;Wade et al., 2002).This representation has proved successful in temperate (Whitehead et al., 1998;Wade et al., 2002), montane (Ranzini et al., 2007;Futter et al., 2009) and boreal (Rankinen et al., 2004) conditions.It has been less successful in Mediterranean conditions (Bernal et al., 2004).However, PERSiST can be set up to include alternative model structures proposed for the simulation of Mediterranean hydrology (e.g.Medici et al., 2008) or diverse montane catchments (Kavetski and Fenicia, 2011).
Riparian and upland areas can also be simulated using PERSiST.Representing watershed hydrology using a series of vertically stacked buckets may not be appropriate for simulating the hydrochemistry of riparian zones in forest (Löfgren et al., 2011) or agricultural (Stutter et al., 2009) watersheds.PERSiST also allows for horizontally stacked buckets, whereby riparian and upland can be simulated separately.When simulating riparian and upland areas, some additional parameters must be specified for the different buckets so as to ensure that hydrology is represented correctly.The relative area index (b 8 ) is used to describe the areal contribution of each bucket to a landscape unit type.For example, a landscape unit type consisting of upland and riparian areas might have relative areas of 0.9 and 0.1 respectively for upland and riparian soil water buckets.

Calculating fluxes of water in a landscape unit type
Fluxes of water through a landscape unit type are calculated in the same order as which buckets are identified in a landscape unit type square matrix.The first row in the square matrix should represent a quick flow bucket and the last row a bucket where water drains only to the stream, and not to any other bucket.All fluxes between buckets are calculated in units of m 3 d −1 for a representative 1 km 2 landscape unit type (level Figures First, evapotranspiration is calculated according to Eq. ( 5).Next, outputs from a bucket to other buckets are calculated.The volume of water transferred from bucket i to bucket j (V i ,j ) is calculated as follows.First, the default output volume is estimated by multiplying the appropriate value in the square matrix (m i ,j ) by the area of the depth of water (∆z i ; Eq. 4) able to leave the bucket all multiplied by relative bucket area (b 8 ).
For each transfer, a test is performed to see if the volume of water being transferred is greater than the empty volume of water in the receiving bucket.In the event that the empty volume in the receiving bucket (j ) is smaller than the volume of water leaving the source bucket (i ), the volume leaving the source bucket is reduced accordingly.
After each transfer, the depth of water in buckets i and j is adjusted accordingly.Third, if a bucket is a quick flow type, snowmelt, rainfall and inundation are added.When the temperature is warm enough for snowmelt, the depth of snowmelt (z snowmelt ) is calculated as follows: where l 4 represents the degree day snowmelt factor and (T Air − l 1 ) is the number of degrees above the snowmelt offset.Rainfall (z rain ), which only occurs when the air temperature is above l 1 + s 1 , is equal to the observed depth of precipitation (P ) multiplied by the water shed scale (l 3 ) and sub catchment scale (s 2 ) precipitation multipliers.
Saturation excess flow is simulated when the depth of water in a bucket exceeds the maximum possible water depth.All saturation excess flow must be routed to a quick bucket.Saturation excess runoff added to a quick bucket is treated the same way as other inputs from precipitation or inundation.8649 Figures

Back Close
Full

PERSiST time series
PERSiST generates time series of inputs, outputs and changes in storage for each bucket in each landscape unit type in every sub catchment (Table 3).At each time step, the depth of water in each bucket is recorded, as are all transfers between buckets.Total fluxes of water from subcatchments to reach are also recorded, as are the rates of stream flow and any infiltration or inundation.Estimates of atmospheric exchange (precipitation and evapotranspiration) and snowpack dynamics are also reported.By tracking the time spent by water in each bucket, it is possible to simulate age of water in the catchment.

INCA compatibility
One of the design criteria for PERSiST is to generate input data files for the INCA model.Currently, INCA requires the use of an external rainfall runoff model to generate time series of soil moisture deficits (SMD) and hydrologically effective rainfall (HER), (Whitehead et al., 1998;Wade et al., 2002).Typically, these time series are obtained from either the HBV model (Saelthun, 1996), WSFS (Vehviläuinen, 2007), MORECS (Hough and Jones, 1997) or a routine developed by Durand (2004).Sub catchment and watershed-scale estimates of SMD are produced in PERSiST based on average differences between depth of water in a bucket (z) and its maximum water holding capacity (b 1 ).
HER is an estimate of the precipitation entering a watershed which eventually contributes to runoff.In PERSiST, this can be estimated at the sub catchment scale as precipitation minus actual evapotranspiration.These calculations take into account the effect of soil moisture status on evapotranspiration rates.HER is estimated by working backwards through time series of simulated actual evapotranspiration and precipitation inputs.Starting from the last day in the simulation, evapotranspiration is accumulated and then the precipitation for that day is subtracted.If the value is less than 0 (i.e.precipitation is greater than the cumulative evapotranspiration deficit) then the differences Figures

Back Close
Full is recorded as HER for that day and the cumulative evapotranspiration set to 0. This is repeated until the start of the simulation is reached.The adequacy of this assumption is tested by comparing the total estimated HER to the total runoff.The SMD is defined as the sum of differences between the maximum capacity (b 1 ) and current depth of water in a series of user-specified buckets.There is a possibility to adjust this depth by a user-specified offset so as to obtain SMD time series with a minimum value of 0.

Model application
An application of PERSiST is presented for the Thames River in the UK.The Thames is a regionally important drinking water source, serving approximately 14 million people in the greater London region (Jin et al., 2012;Whitehead et al., 2013;Crossman et al., 2013a).The watershed has an area of ∼ 10 000 km 2 and the main river a length of 218 km.Land use in the watershed is predominantly agricultural in the upper reaches and becomes more urban near the outflow.Bedrock in the catchment is a predominantly permeable chalk, although there are regions underlain by more recent, less permeable clay deposits and less permeable sedimentary bedrock.Base flows throughout the watershed are high.
For the application presented here twelve landscape unit types were used.These represent forest, arable, pasture and urban land uses on chalk, clay and sedimentary bedrock.Land use areas were obtained by generalizing land cover data from CORINE2000.Areas underlain by chalk, clay and sedimentary bedrock were estimated using the British Geological Survey 1 : 62 500 map of surficial geology.The Thames was divided into 8 reaches based on reach areas from Jin et al. (2012).Reach boundaries were based on the location of flow sampling stations (Fig. 4).Sub catchment areas, reach length and proportional cover of different landscape unit types are shown in Table 4.The model was calibrated using a single time series of temperature and precipitation obtained from the UK Met Office (Crossman et al., 2013a).The time series Introduction

Conclusions References
Tables Figures

Back Close
Full was based on a synthesis of observations made at observing stations in the Thames watershed.
Each landscape unit type was simulated as three vertically stacked buckets representing direct runoff, soil water and groundwater.It was assumed that all precipitation entering the direct runoff bucket percolated to the soil water bucket.Water could leave the soil water bucket as runoff to the stream, percolation to the groundwater bucket or as saturation excess flow which was routed immediately to the stream.All water entering the groundwater bucket was assumed to flow to the stream.It was assumed that there were no losses to deep groundwater.
The model application was performed for the period 1 January 1999-31 December 2008.PERSiST was calibrated against observed stream flows collected at 8 locations along the main stem of the Thames River (Fig. 4).The model was calibrated by first manually adjusting bucket parameters so as to improve the fit between modeled and observed stream flow.In all cases, model fit was assessed using the Nash-Sutcliffe (NS) statistic (Nash and Sutcliffe, 1970).The manual calibration first attempted to obtain the best possible fit to the upper most reach.Once this was obtained, attempts were made to fit subsequent reaches.Manual calibrations continued until it was no longer possible to obtain obvious improvements in NS statistics.The subsequent automated calibration was based on maximizing the sum of NS statistics.Thus, it was possible to accept poorer model performance in one reach if the parameter set under evaluation led to better fits in other reaches.The final manual calibration was used as the starting point for a simple Markov Chain Monte Carlo (MCMC) exploration of parameter space.A total of 108 parameters were allowed to vary during the MCMC analysis (Table 5).These included time constants for quick, soil water and groundwater buckets; subcatchments and land cover type rain multipliers; ET related parameters and reach flow velocity parameters.Introduction

Conclusions References
Tables Figures

Back Close
Full

MCMC tool
The MCMC code was based on the Metropolis Hastings algorithm (Chib and Greenberg, 1995) and did not use Gibbs sampling (Smith and Roberts, 1993).The Monte Carlo procedure consisted of 400 instances of MCMC chains consisting of 2500 model runs, for a total of 1 million model runs.The MCMC sampler operated as follows: 1.The manual calibration was used as a basis for identifying credible maximum and minimum ranges for parameters which would be varied during MCMC analysis.Typically, the credible range was defined as ±20 % of the best performing parameter set value.The system is initialized using the parameter set from the manual calibration to specify the best goodness of fit and the best performing parameter set.In the application presented here, the goodness of fit is estimated as Σ(1 − NS) for modeled and observed flows in the 8 subcatchments.The statistic has a maximum value of 0 and a minimum of −∞.

2.
A random starting point is selected from the parameter space, and the goodness of fit recorded.If the goodness of fit is better than the best goodness of fit, the starting point is accepted and the best goodness of fit and best performing parameter set are updated.If the goodness of fit from the random starting point is worse than the best goodness of fit, the ratio between the two values is compared to a random number between 0 and 1.If the ratio exceeds the random threshold, the starting point is accepted, otherwise a new starting point is drawn and its goodness of fit evaluated.
3. A jump is defined which randomly perturbs the parameter values.
4. The goodness of fit of the parameter set obtained for the jump is assessed.If the fit is better than that obtained from the previous parameter set, the jump is repeated.If the fit is better than the best observed goodness of fit, the best goodness of fit and best performing parameter set are updated.This process is repeated until no Introduction

Conclusions References
Tables Figures

Back Close
Full further improvement in goodness of fit is obtained or until a jump would cross the maximum or minimum value identified in (1).If the jump would cross the maximum or minimum threshold, a new value is randomly selected within the range of that parameter.
5. If a jump does not lead to an improvement in goodness of fit, the jump may still be accepted if the ratio of new and old goodness of fit exceeds a random number between 0 and 1.If the jump is rejected, a counter is incremented.If the counter exceeds a user-specified threshold (50 in the application presented here), control returns to (2) and a new random starting point is defined.If the threshold is not exceeded, control returns to (3).
6.This process is repeated 2500 times and the best performing parameter set retained for further analysis.
This process was repeated 40 times to generate an ensemble of credible parameter sets.The credible parameter sets are not globally optimal, but are the result of the non-exhaustive MCMC exploration of parameter described above.

Results
The best performing parameter sets obtained during the MCMC analysis had NS statistics for individual reaches between 0.47 and 0.86 (Fig. 5).These parameter sets were re-run to generate ensembles of predicted values (Figs. 6 and 7).Parameter sensitivity was assessed by comparing the cumulative distribution of parameters (Fig. 8a, b and Table 6) to a rectangular distribution which would be indicative of parameter randomness.
There was a general trend towards better model fits lower in the catchment.The ensemble of best performing parameter sets were able to reproduce the flow dynamics in both upper (Fig. 6a, b) and lower (Fig. 7a, b) reaches of the Thames.Modelled high flows in the upper part of the catchment were often higher than observed, base flows Introduction

Conclusions References
Tables Figures
Groundwater time constants were skewed towards the lower end of the parameter range (Fig. 8) and degree day evapotranspiration rates towards the upper end of the parameter range.These land cover types are found primarily in the lower reaches of the catchment (Table 4).Reach specific parameters were only sensitive in the upper parts of the catchment.It is notable that while the distribution of groundwater time constants and evapotranspiration rates were strongly non-rectangular, they still spanned the whole parameter range.

Discussion
Here  , 7a) and peak flows were better simulated in the lower reaches.No flows higher than 70 m 3 s −1 were reported for the uppermost reach (Fig. 6b).While it is unlikely, it is not impossible that there are significant uncertainties in flows estimated at this site.Beven (2012) inter alia has noted that it can be difficult to accurately estimate high flows from stage height measurements.It is also possible that the single precipitation time series used in the simulations presented here did not adequately represent the amount of precipitation falling in the upper parts of the catchment.Monte Carlo analysis consistently produced better fitting parameter sets than those obtained through manual calibration.In some cases, the Monte Carlo sampler was able to improve on the calibrated parameter set in as few as 10 model runs, better performance than obtained through manual calibration was always observed within 100 jumps.This is somewhat concerning as it suggests that even expert manual calibration may fail to identify best-performing regions of parameter space.
The existence of a large number of equally credible parameter sets clearly demonstrates the existence of equifinality due to parameter uncertainty.Model performance was clearly sensitive to parameter values in urban land cover types.The relatively low groundwater time constants are consistent with a more flashy response in urban catchments, as has been observed elsewhere (Oni et al., 2013), however the high degree day evapotranspiration rates are somewhat surprising as it might be expected that urban areas would have less vegetation cover and hence lower evapotranspiration rates.
The PERSiST application presented here uses a fairly simple structural model where precipitation is routed vertically through three buckets and to the river.This perceptual model of runoff generation has been widely used for simulating surface water nutrient dynamics in the Thames (Jin et al., 2012;Crossman et al., 2013a;Whitehead et al., 2013) and elsewhere (Whitehead et al., 1998;Wade et al., 2002;Futter et al., 2009).
Performing a similar analysis to that of Fenicia et al. (2013) and evaluating alternate model structures, especially in the upper reaches, might improve model performance.
There are a number of different schools of thought about Monte Carlo analysis, with arguments made for both formal (Vrugt et al., 2009) (Beven, 2006).The approach presented here takes an informal approach that may appear overly simplistic but which should be robust to high dimensionality and potentially non-smooth goodness of fit response surfaces.It appeared that a chain length of 2500 was sufficient for identifying local optima in parameter space.Typically, the parameter set with the highest sum of NS statistics was identified within 1500 model runs and it was extremely rare to obtain better model fits after 2000 or more runs.As the goal of each MCMC run was to identify a locally optimal parameter set, this chain length was deemed adequate.The approach presented here, with an ensemble of locally optimal parameter sets, recognizes the equifinality inherent in catchment-scale hydrological modeling.There is a potentially infinite number of parameter sets all capable of providing credible fits to the observed data.As noted by Beven (2006), this is a common outcome in hydrological modeling.Unlike Laplace's demon, our understanding of the runoff generation process will always be incomplete and multiple competing hypotheses, expressed as combinations of model structures, parameter sets and environmental data may have equal credibility.
There has been considerable discussion in the hydrological literature about the pros and cons of multi-site calibration (Lerat et al., 2010;Cao et al., 2006;Zhang et al., 2005;Gong et al., 2012).Often, modelers are forced to use single site calibrations out of necessity as monitoring agencies rarely collect data at internal points in a catchment.The Thames is a fortunate exception to this trend.While we present a model application based on data from 8 gauging station on the main stem of the river, the UK National River Flow Archive maintains data collected at 122 sites throughout the catchment http://www.ceh.ac.uk/data/nrfa/index.html, visited 3 May 2013).Model simulations consistently showed higher NS statistics at Teddington, the most downstream gauging site, than at Pinkhill in the headwaters of the catchment (Fig. 5).There are a number of possible reasons for this.The model structure used here may be inappropriate for simulating flows in the upper reaches of the Thames, or there may be problems with the input and calibration data.The single rainfall time series used here may be more representative of precipitation patterns lower in the Thames catchment Figures than in the upper reaches.As the PERSiST framework allows each sub catchment to be linked to a unique rainfall time series, future work could attempt to address this question.
The PERSiST modeling framework has been designed to be as simple as possible but no simpler.Both snowmelt and evapotranspiration are dependent on energy balances, which may be poorly represented by temperature index approximations.Despite this shortcoming, simple temperature index snowmelt models are widely used and accepted.Hock (2003) notes that, at a catchment scale, temperature index methods can outperform more data hungry energy balance snowmelt models.However, temperature index methods can have limited temporal resolution and spatial accuracy (Hock, 2003).The temperature index approach to estimating evapotranspiration is similar to the one used in HBV (Saelthun, 1996), which has been the subject of some criticism (Andréasson et al., 2004;Lawrence and Haddeland, 2011).However, the temperature and soil moisture dependent approach to estimating actual evapotranspiration as the difference between precipitation and runoff appears to be robust and has been successfully defended elsewhere (Crossman et al., 2013b).The PERSiST framework assumes that water is transferred instantaneously between buckets.The developers of the SUPERFLEX framework have shown that adding in lag time for water transfers between buckets can improve model performance (Kavetski and Fenicia, 2011;Fenicia et al., 2013).Unlike the INCA approach, which is semi-distributed but does not route water laterally within the terrestrial environment, PERSiST facilitates the simulation of upslope and riparian areas and fluxes of water between them.As PERSiST has been designed primarily for simulating surface water fluxes, it uses a fairly simplistic representation of groundwater and does not currently include the capacity for loss of input of deep groundwater from outside the catchment.
Flexible modeling frameworks are needed in hydrology.Refsgaard (2007) presents a deceptively simple flowchart for decision making in the modeling process.A modeler develops a conceptual model based on field observations and, depending on the correspondence between the conceptual model and available models, either selects an ex-Introduction

Conclusions References
Tables Figures

Back Close
Full isting model or develops a new one from scratch.Unfortunately, the development of an entirely new model is not often an option, and the model user is forced to use someone else's perceptual model.This is not always ideal as it can introduce a known but unavoidable model structural uncertainty.The approach presented here is similar to that advocated by the developers of the SUPERFLEX modeling framework (Kavetski and Fenicia, 2011;Fenicia et al., 2011) insofar as it helps to overcome this problem.One potential difference between the PERSiST and SUPERFLEX approaches is that most published SUPERFLEX applications use a single, lumped representation of a catchment whereas PERSiST has been designed from the outset to be semi-distributed.
As stated in the name, one goal of PERSiST is improved hydrological simulations for solute transport.It is widely accepted by the hydrological modeling community that tracer data can improve the fit and credibility of hydrological models (Tetzlaff and Soulsby, 2008).In some ways, the modeling of pollutant fate and transport can be conceptualized as the use of non-conservative tracers to improve hydrological understanding.Simulating pollutant transport can help to constrain hydrological model predictions, but more importantly, use of appropriate hydrological model structures can aid in understanding the mechanisms behind pollutant transport.It is hoped that the PERSiST model makes an -albeit small -contribution to achieving this goal.Introduction

Conclusions References
Tables Figures

Conclusions References
Tables Figures

Back Close
Full  Full  Full  Full  Full  Full  Full   2002).A watershed (level 1) is represented as one or more reach/subcatchments (level 2).Within each sub catchment, there are one or more landscape units (level 3).Each landscape unit is made up of one or more buckets through which water is routed (level 4).Introduction Full depth (right).The maximum depth of water in a bucket (1; b1) can be partitioned into freely draining 582 water (2) and water that may contribute to evapotranspiration but not drainage (3; b2).The rate of 583 evapotranspiration (4) is not constrained when there is freely draining water in the bucket.When the 584 water depth drops below the freely draining depth, evapotranspiration is limited as a power function of 585 water depth.586 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | presents a deceptively simple flowchart for decision making in the modeling process.A modeler develops a conceptual model based on field observations and, depending on the correspondence between the conceptual model and available models, either selects an existing model or develops a new one from scratch.Discussion Paper | Discussion Paper | Discussion Paper | step, semi-distributed model designed primarily for simulating water and solute fluxes in rivers and streams based on solving the equation R = P − E .Key features of the model include (i) a user-specified model structure suitable for simulating multiple perceptual models of catchment water stores and flow pathways; (ii) semi-distributed flow routing incorporating runoff production from multiple land cover types and an ability to simulate flows at multiple points in a river; (iii) capacity to simulate inundation and infiltration of riparian areas; (iv) simple temperature index snowmelt and evapotranspiration routines; Discussion Paper | Discussion Paper | Discussion Paper | The conceptual hydrological model embodied in PERSiST draws heavily on both HBV and INCA.PERSiST uses the temperature index representations of snow dynamics and evapotranspiration that are incorporated in HBV and the semi-distributed landscape representation from INCA.A key feature differentiating PERSiST from both HBV and INCA is that it has a much more flexible representation of terrestrial hydrology which gives greater flexibility in model structure and an ability to simulate a much wider range of terrestrial hydrologic conditions.At its core, PERSiST is a conceptual, bucket-type model.A watershed is represented as a series of subcatchments made up of one or more landscape units (Fig. 1).Each Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | interception losses and transpiration.The rate of evapotranspiration from a bucket (E (b)) is determined by multiplying the landscape-scale maximum possible rate (E (l )) by the relative evapotranspiration index (b 4 ).Note that the relative evapotranspiration indices should sum Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | also tended to be over-estimated.The model failed to capture flows in the early part of the record in the upper part of the catchment and in the winter of 2002 (Fig. 6a).At Teddington, in the lowest part of the catchment, there was no evidence of bias in reproducing high or low flows.The most sensitive parameters included evapotranspiration and groundwater time constants in urban chalk and urban bedrock land cover types.(Fig. 8a, b and Table , we present a new model for simulating fluxes of water and solutes through heterogeneous landscapes and a simple Monte Carlo tool for model calibration and sensitivity analysis.The PERSiST model was able to reproduce the observed runoff dynamics at 8 sites in the Thames River using time series of air temperature and precipitation with data on areas of different landscape unit types in the catchment.Model performance was better in the lower reaches than the uppermost reach of the Thames.The median NS statistic in the uppermost reach was 0.7, versus 0.82 in the lower most reach.Inspection of the plots of modeled and observed flow values(Figs.6b, 7b)  show that there was less bias in model fit to the lowermost reach.This suggests some data uncertainty in the model results presented here.In some years at the uppermost reach (Fig.6a), the model over-estimated maximum flows, while in other years they were underestimated.Base flows were almost always within the range of predicted flows.In general, the range of predicted flows were tighter in the lower than upper reaches Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 :Fig. 1 .
Figure 1: Conceptual representation of the landscape in PERSiST adapted fromWade et al, (2002).A 566 watershed (level 1) is represented as one or more reach/subcatchments (level 2).Within each sub 567 catchment, there are one or more landscape units (level 3).Each landscape unit is made up of one or 568 more buckets through which water is routed (level 4).569 570

Figure 2 :Fig. 2 .Figure 3 :
Figure 2: Simple landscape unit comprised of three buckets representing direct runoff, soil water and 573 groundwater.The arrow labels (a-f) identify different fluxes and can be linked to the flow partitioning 574 matrix shown in Table 2. (a) represents direct runoff to the river; (b) is infiltration from the direct runoff 575 to soil water bucket; (c) represents saturation excess return flow from the quick bucket to the soil 576 surface while (d) is flow from the soil water bucket to the river.Flux (e) is from the soil water to 577 groundwater bucket while flux (f) represents water flow from the groundwater bucket to the river.578

Fig. 3 .Figure 4 :Fig. 4 .
Fig.3.Generic bucket structure (left) and relative evapotranspiration rate as a function of water depth (right).The maximum depth of water in a bucket (1; b 1 ) can be partitioned into freely draining water (2) and water that may contribute to evapotranspiration but not drainage (3; b 2 ).The rate of evapotranspiration (4) is not constrained when there is freely draining water in the bucket.When the water depth drops below the freely draining depth, evapotranspiration is limited as a power function of water depth.

Figure 5 :Fig. 5 .
Figure 5: Nash Sutcliffe goodness of fit statistics obtained for each reach from ensemble of 400 best 594 performing model calibrations.595 NS by Reach

Figure 8a :
Figure8a: Cumulative distribution function of maximum, minimum and most sensitive rates of 612 evapotranspiration (l 5 ).Note that the "urban bedrock" landscape unit type has the maximum rates of 613 evapotranspiration for all types.614

Table 2 .
Example square matrix with fluxes identified in Fig.2(upper) and candidate values (lower).

Table 3 .
PERSiST Internal Time Series.
V (i , j ) Volume of water transferred from bucket i to bucket j

Table 4 .
Sub catchment and reach dimensions along with proportional cover of different landscape element types.

Table 5 .
Parameter ranges used during MCMC analysis.
• C Land Cover Type Degree Day ET Rate 0.05 0.13 12 mm • C −1 d −1