SPHY v2.0: Spatial Processes in HYdrology

Introduction Conclusions References


Introduction
The number and diversity of water-related challenges are large and are expected to increase in the future. Even today, the ideal condition of having the appropriate amount of good-quality water at the desired place and time, is most often not satisfied (Biswas and change will intensify food insecurity by water shortages (Wheeler and von Braun, 2013), and loss of access to drinking water (Rockström et al., 2012). Current and future water related challenges are location and time specific and can vary from impact of glacier dynamics (Immerzeel et al., 2011), economic and population growth (Droogers et al., 2012), floods or extended and more prolonged droughts (Dai, 2011), amongst 5 others.
In response to these challenges hydrologists and water resources specialists are developing modeling tools to analyze, understand and explore solutions to support decision makers and operational water managers. Despite difficulties to connect the scientific advances in hydrological modeling with the needs of decision makers and 10 water managers, progress has been made and there is no doubt that modeling tools are indispensable in what is called good "water governance" (Droogers and Bouma, 2014;Liu et al., 2008).
The strength of models is that they can provide output on ultimate high temporal and spatial resolutions, and for difficult to observe sub-processes (Bastiaanssen et al., 15 2007). The most important aspect of applying models, is in their use to explore different scenarios, expressing for example, possible effects of changes in population and climate on the water cycle (Droogers and Aerts, 2005). Such scenarios are often referred to as projections. Models are also applied at the operational level to explore interventions (management scenarios) to be used by water managers and policy makers. over the last decade (Smith et al., 2004(Smith et al., , 2012. DMIP, and similar initiatives (Moreda et al., 2006;Khakbaz et al., 2012), have analysed the capacity of models to reproduce observed streamflow. It is however clear that the need of decision makers and water managers goes far beyond streamflow only. The notion that not only water in streams is relevant, but a better understanding of the full hydrological cycle is required, has led 10 to the popular division of "blue" and "green" water (Falkenmark and Rockström, 2010). Moreover, there is a clear need that not only rainfall-runoff processes are relevant, but an integrated approach encompassing processes as glacier and snowmelt, evapotranspiration, reservoirs, changes in landcover and landuse should be better integrated in our models (Bell et al., 2014;Liu et al., 2008). Finally, the linkage with satellite remote 15 sensing data has been shown to be very relevant to increase accuracy and usefulness of hydrological models, especially in data scarce areas (Droogers and Kite, 2002;Ines et al., 2006;Kite and Droogers, 2000;Wanders et al., 2014).
Based on the discussions above, there is a clear need for a hydrological model that combines the strengths of existing modelling approaches such that: (i) most relevant

Introduction
The SPHY model has been developed by combining the best components of existing 5 and well tested simulation models: HydroS (Droogers and Immerzeel, 2010), SWAT (Neitsch et al., 2009), PCR-GLOBWB (Sperna Weiland et al., 2012), SWAP (van Dam et al., 1997) and HimSim (Immerzeel et al., 2011). SPHY was developed with the explicit aim to simulate terrestrial hydrology at flexible scales, under various land use and climate conditions. SPHY is a spatially distributed leaky bucket type of model, and is ap- 10 plied on a cell-by-cell basis. The main terrestrial hydrological processes are described in a physically consistent way so that changes in storages and fluxes can be assessed adequately over time and space. SPHY is written in the Python programming language using the PCRaster (Karssenberg et al., 2001;Karssenberg, 2002;Karssenberg et al., 2010;Schmitz et al., 2009Schmitz et al., , 2013 dynamic modelling framework. 15 SPHY is grid-based and cell values represent averages over a cell (Fig. 1). For glaciers, sub-grid variability is taken into account: a cell can be glacier-free, partially glacierized, or completely covered by glaciers. The cell fraction not covered by glaciers consists of either land covered with snow or land that is free of snow. Land that is free of snow can consist of vegetation, bare soil, or open water. The dynamic vegetation 20 module accounts for a time-varying fractional vegetation coverage, which affects processes such as interception, effective precipitation, and potential evapotranspiration. Figure 2 provides a schematic overview of the SPHY model concepts.
The land compartment is divided in two upper soil stores and a third groundwater store, with their corresponding drainage components: surface runoff, lateral flow Introduction depending on the temperature. Any precipitation that falls on land surface can be intercepted by vegetation and in part or in whole be evaporated. The snow storage is updated with snow accumulation and/or snow melt. A part of the liquid precipitation is transformed in surface runoff, whereas the remainder infiltrates into the soil. The resulting soil moisture is subject to evapotranspiration, depending on the soil properties 5 and fractional vegetation cover, while the remainder contributes to river discharge by means of lateral flow from the first soil layer, and baseflow from the groundwater layer.
Melting of glacier ice contributes to the river discharge by means of a slow and fast component, being (i) percolation to the groundwater layer that eventually becomes baseflow, and (ii) direct runoff. The cell-specific runoff, which becomes available for 10 routing, is the sum of surface runoff, lateral flow, baseflow, snow melt and glacier melt.
If no lakes are present, then the user can choose for a simple flow accumulation routing scheme: for each cell the accumulated amount of water that flows out of the cell into its neighboring downstream cell is calculated. This accumulated amount is the amount of water in the cell itself plus the amount of water in upstream cells of the 15 cell. For each cell, the following procedure is performed: using the local drain direction network the catchment of a cell is determined which is made up the cell itself and all cells that drain to the cell. If lakes are present, then the fractional accumulation flux routing scheme is used: depending on the actual lake storage, a fraction of that storage becomes available for routing and is extracted from the lake, while the remaining part Introduction the Leaf Area Index (LAI) in order to estimate the growth-stage of land cover. For setting-up the model, data on streamflows are not necessary. However, to undertake a proper calibration and validation procedure flow data are required. The model could also be calibrated using actual evapotranspiration, soil moisture contents, and/or snow coverage.

5
The SPHY model provides a wealth of output data that can be selected based on the preference of the user. Spatial output can be presented as maps of all the hydrological processes. Maps that can be displayed as output include actual evapotranspiration, runoff generation (separated by its components), and groundwater recharge. These maps can be generated on daily base, but most users prefer to get those at monthly or 10 annual aggregated time periods. Time-series can be generated for each location in the study area. Time-series often used are stream flow under current and future conditions, actual evapotranspiration and recharge to the groundwater.

Modules
SPHY enables the user to turn on/off modules that are not required. This concept is 15 very useful if the user is studying hydrological processes in regions where not all hydrological processes are relevant. A user may for example be interested in studying irrigation water requirements in central Africa. For this region glacier and snow melting processes are less relevant, and can thus be switched off. Another user may only be interested in simulating moisture conditions in the first soil layer, allowing the possibility 20 to switch off the routing and groundwater modules. The advantages of turning off irrelevant modules are two-fold: (i) decrease model run-time, and (ii) decrease the amount of required model input data. Figure 3 represents an overview of the six modules available: glaciers, snow, groundwater, dynamic vegetation, simple routing, and lake/reservoir routing. All modules can 25 run independently from each other, except for the glaciers module. If glaciers are present, then snow processes are relevant as well (Verbunt et al., 2003;Singh and Kumar, 1997) module is turned on. Since melting glacier water percolates to the groundwater layer, the glaciers module cannot run with the groundwater module turned off. For routing two modules are available, being (i) a simple flow accumulation routing scheme, and (ii) a fractional flow accumulation routing scheme used when lakes/reservoirs are present. The user has the option to turn off routing, or to choose between one of these two rout-5 ing modules. All hydrological processes incorporated in the SPHY model are described in detail in the following sections.

Reference and potential evapotranspiration
Various methods exist to calculate the reference evapotranspiration (ET r ). From the many existing methods, the FAO-56 application of the Penman-Monteith equation 10 (Allen et al., 1998) has been most widely used, and can be considered as a sort of standard (Walter et al., 2001). Despite the good physical underlying theory of this method, its major drawback is its relatively high data demand. The Penman-Monteith equation requires air temperature, windspeed, relative humidity, and solar radiation data. In many regions around the world (e.g. Africa), these data are not readily measured 15 at meteorological stations, and therefore not suitable. Another well-known method to calculate the ET r is the Hargreaves method (Hargreaves and Samani, 1985). Droogers and Allen (2002) showed that in many regions around the globe, the Hargreaves and Penman-Monteith method showed comparable results. A lack of meteorological data brought (Hargreaves and Samani, 1985)  daily maximum and minimum air temperature. The constant 0.408 is required to convert the radiation to mm evaporation equivalents, and Ra can be obtained from tables (Allen et al., 1998) or equations using the day of the year and the latitude of the area of interest. According to Allen et al. (1998), the ET r is the evapotranspiration rate from a ref- 5 erence surface, not short of water. The reference surface is a hypothetical grass reference crop with specific characteristics. The potential evapotranspiration ET p has no limitations on crop growth or evapotranspiration from soil water and salinity stress, crop density, pests and diseases, weed infestation or low fertility. Allen et al. (1998)  transpiration and soil evaporation are integrated into the Kc. If the dynamic vegetation module in SPHY is not used, then the user can opt to use a single constant Kc throughout the entire simulation period or use a time-series of crop coefficients as model input. However, vegetation is generally very dynamic throughout the year, and therefore it is not very realistic to use a single constant Kc throughout 20 the entire simulation period. It is therefore more realistic to use a time-series of crop coefficients as model input, or use the dynamic vegetation module instead. This can be adjusted according to the user's preferences.
The Kc can be estimated using remotely sensed data (Rafn et al., 2008). In the dynamic vegetation module the Kc is scaled throughout the year using the NDVI and Introduction can easily be obtained from Allen et al. (1998) as input to improve model accuracy.

Canopy storage
SPHY allows the user to use the dynamic vegetation module in order to incorporate a time-variable vegetation cover and corresponding rainfall interception. In order to 10 calculate the rainfall interception, the canopy storage needs to be calculated, using a time-series of the Normalized Difference Vegetation Index (NDVI) (Carlson and Ripley, 1997). The first step involves the calculation of the fraction photosynthetically active radiation (FPAR). The FPAR can be calculated using a relation between the NDVI and the FPAR, which was found by Peng et al. (2012) and described by Sellers et al. (1996), 15 according to: with: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 0.001 is equivalent to a minimum LAI. In order to calculate the FPAR, an NDVI timeseries is required. The second step is the calculation of the Leaf-Area-Index (LAI), which is eventually required to calculate the maximum canopy storage (Scan max ). According to Monteith (1973), the LAI for vegetation that is evenly distributed over a surface can be calculated 5 using a logarithmic relation between the LAI and FPAR, according to:  Goward and Huemmrich (1992) is often used. However, since SPHY is generally applied using grid cell resolutions between 250 m and 1 km, we can assume that the effect of having vegetation concentrated in clusters is neglectable. Therefore, the calculation of the LAI 15 in SPHY is done using the logarithmic relation of Monteith (1973) (Eq. 6).
The next step involves the calculation of the maximum canopy storage (Scan max [mm]). Many different relations between Scan max and the LAI can be found in literature depending on the vegetation type (de Jong and Jetten, 2010). The best results for crop canopies are shown by Kozak et al. (2007) and are archived by Von Hoyningen-Huene 20 (1981) who derived the following relation between Scan max and the LAI: Note that even with an LAI of zero, Scan max is still close to 1 mm. Introduction

Interception
Interception is calculated on a daily basis if the dynamic vegetation module is used, and consists of the daily precipitation plus the intercepted water remaining in the canopy storage from the previous day. First of all the canopy storage is updated with the amount of precipitation of the current day: with Scan t [mm] the canopy storage on day t, Scan t−1 [mm] the canopy storage on day t − 1, and P t [mm] the amount of precipitation on day t. The portion of precipitation that cannot be stored in the canopy storage is known as precipitation throughfall, or effective precipitation, according to: with Pe t [mm] the effective precipitation on day t, and Scan t [mm] the canopy storage on day t. This equation shows that precipitation throughfall only occurs if the water stored in the canopy exceeds the maximum canopy storage. After the effective precipitation has been calculated, the canopy storage is updated as: The remaining amount of water stored in the canopy is available for interception, and the amount of water that will be intercepted depends on the atmospheric demand for open water evaporation. A commonly used value for the atmospheric demand for open water evaporation is 1.5 (Allen et al., 1998), which is derived from the ratio between 20 one and the mean pan evaporation coefficient Kp (∼ 0.65). If desired this value can easily be altered in the model code. The interception can now be calculated using: with Int t [mm] the intercepted water on day t, and ET r,t [mm] the reference evapotranspiration on day t. Finally, the canopy storage is updated by subtracting the interception: 2.5 Snow processes 5 For each cell a dynamic snow storage is simulated at a daily time step, adopted from on the model presented by Kokkonen et al. (2006). The model keeps track of a snow storage, which is fed by precipitation and generates runoff from snow melt. Refreezing of snow melt and rainfall within the snowpack are simulated as well. 10 Depending on a temperature threshold, precipitation is defined to fall in either solid or liquid state. Daily snow accumulation, which is defined as solid precipitation, is calculated as: for precipitation to fall as snow. The precipitation that falls as rain is defined as liquid precipitation, and is calculated as:

Snow and rainfall
with P l,t [mm] being the amount of rainfall on day t.

Snow melt, refreezing, and storage
To simulate snow melt, the well-established and widely used degree day melt modeling approach is used (Hock, 2003). The application of degree-day models is widespread in cryospheric models and is based on an empirical relationship between melt and air temperature. Degree-day models are easier to set up compared to energy-balance 5 models, and only require air temperature, which is mostly available and relatively easy to interpolate (Hock, 2005). Using a degree-day modeling approach, the daily potential snow melt is calculated as follows: with A pot,t [mm] the potential snow melt on day t, and DDF s [mm • C −1 d −1 ] a calibrated 10 degree day factor for snow. The actual snow melt is limited by the snow store at the end of the previous day, and is calculated as: with A act,t [mm] the actual snow melt on day t, and SS t−1 [mm] the snow store on day [t − 1]. The snow store from day [t − 1] is then updated to the current day t, using the 15 actual snow melt (A act,t ) and the solid precipitation (P s,t ). Part of the actual snow melt refreezes within the snow pack and thus does not runoff immediately. When temperature is below the melting point, melt water that has refrozen in the snow pack during [t − 1] is added to the snow store as: The capacity of the snow pack to refreeze snow melt is characterized by introducing a calibrated water storage capacity (SSC [mm mm −1 ]), which is the total mm water equivalent of snow melt that can refreeze per mm water equivalent of snow in the snow store. The maximum of melt water that can refreeze (SSW max [mm]) is thus limited by the thickness of the snow store: Then the amount of melt water stored in the snowpack, and that can refreeze in the next time-step, is calculated as: with SSW t the amount of melt water in the snow pack on day t, SSW max,t the maximum 10 of melt water that can refreeze on day t, SSW t−1 the amount of refrozen melt water on day [t − 1], P l,t the amount of rainfall on day t, and A act,t the actual snow melt on day t. The units of all terms are in mm. The total snow storage (SST [mm]) consists of the snow store and the melt water that has refrozen within it, according to: with (1−GlacF) [-] the gridcell fraction not covered with glaciers. In SPHY it is therefore assumed that snow accumulation and melt can only occur on the gridcell fraction determined as land surface. Snow falling on glaciers is incorporated in the glacier module.

Snow runoff
Runoff from snow (SRo [mm]) is generated when the air temperature is above melting point and no more melt water can be refrozen within the snow pack, according to: with ∆SSW [mm] the change in melt water stored in the snowpack according to:

Glacier processes
Since the SPHY model usually operates at a spatial resolution between 250 m and 1 km, the dynamics of glaciers such as ice flow cannot be resolved explicitly. Therefore, glaciers in SPHY are considered as melting surfaces which can completely or partly 10 cover a grid cell.

Glacier melt
Glacier melt is calculated with a degree day modeling approach as well (Hock, 2005). Because glaciers covered with debris melt at different rates than debris-free glaciers (Reid et al., 2012), a distinction can be made between different degree day factors for 15 both types. The daily melt from debris free glaciers (A CI [mm]) is calculated as: is a degree day factor for debris covered glaciers and F DC [-] is the fraction of debris covered glaciers within the fractional glacier cover of a grid 5 cell. The total glacier melt per grid cell (A GLAC [mm]) is then calculated by summing the melt from the debris-covered and debris-free glacier types and multiplying by the fractional glacier cover, according to: 2.6.2 Glacier runoff 10 In SPHY a fraction of the glacier melt percolates to the ground water while the remaining fraction runs off. The distribution of both is defined by a calibrated glacier melt runoff factor (GlacROF [-]) which can have any value ranging from 0 to 1. Thus the generated runoff GRo [mm] from glacier melt is defined as:

Glacier percolation
The percolation from glacier melt to the groundwater (G perc [mm]) is defined as: The glacier water that percolates to the groundwater is added to the water percolating from the soil layers of the non-glacierized part of the grid cell (Sects. 2.7.1 and 2.7.6) 20 which will then recharge the groundwater.

Soil water balances
The soil water processes in SPHY are modelled for three soil layers (Fig. 2), being (i) the first soil layer (rootzone), (ii) second soil layer (subzone), and (iii) third soil layer (groundwater layer). The water balance of the first soil layer is: with SW 1,t and SW 1,t−1 the water content in the first soil layer on day t and [t − 1], respectively, Pe t the effective precipitation on day t, Ea t the actual evapotranspiration on day t, RO t the surface runoff on day t, LF 1,t the lateral flow from the first soil layer on day t, Perc 1,t percolation from the first to the second soil layer on day t, Cap t capillary 10 rise from the second to the first soil layer on day t. The second soil layer water balance is: with SW 2,t and SW 2,t−1 the water content in the second soil layer on day t and [t − 1], respectively, and Perc 2,t percolation from the second to the third soil layer on day t. 15 The third soil layer water balance is given as: with SW 3,t and SW 3,t−1 the water content in the third soil layer on day t and [t − 1], respectively, Gchrg t groundwater recharge from the second to the third soil layer on day t, and BF t baseflow on day t. If the glacier module is used, then groundwater 20 recharge consists of percolation from the second soil layer and percolated glacier melt, otherwise only percolation from the second soil layer is taken into account. The user can opt to run SPHY without the third soil layer (groundwater). This may be desirable if the user for example is mainly interested in simulating soil moisture 1704 Introduction conditions in the rootzone, instead of evaluating e.g. the contribution of baseflow to the total routed river flow. In that case only the two upper soil layers are used where the bottom boundary of soil layer two is controlled by a seepage flux (pos. outward), and instead of baseflow from the third soil layer, water leaves the second soil layer through lateral flow. With the groundwater module turned off, the water balance for the second 5 soil layer is: with LF 2,t lateral flow from the second soil layer, and Seep seepage in or out of the second soil layer (pos. is outgoing). The units for all water balance terms are mm. 10 Evapotranspiration refers to both the transpiration from vegetation and the evaporation from soil or open water. As was mentiond in Sect. 2.3, the Kc accounts for both the crop transpiration and soil evaporation. The additional use of the dynamic vegetation module accounts for a time-variable vegetation cover, meaning that the role of evaporation becomes more dominant as soons as vegetation cover decreases. 15 Many limiting factors (e.g. salinity stress, water shortage, water excess, diseases) can cause a reduction in potential evpotranspiration (ET p ), resulting in the actual evapotranspiration rate (ET a ). Since SPHY is a water-balance model, SPHY only accounts for stresses related to water shortage or water excess. If there is too much water in the soil profile, then the plant is unable to extract water because of oxygen stress 20 (Bartholomeus et al., 2008). The calculation of evapotranspiration reduction due to water excess (oxygen stress) is quite complex and requires a vast amount of parameters that is generally not available for the spatial scale that SPHY is applied on. Therefore, SPHY uses a evapotranspiration reduction parameter (ETred wet ) that has a value of 0 if the soil is saturated, and otherwise it will have a value of 1. This parameter is used in 25 the following equation to calculate the actual evapotranspiration: with ET a,t [mm] the actual evapotranspiration on day t, ET p,t [mm] the potential evapotranspiration on day t, and ETred wet and ETred dry being the reduction parameters for water excess and water shortage conditions, respectively. The ETred dry is calculated using the Feddes equation (Feddes et al., 1978), which assumes a linear decline in rootwater uptake if the water pressure head drops below a critical value. This critical 5 value can be determined using the soil water retention curve (pF-curve), which relates the soil type to its water binding capacity. This binding capacity is a suction force, and is often expressed in cm negative water column. The pF is simply calculated as:

Surface runoff
Within the field of hydrology two types of surface runoff can be distinguished: Hewlettian (Hewlett, 1961) and Hortonian (Beven, 2004;Corradini et al., 1998)  runoff is often referred to as saturation excess overland flow, while Hortonian runoff is also known as infiltration excess overland flow. Hewlettian runoff occurs when the soil is saturated, and more rainfall is added, resulting in direct surface runoff. Hortonian runoff occurs when the rainfall intensity exceeds the infiltration capacity of the soil: during wetting of the soil, the soil infiltration capacity decreases over time, and if at a certain 5 moment the rainfall intensity exceeds the soil infiltration capacity, surface runoff occurs. Since the SPHY model runs on a daily time-step, the model does not account for sub-daily variability in rainfall intensities, meaning that Hortonian runoff processes can be considered as less important. For this reason SPHY uses the Hewlettian runoff concept to calculate surface runoff. Surface runoff is calculated from the first soil layer: with RO [mm] surface runoff, SW 1 [mm] the water content in the first soil layer, and SW 1,sat [mm] the saturated water content of the first soil layer.

Lateral flow
Lateral flow is substantial in catchments with steep gradients and soils with high hy-15 draulic conductivities (Beven, 1981;Beven and Germann, 1982;Sloan and Moore, 1984). In SPHY it is assumed that only the amount of water exceeding field capacity can be used for lateral flow. Therefore, first the drainable volume of water (excess water) needs to be calculated:  Sloan and Moore (1984), the lateral flow at the hillslope outlet can be calculated as: with LF * l [mm] lateral flow from soil layer l , W l ,excfrac [-] the drainable volume of water as fraction of the saturated volume, and v lat,l [mm d −1 ] the flow velocity at the outlet. In SPHY, the drainable volume as fraction of the saturated volume is calculated as: The slope (slp) in SPHY is calculated for each gridcell as the increase in elevation per 10 unit distance. According to Neitsch et al. (2009), only a fraction of lateral flow will reach the main channel on the day it is generated if the catchment of interest is large with a time of concentration is greater than 1 day. This concept is also implemented in the SPHY model, and uses a lateral flow travel time TT lag,l [d] to lag a portion of lateral flow 15 release to the channel: A larger lateral flow travel time will result in a smoother streamflow hydrograph.

5
If the groundwater module is used, then water can percolate from the first to the second soil layer, and from the second to the third soil layer. If the user decides to run SPHY without the groundwater module, percolation only occurs from the first to the second soil layer. In SPHY water can only percolate if the water content exceeds the field capacity of that layer, and the water content of the underlying layer is not saturated.

10
A similar approach has been used in the SWAT model (Neitsch et al., 2009). The water volume available for percolation to the underlying layer is calculated as: . This approach follows the storage routing methodology, which is also implemented in the SWAT model (Neitsch et al., 2009): with w l ,perc [mm] the amount of water percolating to the underlying soil layer. The travel time for percolation is calculated in the same way as the travel time for lateral flow (Eq. 41).

Groundwater recharge
Water that percolates from the second to the third soil layer will eventually reach the 5 shallow aquifer. This process is referred to as groundwater recharge hereafter. If the glacier module is used as well, then also glacier melt that percolates contributes to the groundwater recharge. Groundwater recharge often does not occur instantaneous, but with a time lag that depends on the depth of the groundwater table and soil characteristics. SPHY uses the same exponential decay weighting function as proposed 10 by Venetis (1969) and used by Sangrey et al. (1984) in a precipitation groundwater response model. This approach has also been adopted in the SWAT model (Neitsch et al., 2009), using: with Gchrg t [mm] and Gchrg t−1 [mm] the groundwater recharge on day t and t − 1, 15 respectively. δ gw [d] is the delay time and w 2,perc [mm] is the amount of water that percolates from the second to the third layer on day t.

Baseflow
After groundwater recharge has been calculated, SPHY calculates baseflow which is defined as the flow going from the shallow aquifer to the main channel. Baseflow only 20 occurs when the amount of water stored in the third soil layer exceeds a certain threshold (BF thresh ) that can be specified by the user. Baseflow calculation in SPHY is based on the steady-state response of groundwater flow to recharge (Hooghoudt, 1940) and the water table fluctuations that are a result of the non-steady response of groundwater flow to periodic groundwater recharge (Smedema and Rycroft, 1983 (Neitsch et al., 2009) assumes a linear relation between the variation in groundwater flow (baseflow) and the rate of change in water table height, according to: 10 with BF t [mm] the baseflow into the channel on day t, and BF t−1 [mm] the baseflow into the channel on day t − 1. Since this equation has proven its succes in the SWAT model (Neitsch et al., 2009) throughout many applications worldwide, this equation has been adopted in the SPHY model as well. The baseflow recession coefficient (α gw ) is an index that relates the baseflow re-15 sponse to changes in groundwater recharge. Lower values for α gw therefore correspond to areas that respond slowly to groundwater recharge, whereas higher values indicate areas that have a rapid response to groundwater recharge. The baseflow recession coefficient is generally used as a calibration parameter in the SPHY model, but a good first approximation of this coefficient can be calculated using the number of 20 baseflow days (Neitsch et al., 2009):

Streamflow routing
After calculating the different runoff components, the cell specific total runoff (QTot) is calculated by summarizing these different runoff components. Depending on the modules being switched on, the different runoff components are (i) rain runoff (RRo), (ii) snow runoff (SRo), (iii) glacier runoff (GRo), and (iv) baseflow (BF). If the groundwater 5 module is not used, then baseflow is calculated as being the lateral flow from the second soil layer. Rainfall runoff is the sum of surface runoff (RO, Sect. 2.7.3) and lateral flow from the first soil layer (LF 1 , Sect. 2.7.4). Then QTot is calculated according to: methods require a flow direction network map, which can be obtained by delineating a river network using PCRaster or GIS software in combination with a Digital Elevation Model (DEM).

Routing
In hydrology streamflow routing is referred to as being the transport of water through an 20 open channel network. Since open-channel flow is unsteady, streamflow routing often involves the solving of complex partial differential equations. The St. Venant equations (Brutsaert, 1971;Morris and Woolhiser, 1980) are often used for this, but these have high data requirements related to the river geometry and morphology, which is unavailable for the spatial scale SPHY is generally applied on. Additionally, solving these equations require the use of very small time steps, which result in large model calculation times. Other models, such as e.g. SWAT (Neitsch et al., 2009), use the Manning equation (Manning, 1989) to define the rate and velocity of river flow in combination with the variable storage (Williams, 1975) or Muskingum (Gill, 1978) routing methods to obtain river streamflow. But also the Manning equation requires river bed dimen-5 sions, which are generally unknown on the spatial scale that SPHY generally is applied on. Therefore, SPHY calculates for each cell the accumulated amount of water that flows out of the cell into its neighbouring downstream cell. This can easily be obtained by using the accuflux PCRaster builtin function. If only the accuflux function would be used, 10 then it is assumed that all the specific runoff generated within the catchment on one day will end up at the most downstream location within one day, which is not plausible. Therefore, SPHY implements a flow recession coefficient (kx [-]) that accounts for flow delay, which can be a result of channel friction. Using this coefficient, river flow in SPHY is calculated using the three equations shown below: snow melt to the total routed runoff. However, this increases model run-time substantially, because the accuflux function, which is a time consuming function, needs to be called multiple times depending on the number of flow contributors to be routed.

Lake/reservoir routing
Lakes or reservoirs act as a natural buffer, resulting in a delayed release of water from 5 these water bodies. SPHY allows the user to choose a more complex routing scheme if lakes/reservoirs are located in their basin of interest. The use of this more advanced routing scheme requires a known relation between lake outflow and lake level height (Q(h) relation) or lake storage. To use this routing scheme, SPHY requires a nominal map with cells identified as 10 lakes each have a unique ID, while the non-lake cells have a value of zero. The user can supply a boolean map with True for cells that have measured lake levels, and False for lake cells that do not have measured lake levels. This specific application of SPHY is discussed in detail in Sect. 3.3. Four different relations can be chosen to calculate the lake outflow from the lake 15 level height or lake storage, being: (i) an exponential relation, (ii) a first-order polynomial function, (iii) a second-order polynomial function, and (iv) a third-order polynomial function. The user needs to supply maps containing the coefficients used in the different functions. The lake/reservoir routing scheme simply keeps track of the actual lake storage, 20 meaning that an initial lake storage should be supplied. Instead of the simple accuflux function described in the previous section, the lake/reservoir routing scheme uses the PCRaster functions accufractionstate and accufractionflux. The accufractionflux calculates for each cell the amount of water that is transported out of the cell, while the accufractionstate calculates the amount of water that remains stored in the cell. For 25 non-lake cells the fraction that is transported to the next cell is always equal to one, while the fraction that is transported out of a lake/reservoir cell depends on the actual lake storage. Each model time-step the lake storage is updated by inflow from 1714 Introduction upstream. Using this updated storage, the lake level and corresponding lake outflow can be calculated using one of the four relations mentioned before. The lake outflow can then be calculated as a fraction (Q frac [-]) of the actual lake storage. Instead of using Eq. (50), the Q frac is then used in Eqs. (52) and (53)  non-lake cells, the accufractionflux function becomes equal to the accuflux function used in the previous section. This actually means that for the river network the same routing function from Sect. 2.8.1 is used, and that Eqs. (52) and (53) only apply to lake/reservoir cells. In order to account for non-linearity and slower responding catchments, the same kx 15 coefficient is used again. This involves applying Eq. (51) as a last step after Eq. (52) and converting the units from m 3 d −1 to m 3 s −1 . Since the accufractionflux and accufraction state functions are more complex to compute, the use of these functions increases model run-time.

20
The SPHY model has been applied and tested in various studies ranging from realtime soil moisture predictions in lowlands, to operational reservoir inflow forecasting applications in mountainous catchments, to irrigation scenarios in the Nile Basin, to detailed climate change impact studies in the snow-glacier-rain dominated Himalayan region. Some typical applications will be summarized in the following sections.

Irrigation management in lowland areas
As SPHY produces spatial outputs for the soil moisture content in the root zone and the potential and actual evapotranspiration (ET), it is a useful tool for application in agricultural water management decision support. By facilitating easy integration of remote sensing data, crop growth stages can be spatially assessed at different moments in 5 time. The SPHY dynamic vegetation module ensures that all relevant soil water fluxes correspond with crop development stages throughout the growing season. Spatially distributed maps of root water content and ET deficit can be produced, enabling both the identification of locations where irrigation is required and a quantitative assessment of crop water stress.

10
SPHY has been applied with the purpose of providing field-specific irrigation advice for a large-scale farm in western Romania, comprising 380 individual fields and approximately ten different crops. Contrary to the other case studies highlighted in this paper, a high spatial resolution is very relevant for supporting decisions on variable-rate irrigation. The model has therefore been set up using a 30 m resolution, covering the 2013 15 and 2014 cropping seasons on a daily time step. Optical satellite data from Landsat 8 (USGS, 2013) were used as input to the dynamic vegetation module. Soil properties were derived from the Harmonized World Soil Database (Batjes et al., 2012), which for Romania contains data from the Soil Geographical Database for Europe (Lambert et al., 2003). Using the Van Genuchten equation (Van Genuchten, 1980), soil saturated 20 water content, field capacity, and wilting point were determined for the HWSD classes occurring at the study site. Elevation data was obtained from the EU-DEM dataset (EEA, 2014), and air temperature was measured by two on-farm weather stations.
In irrigation management applications, a model should be capable of simulating the moisture stress experienced by the crop due to insufficient soil moisture contents, 25 which manifests itself by an evapotranspiration deficit (potential ET − actual ET > 0). Field measurements of soil moisture and/or actual ET are therefore desired for calibration purposes. In this case study, a capacitance soil moisture sensor was installed Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | on a soybean field to monitor rootzone water content shortly after 1 May 2014, which is the start of the soybean growing season. The sensor measures volumetric moisture content for every 10 cm of the soil profile up to a depth 60 cm. It is also equipped with a rain gauge measuring the sum of rainfall and applied irrigation water, which was used as an input to SPHY. Soil moisture measured over the extent covered by the crop root 5 depth was averaged and compared to simulated values (Fig. 4).
As can be seen in Fig. 4, the temporal patterns as measured by the soil moisture sensor are well simulated by the SPHY model. Based on daily soil moisture values, a Nash-Sutcliffe (Nash and Sutcliffe, 1970) model efficiency coefficient of 0.6 was found, indicating that the quality of prediction of the uncalibrated SPHY model is "good" (Foglia et al., 2009).
The performance of the SPHY model can be improved by calibration using remotely sensed evapotranspiration (Immerzeel and Droogers, 2008), although such data is often not available on these small scales as ET is a very complex parameter to assess (Samain et al., 2012). It should also be noted that soil moisture content is typically 15 highly variable in space; a very high correlation between point measurements and grid cell simulations of soil moisture may therefore not always be feasible (Bramer et al., 2013).

Snow and glacier-fed river basins
SPHY is being used in large Asian river basins with significant contribution of glacier- 20 and snow melt to the total flow (Lutz et al., 2012b. The major goals of these applications are two-fold: -Assess the current hydrological regimes at high-resolution; e.g. assess spatial differences in the contributions of glacier melt, snow melt and rainfall-runoff to the total flow Rivers originating in the high mountains of Asia are considered to be the most meltwater dependent river systems on Earth (Schaner et al., 2012). In the regions surrounding the Himalayas and the Tibetan Plateau large human populations depend on the water supplied by these rivers (Immerzeel et al., 2010). However, the dependency on melt water differs strongly between river basins as a result of differences in climate 5 and differences in basin hypsometry (Immerzeel and Bierkens, 2012). Only by using a distributed hydrological modelling approach that includes the simulation of key hydrological and cryospheric processes, and inclusion of transient changes in climate, snow cover, glaciers and runoff, appropriate adaptation and mitigation options can be developed for this region (Sorg et al., 2012). The SPHY model is very suitable for such 10 an approach, and has therefore been widely applied in the region.
For application in this region, SPHY was setup at a 1 km spatial resolution using a daily time-step, and forced with historical air temperature (Tavg, Tmax, Tmin) and precipitation data, obtained from global and regional datasets (e.g. APHRODITE (Yatagai et al., 2012), Princeton (Sheffield et al., 2006), TRMM (Gopalan et al., 2010)) or 15 interpolated WMO station data from a historical reference period. For this historical reference period SPHY was calibrated and validated using observed streamflow. For the future period, SPHY was forced with downscaled climate change projections obtained from General Circulation Models (GCMs), as available through the Climate Model Intercomparison Projects (e.g. CMIP3, Meehl et al., 2007, CMIP5, Taylor et al., 2012, 20 and which were used as basis for the Assessment Reports prepared by the Intergovernmental Panel on Climate Change (IPCC).
In Central Asia, SPHY was succesfully applied in a study (ADB, 2012;Lutz et al., 2012b, a) that focused on the impacts of climate change on water resources in the Amu Darya and Syr Darya river basins. SPHY was used to quantify the hydrological 25 regimes in both basins, and subsequently to project the outflow from the upstream basins to the downstream areas by forcing the model with an ensemble of 5 CMIP3 GCMs. The SPHY model output fed into a water allocation model that was setup for the downstream parts of the Amu Darya and Syr Darya river basins. In the Himalayan Climate Change Adaptation Programme (HICAP), led by the International Centre for Integrated Mountain Development (ICIMOD), SPHY has been successfully applied in the upstream basins of the Indus, Ganges, Brahmaputra, Salween and Mekong rivers (Lutz et al., 2013. In this study the hydrological regimes of these five basins have been quantified and the calibrated and validated model (Fig. 5) 5 was forced with an ensemble of eight GCMs to create water availability scenarios until 2050. SPHY allowed the assessment of current and future contribution of glacier melt and snow melt to total flow (Fig. 6), and how total flow volumes and the intra-annual distribution of river flow will change in the future. 10 In harsh environments like Chile, hydropower companies have been using forecasting solutions such as statistical regressions of flow observations, or mathematical simulations with a minimum content of real-time physical information. These solutions are generally based on long historical data sets of flow observations, and do not take into account the current hydrological state of the basin and the climate for the coming 15 weeks/months. During the INTOGENER project, the SPHY model was used to demonstrate the operational capabilities of a water flow monitoring and prediction system aimed at hydropower production and water management organizations.

Flow forecasting
The INTOGNER system uses measurements derived from both Earth Observation (EO) satellites and in-situ sensors. Data that were retrieved from EO satellites con-20 sisted of a DEM map and a map-series of snow cover. Snow cover images were retrieved on a weekly basis, using RADARSAT and MODIS (Parajka and Blöschl, 2008;Hall et al., 2002) imagery. These images were used to update the snow store (SS [mm]) in the model in order to initialize the model for the forecasting period. Discharge, precipitation and temperature data were collected using in-situ meteorological stations. 25 In order to calculate the lake outflow accurately, the SPHY model was initialized with water level measurements retrieved from reflected Global Navigation Satellite System (GNSS) signals in Laja Lake. Static data that was used in the SPHY model consisted 1719 Introduction  jes et al., 2009) and land use data obtained from the GLOBCOVER (Bontemps et al., 2011) product. The SPHY model was setup to run at a spatial resolution of 200 m. Figure 8 shows the observed vs. simulated daily streamflow for two locations within the Laja River Basin for the historical period 2007-2008. It can be seen that model per-5 formance is quite satisfactory for both locations, with volume errors of −4 and −9.4 % for Canal Abanico (downstream of Laja Lake) and Rio Laja en Tucapel, respectively. The NS-coefficient (Nash and Sutcliffe, 1970), which is especially useful to assess the simulation of high discharge peaks, is less satisfactory for these locations. Hydropower companies, however, have more interest in expected flow volumes for the coming weeks/months than for accurate day-to-day flow simulations, and therefore the NS-coefficient is less important in this case. If the NS-coefficient is calculated for the same period on a monthly basis, then the NS-coefficients are 0.53 for Canal Abanico and 0.81 for Rio Laja en Tucapel. It is likely that SPHY model performance would even be better if a full model calibration would have been performed. 15 The hydropower company's main interest is the model's capacity to predict the total expected flow for the coming weeks during the melting season (October 2013 through March 2014). Figure 9 shows the bias between the total cumulative forecasted flow and observed flow for the 23 model runs that were executed during operational mode. Although there are some bias fluctuations in the Rio Laja en Tucapel model runs, it can 20 be concluded that the bias decreases for each next model run for both locations, which is a logical result of a decreasing climate forcing uncertainty as the model progresses in time. It can be seen that the SPHY model streamflow forecasts for Canal Abanico, which is downstream of Laja Lake, are substantially better than for Rio Laja en Tucapel (most downstream location). The reason for this has not been investigated during the 25 demonstration study, but since model performance for these two locations was satisfactory during calibration, a plausible explanation could be the larger climate forecast uncertainty in the higher altitude areas (Hijmans et al., 2005;Rollenbeck and Bendix, 2011) in the northeastern part of the basin that contributes to the streamflow of Rio Laja en Tucapel. Additionally, only two in-situ meteorological stations were available during operational mode, whereas during calibration 20+ meteorological stations were available. Moreover, these operational meteorological stations were not installed on higher altitudes, where precipitation patterns tend to be spatially very variable (Wagner et al., 2012;Rollenbeck and Bendix, 2011).

4 Future outlook
Further development and refinement of the SPHY model are foreseen in six areas, being (i) evapotranspiration processes, (ii) depression storage, (iii) representation of lakes/reservoirs, (iv) streamflow routing, (v) cryospheric processes, and (vi) the development of a graphical user-interface. 10 The current version of the SPHY model calculates the reference evapotranspiration using the modified Hargreaves equation (Hargreaves and Samani, 1985), which requires three meteorological variables: the average, maximum and minimum daily air temperature. Although this method is less data demanding than the well-known Penman-Monteith (Allen et al., 1998) equation, additional methods for the calculation 15 of the reference evapotranspiration that are even less data demanding should be explored for future versions of the SPHY model.
Besides the dynamic simulation of a time-varying fractional vegetation coverage, using NDVI time-series as input, it is currently assumed that the root depth of a growing crop is constant throughout the simulation period, and that it is equal to the depth of the 20 first soil layer. In reality a crop is seeded, starts to develop its root system, and finally gets harvested. Although this root development process is less dynamically for forests and therefore more relevant for agricultural crops, it can be seen as an improvement to be included in the SPHY model. A similar approach as in the SWAP model (van Dam et al., 1997)  and that SPHY is generally applied on larger spatial scales where cropping calendars, and thus corresponding root depth information is harder to obtain. Surface runoff in SPHY occurs whenever excess rainfall results in saturation of the first soil layer. Subsequently, this excess rainfall leaves the grid-cell as surface runoff without any obstacles or friction, and enters the river system within the same day. In 5 reality, surface runoff can be delayed because of surface friction, or stored in local depressions or man-made ponds (surface irrigation), before it is released to the river system. The concept of having a ponding layer is mainly beneficial for (i) agricultural applications where surface irrigation in combination with a ponding layer plays a major role, e.g. rice irrigation, and/or (ii) if the model is run on very high spatial resolution where the role of local depressions becomes more prominent; water can be stored for hours/days in these depressions, and water may be partly or completely evaporated before it has the chance to flow to the channels. The effect of delayed surface runoff can be related to dense vegetation or agricultural land management practices (Hunink et al., 2013;Kauffman et al., 2014), with bench terraces and contour tillage being ex- 15 amples of these practices. The effect of delayed or stored surface runoff is currently not implemented in the SPHY model, but should be further explored for future SPHY versions.
As was mentioned in Sect. 2.8.2, each individual lake is represented by maximum one grid-cell in SPHY. This is related to the PCRaster flow direction calculation: each 20 cell needs to flow to a downstream cell, and the total flow in the most downstream cell should be equal to the accumulated flow of all upstream cells plus the flow generated within that cell. This means that large lakes, consisting of multiple neighboring cells, can be seen as an enormous pit in the flow direction network, meaning that streamflow routing from the upstream lake cells is interrupted at these pit cells. Since large lakes 25 can evaporate a substantial amount of water (Gat et al., 1994), a lake-cell merging procedure should be developed for future versions of the SPHY model, that allows to automatically re-create the flow direction network based on the locations of the lake inflow cells and outflow cell, such that all the accumulated water of the inflow cells is allocated and stored in the downstream cell, where it can continue its way downstream. Although the currently implemented routing scheme in SPHY has proven its success during various applications (see Sects. 3.2 and 3.3), improvements in the routing scheme are foreseen for future SPHY versions. Using the current approach it is as- 5 sumed that the open channel surface area equals the grid cell area. If coarse model resolutions are used, then this assumption becomes less plausible, because in reality the surface area of a river would only cover a small fraction of this grid cell area. Therefore, it would be interesting to explore using more advanced routing schemes that take into account the river's cross section dimensions and corresponding velocity 10 flow rates. An example of a more advanced routing scheme, that takes into account the channel dimensions, and has been found to yield good results under a wide range of conditions in natural rivers (Brutsaert, 2005) is the Muskingum method (Gill, 1978). Considering the larger spatial resolutions that SPHY is generally applied for, it seems almost infeasible to determine the channel dimensions. However, some emperical re- 15 lations between the bank-full dischargen and channel dimensions are available (Park, 1977), and should be explored for implementing in SPHY.
Sublimation of snow is an important component of the water balance in areas experiencing snow cover during significant time of the year (Strasser et al., 2008). Different approaches have been developed to estimate this flux in hydrological models (Bowl-20 ing et al., 2004;MacDonald et al., 2009;Lenaerts et al., 2010;Groot Zwaaftink et al., 2013). Development of a parameterization of sublimation using the low data requirement of SPHY is foreseen.
In SPHY, glaciers are considered as entities generating melt using a temperatureindex model. Improvements in the representation of glaciers are foreseen to make them 25 mass-conserving, e.g. considering the precipitation falling in their accumulation areas. In the current SPHY version, all surface runoff that is generated on glacier covered areas is considered as glacier melt, and inclusion of further specification to seasonal snow melt and glacier ice melt is foreseen. Besides, the temperature-index model can be improved by including the incoming radiation in addition to air temperature in the temperature-index model (Hock, 2003;Pellicciotti et al., 2005;Heynen et al., 2013). Additionally, the glacier module will be extended with a parameterization for modeling glacier dynamics, which would enable quantifying the retreat or advance of glaciers as a result of climate perturbations. 5 When SPHY is run at a spatial resolution of 1 km or coarser for mountainous regions, improvement of the respresentation of sub-grid processes becomes useful. For example, when the daily average air temperature for a given 1 km grid cell is 0 • C, no melt will be simulated for that grid cell, although the subgrid variability in elevation shows that part of the grid cell is lower than the average elevation of the grid cell and thus has 10 temperatures above 0 • C, resulting in the generation of melt water. Although the model is relatively easy to understand and applicable by hydrologists and scientists with basic skills, it would even be better if the model could be used by a wider group of users with basic hydrological and computer skills. Therefore, the objective of an ongoing project is to develop a Graphical User Interface (GUI) for the 15 SPHY model. This GUI will be developed as plugin for the open-source QGIS environment, making it easy setup the model and analyse model input and output spatially and temporally.

Conclusions
In response to the diversity of water-related challenges, hydrologists have been devel-20 oping tools to analyze, understand and explore solutions to support decision makers and operational water managers. Although progress has been made, hydrologists have put a strong emphasis on streamflow analysis and forecasting, while ignoring other hydrological processes. It is clear that not only rainfall-runoff processes are relevant, but that an integrated approach where processes as glacier and snowmelt, evapotranspira- The objective of this paper is to introduce and present the SPHY model, its development background, and demonstrate some typical applications. Compared to other hydrological models, that typically focus on the simulation of streamflow only, the SPHY model has several advantages: it (i) integrates most relevant hydrological processes, (ii) is setup modular, (iii) is easy adjustable and applicable, (iii) can easily be linked 5 to remote sensing data, and (iv) can be applied for operational as well as strategic decision support.
The most relevant hydrological processes that are integrated in the SPHY model are rainfall-runoff processes, cryosphere processes, evapotranspiration processes, the simulation of dynamic vegetational cover, lake/reservoir outflow, and the simulation of 10 rootzone moisture contents. The capability of SPHY to successfully simulate rainfallrunoff and cryosphere processes was proven during its applications in the snow and glacier-fed river basins in Asia, and in Chile, where it was used to forecast streamflow during the snow melting season. Both the applications in Chile and in Romania show the strength of SPHY to include remotely sensed data as dynamic model input: in Chile 15 remotely sensed snow cover was used to implement a time-variable fractional snow coverage, and in Romania the NDVI was used to simulate dynamic vegetational cover. Whereas the application in the glacier-fed river basins in Asia to study the effects of climate change on the available water resources is a typical example of the use of SPHY to support strategic decision makers, the application in Romania for real-time irrigation 20 support, and the application in Chile to support hydropower companies for their reservoir management, are typical examples of the use of SPHY for operational purposes. The application of real-time irrigation support in Romania emphasizes the strength of SPHY to succesfully simulate other hydrological processes than rainfall-runoff as well, being evapotranspiration and rootzone soil moisture simulations. The different spatial 25 resolutions and orograhically and climatologically regions in which SPHY was applied, demonstrate its flexibility in scaling and applicability: in Romania it was applied at the farm-field level, requiring a spatial resolution of 30 m, whereas in Chile and Asia it was applied at resolutions varying from 200 m to 1 km. The order in which the model algorithms are executed are defined in the sphy.py file, in which the required modules are imported depending on the settings in sphy_config.cfg. Model settings (parameters, input and output) can be modified in the sphy_config.cfg configuration file. It is mandatory to have this file and the source code (or model executable) in the same folder on the pc's harddrive. If the user opts to run 10 the SPHY model using the model's source code, then the SPHY model is executed by entering "python.exe sphy.py" in the windows command prompt. Otherwise, the model can be executed by entering "sphy.exe" in the windows command prompt. Both the model source code and it's executable can be obtained from the SPHY model website (http://www.sphy.nl). 15 In order to run the SPHY model v2.0, it is required to have the following software installed on your pc: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Droogers, P., Immerzeel, W. W., Terink, W., Hoogeveen, J., Bierkens, M. F. P., van Beek, L. P. H., and Debele, B