HYDROGEIOS : a semi-distributed GIS-based hydrological model for modified river basins

The HYDROGEIOS modelling framework represents the main processes of the hydrological cycle in heavily modified catchments, with decision-depended abstractions and interactions between surface and groundwater flows. A semi-distributed approach and a monthly simulation time step are adopted, which are sufficient for water resources management studies. The modelling philosophy aims to ensure consistency with the physical characteristics of the system, while keeping the number of parameters as low as possible. Therefore, multiple levels of schematization and parameterization are adopted, by combining multiple levels of geographical data. To optimally allocate human abstractions from the hydrosystem during a planning horizon or even to mimic the allocation occurred in a past period (e.g. the calibration period), in the absence of measured data, a linear programming problem is formulated and solved within each time step. With this technique the fluxes across the hydrosystem are estimated, and the satisfaction of physical and operational constraints is ensured. The model framework includes a parameter estimation module that involves various goodness-of-fit measures and state-of-the-art evolutionary algorithms for global and multiobjective optimization. By means of a challenging case study, the paper discusses appropriate modelling strategies which take advantage of the above framework, with the purpose to ensure a robust calibration and reproduce natural and human induced processes in the catchment as faithfully as possible. Correspondence to: A. Efstratiadis (andreas@itia.ntua.gr)


Introduction
A central goal of the recent Water Framework Directive (2000/60/EU) is the establishment of river basin management plans.The plans will account, in detail, how the objectives set for each basin regarding its ecological, quantitative and chemical status are to be reached within the time horizon required, reviewing, among others, the impacts of human activities on the status of waters.Hence, the Directive emphasizes modified watersheds that are significantly affected by man-made interventions, structural or non-structural.In such basins, the inherent natural complexity of hydrological processes is amplified due to interactions between natural and artificial water bodies, on one hand, and structures regulated by man, on the other hand.Another source of complexity is the interaction between surface and ground waters, particularly when both are modified by withdrawals and disposal of previously used water.It becomes evident that a modelling attempt must account for the aforementioned interactions, through a simultaneous representation of the key physical processes (hydrological, hydrogeological, hydrodynamic, hydrochemical) and the water management practices.
Despite the wide expansion of hydrological simulation tools, the vast majority of them are applicable to natural basins.Hence, in modified watersheds, the usual practice consists in dividing the system under study into parts or sub-basins that can be modelled as natural ones.At a second stage, the outputs (i.e.river flows) are brought together, by adding man-made structures and their related processes.However, this two-stage procedure is not always feasible without drastic and, to a certain extent, unrealistic assumptions.For example, in case of conjunctive use of surface and ground water, the impacts of abstractions (e.g.baseflow reduction due to pumping) affect the entire downstream Published by Copernicus Publications on behalf of the European Geosciences Union.system, especially in leaky basins, with complex interactions between soil and aquifer.Besides, detailed historical data regarding abstractions are usually missing, especially when these are implemented through small private works (e.g.wells).Therefore, additional modelling is required to estimate abstractions, on the basis of theoretical water demands.
In the last years, several approaches have appeared on coupling hydrological models for surface and sub-surface flows (e.g., Singh and Bhallamudi, 1998;Panday and Huyacorn, 2004).Typically, however, these approaches do not represent all aspects of an operational water management problem at the river basin scale.Moreover, they preclude using stochastic simulation and forecasting, which are efficient methods to support water management.On the other hand, in most decision support systems (DSSs) for water resources management water flows are represented through network-type hydrosystems, thus ignoring the distributed regime of hydrological processes.In some cases, the problem is tackled by means of special elements in the hydrosystem, based on some form of elementary lumped models, such as the "groundwater reservoir node" in the RIBASIM software (Waterloopkundig Laboratorium, 1991).However, these elements, apart from being too simplified, contain parameters that normally require calibration (Nalbantis et al., 2002).Attempts to bridge this gap are very few, such as the MODSIM package (Fredericks et al., 1998;Dai and Labadie, 2001).
HYDROGEIOS is a new, GIS-based software system that provides a holistic framework, aiming at combining hydrological and hydrogeological simulation in modified basins.The model represents the governing interactions between surface flows, groundwater flows and man-made interventions, on the basis of a semi-distributed configuration.It integrates ideas from previous approaches (Nalbantis et al., 2002;Rozos et al., 2004;Rozos and Koutsoyiannis, 2006), whereas some components (e.g. the GIS module) are entirely new.HYDROGEIOS uses historical hydrological data for calibration and validation, as well as synthetic data for stochastic forecasting; in the last case, it co-operates with the DSS HYDRONOMEAS, which implements the optimization of the hydrosystem operation policy (Koutsoyiannis et al., 2002;Koutsoyiannis et al., 2003;Efstratiadis et al., 2004).Regarding the conceptualization, model parameters retain some physical consistency, since they are assigned on the basis of distributed data.For their estimation, the software encompasses a specific module, containing multiple fitting measures, statistical and empirical, and evolutionary algorithms for single-objective and multiobjective optimization.

Parameter uncertainty and calibration
Uncertainty is an inherent characteristic of all hydrological processes, further amplified by the weaknesses of determin-istic conceptual watershed models, whose parameters are estimated through calibration.Uncertainty is due to system complexity and multiple error sources, which are interacting in an unknown manner, thus making the traditional automatic calibration approach behave like a black-box mathematical game.Apart from evident errors in raw measurements and data-processing, typical sources of uncertainty are the inadequate representation of processes or, in the opposite, the formulation of too complex representations, unable to be supported by the existing knowledge about the physical system (Refsgaard, 1997;Wagener et al., 2001;Butts et al., 2004), the poor representation of the temporal and spatial variability of model forcing (Paturel et al., 1995;Chaubey et al., 1999;Beven, 2000;Andréassian et al., 2001), the nonrepresentativeness of calibration data and the use of statistically inconsistent fitting criteria (Sorooshian and Dracup, 1980;Kuczera, 1982;Sorooshian et al., 1983;Yapo et al., 1996;Gan et al., 1997), the poor identification of initial and boundary conditions (Kitanidis and Bras, 1980), the weaknesses of most optimization methods to handle response surfaces of irregular topography (Duan et al., 1992), as well as the temporal changes of natural and anthropogenic processes (Nandakumar and Mein, 1997;Brath et al., 2006;Ewen et al., 2006).The above problems have been thoroughly examined for more than three decades, concluding that uncertainty is inherent, thus unavoidable, and increases with model complexity.Uncertainty is also strongly related to the "equifinality" concept (Beven and Binley, 1992), practically identified as the existence of multiple acceptable parameter sets, on the basis of different model configurations, calibration data and fitting criteria.As a consequence of equifinality, it is impossible to detect a "global" optimal model structure or a "global" optimal parameter set, which definitely better reproduce the entire hydrological regime of a river basin.
In the last years, a variety of mathematical techniques were developed to quantify the uncertainty of conceptual model predictions.Most of them are embedded in the calibration procedure and seek "promising" trajectories of model outputs that correspond to multiple, "behavioural" parameter sets (Beven and Binley, 1992;Freer et al., 1996;Kuczera and Parent, 1998;Thiemann et al., 2001;Vrugt et al., 2002).Yet, their application indicates that, usually, the model predictive uncertainty proves comparable to the statistical uncertainty of the measured outputs.Moreover, many questions arise regarding the practical aspects of uncertainty analysis methods, such as the computational effort for multidimensional applications, the acceptance by policy makers and the public, and the inability to provide a final decision regarding a unique "best-compromise" parameter set (Pappenberger and Beven, 2006).
On the other hand, the applicability of physically-based models, which, in theory, would enable their parameters to be derived from field measurements, is significantly restrained by the heterogeneity of processes and the unknown scaledependencies of parameters (Beven, 1989;Wagener et al., Hydrol. Earth Syst. Sci., 12, 989-1006, 2008 www.hydrol-earth-syst-sci.net/12/989/2008/ 2001).This is the reason why many researchers tend to employ optimization for a small portion of the parameters (Refsgaard, 1997;Beven, 2001;Eckhardt and Arnold, 2001;Madsen, 2003;Vrugt et al., 2004;Muleta and Nicklow, 2005).In any case, integrating physically-based schemes within river basin management models has severe practical drawbacks, regarding the amount of spatial data required and the prohibitive computational effort.
Recent research revealed the advantages of conceptual semi-distributed models for streamflow estimation, in comparison to lumped ones (Boyle et al., 2001;Ajami et al., 2004).Such schemes allow a satisfactory representation of watershed heterogeneities, provide the required level of detail for an engineering application (due to the network-type configuration), while being computationally efficient.However, if interior calibration data are missing (e.g., hydrographs across the river basin), any movement from a lumped to a semi-distributed approach increases model complexity which, in turn, creates more uncertainty in the results.
HYDROGEIOS employs a semi-distributed scheme for the spatial representation of physical processes in modified hydrosystems.Instead of allocating parameters per subbasin, the parameterization of surface hydrological processes is implemented on the basis of hydrological response units (HRUs).The term was introduced by Flügel (1995) to characterize homogeneous areas with similar geomorphologic and hydrodynamic properties.The concept is widely used in distributed models, such as SWAT (Srinivasan et al., 2000), where the river basin is assumed to be an assembly of discrete entities with different characteristics that contribute differently to its responses.While a HRU is defined to serve a particular model conceptualization, commonly it denotes a spatial element of pre-determined geometry, identical to the schematization of the watershed.The main drawback of this approach is the huge number of unknown properties involved, which may be two or three orders of magnitude larger than the number of parameters of a lumped model, as indicated by Refsgaard (1997).To handle this problem in HYDROGEIOS, the HRU concept is used differently; it represents soil and land types, defining partitions of the basin, rather than "units" of contiguous geographical areas.In particular, the HRUs are defined as the product of separate partitions accounting for different properties such as soil permeability, land cover, terrain slope, etc.This product is formally known as common refinement of the partitions, while in the GIS terminology the related procedure is often called "union of layers".Through an appropriate classification of the above properties, one can adjust the number of HRUs and, consequently, the number of the parameters describing the soil hydrological mechanisms.Hence, parameters retain some physical meaning, which also allows a better identification of their prior uncertainty (i.e., the lower and upper bounds, used in calibration).Similar to surface water processes, the groundwater processes are represented through discretising the aquifer, by means of polygonal cells.There are no specific restrictions in the shape of cells, contrary to the majority of groundwater models, which implement a detailed, grid-based schematization.
The flexibility in the definition of the HRUs and the multicell representation of the aquifer allows the formulation of modelling schemes of variable complexity, depending on the available geographical and input data and the desirable parameterization.In contrast, the schematization of the hydrosystem, regarding the configuration of the physical and artificial network, is only restricted by the specified study requirements (i.e. the location of control points for the water balance calculations), and has no practical influence on the number of parameters.This is consistent with the principle of parsimony, which stands as a key point to handle uncertainties in model predictions.The simplest schematization and parameterization is ensured by using as many degrees of freedom as can be explained by the available "knowledge" about the system, regarding the hydrological and geographical data as well as the modeller's experience.
Results from previous research suggest that only five or six parameters can be identified on the basis of a single hydrograph (Jakeman and Hornberger, 1993); otherwise, parameter uncertainty due to poor identifiability may limit significantly the predictive capacity of models (Wagener et al., 2001).Thus, apart from a parsimonious configuration, it is critical to increase the amount of information in calibration, using multiple output variables and multiple performance measures, since a single measure would likely fail to reproduce all essential characteristics of the physical system that are reflected in the observations (Gupta et al., 1998).
Recently, there has been great interest in employing multiobjective optimization to better control the distributed model responses or specific aspects of them (Yapo et al., 1998;Madsen, 2000;Seibert, 2000;Beldring, 2002;Madsen, 2003;Muleta and Nicklow, 2005;Schoups et al., 2005;Tang et al., 2006).However, given the large number of parameters involved, a faithful implementation of such approaches requires extended hydrological measurements.Given also the scarcity of such information in most basins, an alternative is to use "soft" information, namely qualitative criteria indicating the acceptability of parameter values (Seibert and McDonnell, 2002).According to this, HYDROGEIOS implements multiobjective calibration, on the basis of typical statistical measures for the measured outputs (i.e.runoff and groundwater levels), criteria suitable for sparse measurements, and empirical criteria to control the internal model variables.Finally, for detection of a best-compromise parameter set, a hybrid strategy is proposed, based on a combination of automatic and manual techniques, which have proved very effective for complex hydrological models (Boyle et al., 2000;Rozos et al., 2004;Mazi et al., 2004 3 Model overview

Model formulation and input data
Since HYDROGEIOS is focused on the water resources management problem, the entire modelling approach is based on a network-type schematization of the physical and the artificial components of the hydrosystem.The two networks are linked together through water diversions, abstractions and disposals.The aquifer is also modelled as a network of conceptual storage (tanks) and transportation elements (conduits) connecting adjacent tanks.Most components have georeference and are handled through a GIS module, implemented in an ArcGIS environment.This module is also used for processing distributed data, used in the formulation of the HRUs, and the generation of the derivative layers (unions, intersections).

Hydrographic network
The schematization of the hydrographic network is implemented through a two-step procedure.First, an initial network is formulated on the basis of a digital terrain model, by adjusting the flow accumulation parameter.Next, additional control points are added across the network, which correspond to flow monitoring stations, abstraction points, inflow nodes, etc.The sub-basins upstream of each node are then created such that each river segment crosses a unique subbasin.
Hydrological inputs are precipitation and potential evapotranspiration time series, assigned to each sub-basin.For each sub-basin, the model calculates the transformation of precipitation to actual evapotranspiration, deep percolation and surface runoff; the latter is transferred as point inflow to the corresponding downstream node.The related processes are conceptualized via the soil moisture accounting model, described in Sect.3.2.1.All calculations are implemented on a derivative layer, generated as the product of the sub-basin and HRU layers, since meteorological forcing (precipitation and potential evapotranspiration) varies for each sub-basin.

Water management network
The applicability of HYDROGEIOS to modified basins is ensured through a coarse depiction of the major hydraulic works (pipes, channels, wells, etc.), the corresponding water uses and constraints and their interactions with the physical system.All are represented as network components, namely nodes and aqueducts; the latter may conduct water to the hydrographic network or abstract it to satisfy demands.Finally, wells lying on neighbouring locations and serving the same use are conceptualized as clusters, named borehole groups.
The network properties are discharge and pumping capacities, target priorities, demand time series and unit transportation costs.The priorities and costs are assigned to express preferences regarding the allocation of abstractions, in case of multiple water sources and conveyance paths.When a demand can be fulfilled through different abstractions, the user can impose unit costs (actual or hypothetical ones) to the corresponding aqueducts.For instance, a zero and a positive unit cost for surface and groundwater abstractions, respectively, will force the model to abstract water from the river rather than from groundwater.The preservation of target priorities and the minimization of costs are both ensured via the flow allocation model, explained in Sect.3.2.3.

Groundwater network
The aquifer is represented as a multicell network, on the basis of a polygonal discretization of the aquifer.According to Rozos and Koutsoyiannis (2005), the mathematical concept derives from the finite volume method, provided that the cell edges are parallel or normal to the equipotential contours and the line joining the centroids of adjacent cells is perpendicular to their common edge.This approach enables the exploitation of the available piezometric data for the study area.Moreover, the flexible number and shape of cells allows the description of aquifers of complex geometries on the basis of the physical characteristics of the system (e.g., geology) through parsimonious structures.Hence, the parameterization has a physical concept and the computational effort is significantly reduced, when compared to typical finite difference or finite element schemes.
Input properties are the top and bottom elevation of each tank and the water table at the beginning of simulation (initial condition of the model).Regarding boundary conditions, the user can prohibit the exchange of water between neighbouring cells, assuming an impervious common edge.The geometrical properties are automatically computed via the GIS module.These include cell areas, centroid coordinates, distances and common edge lengths between adjacent cells, plus all unions and intersections with the surface geographical layers.Note that some cells may lie out of the watershed bounds, to direct the groundwater sinks to the sea or neighbouring basins.Other components are springs, which represent point outflows that are transferred to the downstream node of the corresponding sub-basin; their properties are the altitude and the interconnected cell.

Mathematical framework
The mathematical representation of the hydrological, hydrogeological and anthropogenic processes is based on the combination of three related models running within a loop, as explained in Sect.3.2.4.The time scale of simulation is monthly, which is sufficient for water management studies.

Surface hydrology model
The hydrological processes above and across the unsaturated zone are modelled through a conceptual water balance Hydrol.Earth Syst.Sci., 12, 989-1006Sci., 12, 989- , 2008 www.hydrol-earth-syst-sci.net/12/989/2008/ scheme, illustrated in Fig. 1.Model inputs are the areal precipitation P t and potential evapotranspiration E P t , while outputs are the soil moisture S t , the surface runoff Q t , the actual evapotranspiration E t , and the percolation G t .At each time step t, the water balance equation is written in the form: For a given value of soil storage at the beginning of simulation, the above formula is solved on the basis of some assumptions regarding the unknown variables Q t , E t and G t .
The ground operates as a filter, transforming precipitation to direct evapotranspiration E Dt , direct runoff Q Dt , and infiltration I t .Direct evapotranspiration represents the amount of precipitation evaporated quickly, before infiltrating, and cannot exceed a retention capacity R, and the theoretical demand E P t .Direct runoff represents the excess of precipitation conducted through the impervious areas of the basin to its outlet within the time interval, and calculated as Q Dt =c(P t −E Dt ), where c is a constant ratio that depends on the physical properties of the ground (soil permeability, vegetation, slope) and the existence of flood-prevention works.
The remaining precipitation infiltrates to the unsaturated zone, represented as a soil moisture accounting tank of capacity S max .The tank is divided into two zones (upper and lower), using a dimensionless parameter κ.The evapotranspiration deficit, i.e. the amount E P t −E Dt , is satisfied by the actual moisture, using different mechanisms for the two zones.Specifically, the whole amount of moisture in the upper zone is assumed available for evapotranspiration, whereas the lower zone moisture is partially available.In the last case, the rate of soil evapotranspiration is taken to be proportional to the ratio S Lt /(κS max ), where S Lt is the moisture depth stored in the lower zone and κS max is the corresponding capacity.The process is mathematically expressed by a declining exponential function, similar to that of the wellknown Thornthwaite model (Thornthwaite, 1948).
Moreover, the soil moisture tank provides options for horizontal and vertical outflow, implemented via two orifices, the one lying at level κS max and the other at the bottom.These represent a time-lagged runoff component (interflow) Q I t and the percolation to deeper zones G t , respectively.The corresponding outflow rates are controlled through the retention coefficients λ and µ.
At the end of the simulation step, the soil moisture excess (if it exists) contributes to the streamflow as quick runoff due to saturation, Q St .Within the time interval the storage is allowed to exceed S max .
Model parameters are the retention capacity R, the direct runoff coefficient c, the soil capacity S max , the ratio of the lower zone capacity to the total one κ, and the retention coefficients λ, µ for interflow and percolation, respectively.These differ for each HRU, as explained in Sect.3.1.1.All variables are integrated to the sub-basin scale, except for the percolation, which is integrated to the groundwater cell scale.The model then adds to the surface flow, the direct, quick and 37   time-lagged runoff components, and the baseflow arriving from the springs located in the sub-basin.A percentage of the total runoff, equal to an infiltration coefficient δ assigned to each river segment crossing the sub-basin, recharges the groundwater system, whereas the rest is conducted to the outlet node, as point inflow to the hydrographic network; the coefficient δ is either pre-specified or estimated through calibration.

Groundwater model
Each groundwater cell is represented by a conceptual tank, whose parameters are the specific yield (dimensionless) and the conductivity, expressed in velocity units.The stress components of groundwater tanks are: (a) areal inflows due to percolation through each sub-basin and HRU combination; (b) inflows due to infiltration underneath each river segment; and (c) point outflows due to pumping from each well.We remind that percolation rates are output of the surface hydrological model, whereas infiltration and pumping rates are output of the water management model.Regarding percolation, the model integrates the equivalent depths from each sub-basin and HRU combination on the corresponding cell area.Regarding infiltration, the model estimates the river segment losses supplying each tank, assuming that: where I j is the sum of infiltration losses through the river segment j , L j is the segment length and L ij is its partial length over cell i.
For given stresses, the flow field problem is solved according to a simplified version of the scheme introduced by Rozos and Koutsoyiannis (2005), which proved suitable for simulating aquifers of high parameter uncertainty (e.g., karst where w i is the water level in tank i, w min i and w max i are the bottom and top absolute levels, respectively (w max i =w min i + b i ), θ is the ratio of specific yield to confined storage coefficient and b i is the layer thickness.The upper equation in (3) corresponds to phreatic conditions, while the lower corresponds to confined conditions, so that b i represents also the threshold between confined and unconfined conditions.The water volume contained in the tank equals the base area of the corresponding cell multiplied by the level and the specific yield; a low value of the latter indicates that a large water level increase is required to store a particular amount of groundwater, and vice versa.
A constant head condition is represented by assigning tanks with very large base, which forces the corresponding water level to remain practically constant and close to the prescribed boundary value.Likewise, springs are modelled assuming such dummy tanks, for which the slight changes of the water level are directly transformed into outflow hydrographs.A similar representation is implemented for simulating groundwater losses, conducted to neighbouring basins or the sea.
Groundwater flows are implemented through conceptual conduits (i, j ), where the indices denote the interconnected tanks.Their properties are the cross sectional area A ij (equal to the common plane area between cells i and j , which is assumed constant within each time interval), the length l ij (centroid distance) and the conductivity K ij , computed as the arithmetic or geometrical mean of the corresponding tank conductivities.The discharge Q ij is calculated using a Darcian formula: where h i and h j are the head values of the adjacent tanks i and j .Equations ( 3) and (4) formulate a system of equations that can be solved via explicit or implicit numerical schemes.HYDROGEIOS implements both schemes, for which a proper time discretization (i.e.number of time intervals within a simulation step) must be defined, to ensure numerical stability.

Water management model
Outputs of the surface and groundwater hydrological models are the sub-basin and spring runoff, both assumed as point inflows to the river network.The water allocation is implemented on the unified network, to define the unknown fluxes through the entire hydrosystem.These include the discharge rates and losses (due to leakage) across the river and the aquifer and, subsequently, the abstractions from surface and groundwater resources.
The modelling concept is based on a linear programming (LP) approach.Similar ones have been used in some water resource planning and management applications, where linear optimization is embedded within simulation to find the least cost flow allocation through hydrosystems of network format (Graham et al., 1986;Kuczera, 1989;Fredericks et al., 1998;Dai and Labadie, 2001).The optimization is based on real economic criteria or artificial costs, assigned to preserve water rights and water use priorities.
HYDROGEIOS implements a simplified version of the scheme described by Efstratiadis et al. (2004).The general idea is to distinguish the hydrosystem variables and optimally allocate them through the hydrosystem, which is represented as a digraph.Apart from real world components, the digraph includes dummy nodes and links where virtual attributes are imposed, namely the conveyance capacity and the unit transportation cost.The latter may be either positive or negative.Particularly, positive unit costs are imposed to penalize non-desirable water fluxes, whereas negative unit costs are imposed to force the model to provide water to fulfil the physical and operational constraints.Specifically: -Leaky river segments are represented by two links, the one carrying the discharge arriving at the downstream node, and the other carrying the infiltration, which is transferred to an accounting node, a fictitious component inserted for mathematical convenience (to ensure that the sum of inflows equals the sum of outflows).
-River segments or aqueducts, where minimum flow preservation targets are imposed, are represented by two parallel links, the one having discharge capacity equal to the actual target value and negative cost, whereas the other has the rest of capacity (infinite in the case of a river segment) and unit cost equal to the real transportation cost (zero in case of river segment).Maximum flow targets are handled in the same way; positive cost is used here to prevent the violation of the discharge bound.
-Borehole groups are represented through a virtual groundwater node, the inflow of which is equal to the total pumping capacity of the wells.Two links are connected to this node, the one carrying the groundwater abstraction to the corresponding downstream node with unit cost equal to the total pumping charge, and the other transferring the rest of inflow to the accounting node, without cost.
-Demand nodes are connected with the accounting node via a virtual link, the capacity of which is set equal to the actual demand rate, while a negative unit cost is imposed to force the model to satisfy the corresponding target.
The transformation of hydrosystem components to digraph components and the assignment of capacity and unit cost values are automatically executed by the program.Most formulations are done once, at the beginning of each simulation.Especially, the assignment of costs is a key part of the model, since this ensures the preservation of both feasibility and economy.Feasibility refers to the strict satisfaction of all physical constraints (nodal water balance equations and arc capacity bounds) and the hierarchical satisfaction of water uses, keeping the user-defined priorities, whereas economy refers to the minimization of total transportation and pumping costs.
The problem is expressed as: where x is the vector of control variables, corresponding to the hydrosystem fluxes; c is the vector of unit costs; A is the incidence matrix, describing the continuity equations, with elements taking values {−1, 1, 0}; is a matrix describing constraints for leaky segments, with elements taking values , 0}; y is the vector of inflows; and u is the vector of link capacities.Due to the particular structure of (5), primarily the sparse format of matrix A, its solution is very fast through appropriate versions of the simplex method, thus ensuring computational efficiency.

Model integration within simulation
Due to the interactions between surface and groundwater resources, as well as the physical and man-made processes, the application of the aforementioned models within simulation requires a looped architecture, as illustrated in Fig. 2. At the beginning of each time step, dynamic input data includes precipitation and potential evapotranspiration depths, assigned at each sub-basin, and water demand values, assigned at specific nodes of the hydrosystem.The remaining hydrological variables are unknown, and for some of them initial guesses are necessary.The simulation procedure is implemented as follows: First, and outside of the loop, the surface hydrological model runs to estimate real evapotranspiration, percolation, surface runoff and soil moisture for each combination of subbasins and HRUs.Runoff is transferred to the outlet node of each sub-basin, after adding baseflow and, then, excluding losses due to infiltration.Baseflow is computed by adding discharge values of springs lying in each sub-basin.Initially, these are assumed equal to the values of the previous time step, but as the loop proceeds, the real ones are assigned, as estimated from the groundwater model.
Next, the water management scheme runs to estimate all hydrosystem fluxes, i.e. discharge and infiltration val- ues across the hydrographic network, and abstractions from surface and groundwater resources.Inflows to the digraph model are runoff values, assigned downstream of each subbasin, already known from the surface hydrological model, as well as external inflows, given as known time series.The implementation of the above models allows the assignment of groundwater stresses at each cell.These are estimated by adding percolation from each sub-basin and HRU combination and infiltration from each river segment, and excluding pumping from each well.Next, the groundwater flow model is solved, to estimate the tank levels, the spring flows and the underground losses.
Based on the actual evaluation of spring discharge, HY-DROGEIOS recalculates baseflow and corrects all runoff estimates.This requires new runs of the water management and groundwater models, until stabilization of baseflow.Practically, this scheme converges after one or two cycles only, thus ensuring both accuracy and efficiency.

Fitting measures
The mathematical framework described herein comprises a variety of parameters, illustrated in Table 1.Specifically, it uses one parameter per river segment, six parameters per HRU, and two parameters per groundwater tank, while the water management module does not contain free variables to optimize.Thus, even for a relatively small hydrosystem, some dozens of control variables would be involved.Generalising the empirical rule mentioned in Sect.2, initially used in lumped rainfall-runoff models, multiple criteria should be introduced, to ensure consistency with the principle of The statistical measures used are the coefficient of efficiency, also known among hydrologists as the Nash-Sutcliffe index (Nash and Sutcliffe, 1970), and the bias in the mean, the standard deviation and the coefficient of variation of measured responses; the latter may refer to both discharge and level time series.Yet, due to the usually rough delineation of the groundwater field, resulting to large cell areas and thus relatively small variability of the corresponding simulated levels, one should employ very carefully calibration on the basis of observed water table data, due to the issue of scale compatibility.
Furthermore, the model uses an empirical penalty term, to better control the intermittencies, since a "zero discharge" state is very frequent, and its reproduction is important for a water management study.Given the observed and simulated time series, y t and y t , respectively, the penalty is calculated as: where z t is an auxiliary variable, computed as: and T 0 the number of time steps for which the model fails to reproduce an observed flow interruption or, in the opposite, erroneously yields zero discharge.
A final measure is used to control the behaviour of the internal model variables, specifically the generation of unreasonable trends regarding groundwater levels, based on the Mann-Kendall rank correlation test (Kottegoda, 1980, p. 32-34).When attempting to calibrate the groundwater parameters (i.e., conductivities) merely on spring hydrographs, without using observed level data, a conjunctive model could easily preserve the water balance of surface flows by leaving some upstream tanks empty and, simultaneously, accumulating the excess of water downstream.This situation is not consistent with the physical behaviour of an aquifer, the level of which follows the typical seasonal and overyear fluctuation of precipitation.However, in heavily modified basins with intensive exploitation of groundwater, a systematic decline of the water table could occur.Thus, even if level data is totally missing, an appropriate use of the trend penalty should significantly improve the identifiability of the groundwater parameters, thus leading to more reliable schemes.
The Mann-Kendall test is implemented as follows: Given a sample (x 1 , x 2 , ..., x N ), the statistic T =r/ σ 2 r is a standard normal variable, where: and P is the number of all pairs {x i , x j , j >i} with x i <x j .For a two-tailed test and for a level of significance a, we reject the null hypothesis of no trend presence if |T |<z a/2 .In that case, a penalty value equal to |T |−z a/2 is assigned.

The evolutionary annealing-simplex algorithm
The evolutionary annealing-simplex algorithm is a probabilistic heuristic global optimization technique, combining the robustness of simulated annealing in rough response surfaces, with the efficiency of hill-climbing methods in convex areas.The version used in HYDROGEIOS differs slightly from the ones presented by Efstratiadis and Koutsoyiannis (2002) and Rozos et al. (2004), which proved effective and efficient for a variety of hydrological applications, including calibration problems.An innovation is the assumption of two parameter ranges; the interior one is used for the generation of the initial Hydrol.Earth Syst.Sci., 12, 989 -1006, 2008 www.hydrol-earth-syst-sci.net/12/989/2008/ population, but can be violated through the evolution, whereas the exterior one is rigid, and represents the feasible space.The corresponding interior bounds represent initial guesses for parameters, based on their physical interpretation, while the exterior ones express their physical bounds.During one generation, the population evolves as follows: First, a simplex-based pattern is formulated, using random sampling.Next, a candidate individual is selected to die, according to a modified objective function of the form: where f is the original objective function, T is the current "temperature" and u is a random number from the uniform distribution.The temperature is gradually reduced, according to an appropriate annealing cooling schedule, automatically adapted during the evolution.Consequently, the probability of replacing individuals with poor performance increases, since the procedure gradually moves from a random walk to a local search.
The recombination operator is based on the well-known downhill simplex transitions (Nelder and Mead, 1965).According to the relative values of the objective function at the vertices, the simplex is reflected, expanded, contracted or shrinks, where quasi-stochastic scale factors are employed instead of constant ones.To ensure more flexibility, additional transformations are introduced, namely multiple expansion towards the direction of reflection, when a downhill path (i.e., the gradient of the function) is located, and similar expansions but on the opposite (uphill) direction, in order to escape from the nearest local minimum.If any of the above transitions improves the function value, the new individual is generated through mutation.The related operator employs a random perturbation scheme outside of the usual range of the population, as determined on the basis of the average and standard deviation values of its coordinates.

Multiobjective version
Recently, Efstratiadis (2008) developed a multiobjective version of the above scheme, suitable for challenging calibration problems, where multiple responses are to fit on multiple criteria (Efstratiadis and Koutsoyiannis, 2008).The multiobjective evolutionary annealing-simplex is also embedded in the software, although its full description is out of the scope of this study.The source code of both the single-and multiobjective optimization algorithms is available on-line at www.itia.ntua.gr/en/docinfo/838/.
The algorithm embodies two phases; in the evaluation phase a fitness measure is assigned to all population members, whereas in the evolution phase new individuals are generated on the basis of their fitness values.The fitness measure aggregates various terms, to guide the search towards non-dominated solutions, to provide well-distributed populations and to penalize non-dominated parameter sets with extreme performance (i.e., too good against some criteria, but   too bad against some other).Thus, the most promising part of the Pareto front is approximated, in an attempt to surround a best-compromise solution.
Regarding the evolving phase, some aspects are similar to the transitions used in the single-optimization approach, whereas some alterations are necessary to prohibit population convergence (e.g., the simplex is not allowed to shrink).Moreover, the mutation operator employs two schemes, with equal probability; one allows small perturbations around the candidate individual to die, while the other ensures the generation of a random solution outside of the average range of the population.

The study area
The Boeoticos Kephisos river basin lies in the Eastern Sterea Hellas, north of Athens, and drains a closed area (i.e., without an outlet to the sea) of 1956 km 2 (Fig. 3).The catchment geology comprises heavily karstified limestone, mostly developed on the mountains, and alluvial deposits, lying in the plain areas.Due to its karst subsurface, the watershed has a considerable groundwater yield.The main discharge points are large springs in the upper and middle part of the basin that account for more than half of the annual catchment runoff.Moreover, an unknown amount of groundwater is conducted to the sea.
The basin is significantly modified, since it serves multiple and contradictory water uses.Specifically, through an extended drainage network, the entire surface resources are diverted to the neighbouring Lake Hylike (one of the major water storage projects of Athens), through a canal and a tunnel.Besides, important supply boreholes are located at the middle course, just upstream of the Mavroneri springs; these are activated in case of emergency, and affect significantly the flow regime of the hydrosystem.In addition to drinking water, significant surface and groundwater resources of the basin are used for irrigation.During the summer period, all surface water is used for irrigation, thus drying the canal at the basin outlet; in addition, part of the demand is satisfied via pumping from Hylike.
The estimation of the water balance of the basin on the basis of runoff data is impossible, because of the groundwater losses, the large amounts of water infiltrating across the upper course of the river (a 25% reduction of discharge is detected, according to a series of flow measurements), and the existence of combined abstractions.Previous attempts, thought a simplified version of the model, with lumped description of the main processes (Rozos et al., 2004), indicated that, due to the unknown distribution between evapotranspiration and sea outflows, the problem is ill-posed, and an infinite number of solutions exist, providing similar performance.In the present approach, we tried to establish a much more "physical" scheme, to enhance the information content in calibration.

Model formulation and data
As illustrated in Fig. 3, the river network comprises a main branch, divided in four segments, and five sub-basins upstream of or between the corresponding nodes; the main properties of the sub-basins are summarized in Table 2.For the delineation of soil processes we defined six HRUs by combining two geographical layers representing three categories of permeability (low, medium, high), and two categories of terrain slope, with threshold 10% (Fig. 4; see also Table 4).These three categories of permeability were extracted through aggregating detailed hydrolithological data.The low-permeable areas of the basin (21.6%) are mainly constituted by flysch, which is waterproof, the mediumpermeable ones (42.8%) are generally captured by quaternary (alluvial) deposits, while the high-permeable areas (35.6%) are constituted by limestones of different phases, including highly karstified formations.To restrict redundant parameters, we didn't incorporate other types of information, e.g., land use, which practically overlaps with slope, since mountainous areas are covered by forests, while plain ones are dominated by crops and, generally, low vegetation.The groundwater flow field is divided in 35 non-rectangular cells (Fig. 5); six of them implement surface outflows through the major karst springs, while two are located outside of the basin to simulate the draining of underground leakages to the sea.The spatial distribution of springs, with regard to the surface and groundwater system delineation, is illustrated in Figs. 3 and 5.
The water management network, sketched in Fig. 6, includes four conceptual nodes that represent extended agricultural areas (totally 348 km 2 ), six borehole groups and a dozen of aqueducts conducting abstractions from the river and the aquifer to the related nodes.Since some demands are fulfilled via multiple sources, virtual costs are assigned to the corresponding aqueducts thus representing the real management policy (i.e.priority in using surface resources, instead of the groundwater ones).Additional targets are water supply through the middle course boreholes that were drilled during the early 1990's.Historical abstractions from Hylike are imported to the network as external inflows, with known values (28.9 hm 3 on annual average).
For the above schematization, the total number of parameters is more than 100.Taking into account the aim of the study (i.e.combined simulation of the most dominant surface, groundwater and man-induced processes at multiple sites) and the available data (discharge measurements and additional observations, as explained below), we consider that this number is justifiable for realistically representing the basin's heterogeneities and uncertainties.Besides, it is essential to embed multiple criteria within calibration, to avoid over-parameterization and to properly represent all important characteristics of the physical system that are reflected in the observations.The latter refer to systematic (daily) discharge measurements at the basin outlet (Karditsa tunnel) and sparse (two per month) measurements downstream of the six springs.These samples were used to construct monthly hydrographs at seven discharge points, for a 10-year period (October 1984-September 1994), which was the control horizon of the study.Plotted data in Figs.7-13 illustrate the irregular behaviour of most hydrographs, which reflects the high complexity of the hydrosystem.Unfortunately, it was impossible to use similarly sparse discharge data along the river, which would provide valuable information about     the allocation of runoff at the sub-basin scale, thus improving the identifiability of parameters of the surface hydrological model.
With respect to groundwater level, about 20 gauges were available for the aforementioned period, mostly located in the vicinity of Boeoticos Kephisos.However, due to the difference of scale between point observations and averages over the cell areas, superimposed to high heterogeneity and uncertainty of the karst aquifer, we preferred not to include this information in calibration, i.e. we didn't attempt any kind of "adjustment" of the simulated levels to the piezometric information.Thus, the estimation of groundwater model parameters was based on the observed discharge data downstream of the six springs, the empirical criteria used to avoid unrealistic trends of the simulated aquifer levels and through rough visual inspection of the groundwater outputs within the hybrid calibration procedure (see next section).The remaining hydrological inputs were monthly precipitation and potential evapotranspiration time series, integrated on the surface of the five sub-basins.The former were constructed via the Thiessen method using point data from 12 rainfall stations, well-distributed through the basin, whereas the latter was estimated via the Penman-Monteith method, on the basis of land-use information.Theoretical irrigation needs were approximated by assuming an average annual value of 6500 m 3 /ha of irrigated land, which corresponds to an annual demand of 226 hm 3 with unknown (a priori) allocation.

Calibration strategy
The model parameters (∼100) were fitted on about 40 criteria, weighted in a single performance measure.These include: (a) the efficiency index and the average bias, to calibrate the hydrographs at the basin outlet (Karditsa tunnel) and downstream of the six karst springs; (b) the zero-flow penalty to better fitting the discharge at the outlet (systematically going to zero during each irrigation period) and downstream of the Mavroneri springs (with zero flow, in case of intensive pumping); and (c) the trend penalty to realistically represent the variability of groundwater cell levels, except for those lying in the neighbourhood of springs, the behaviour of which is better controlled though the hydrographs.The first two statistical criteria ensured satisfactory spread around the measured values and avoidance of systematic errors, while the empirical measures enhanced the information embedded in calibration (thus ensuring compatibility between model free variables and criteria involved) and helped to control much larger number of model responses than the observed ones.The 10-year control period was split in a six-year calibration period (October 1984-September 1990) and a fouryear validation period (October 1990-September 1994).
Given the intricacy of the derived optimization problem, it was impossible to obtain a reliable solution in a single run.Apart from the vast number of local optima and the irregularities of the response surface, an additional complexity factor was the different order of magnitude between the groundwater conductivities, taking values in the range 10 −8 to 10 −1 , and the rest of parameters, most of them being dimensionless.This was partly remedied using logarithms of conductivities.
To deal with the multiple puzzles of the calibration problem, while ensuring a satisfactory predictive capacity of the model, a hybrid calibration strategy was employed, through progressive improvements of relatively small groups of parameters.At a preliminary phase, we used extended bounds for the search space and tried various combinations of weighting coefficients, to obtain a general overview of the problem, regarding the multiple criteria interactions and their feasible range.
At the second phase, we attempted to optimize the HRU parameters, as well as the most important parameters of the groundwater model, specifically the conductivities of the springs and their adjacent cells.Moreover, the interior bounds of parameters were restricted to be consistent with their broad physical interpretation.The main objective was to guarantee a good fitting of the hydrograph at the outlet, especially its parts related to high flows, and a satisfactory fitting of the spring flows.At the end of this phase, we removed most trend penalties, since we ensured a "regular" behaviour of the groundwater model, by appropriately adjusting the corresponding parameters.
In the last phase, starting from a relatively good solution, we focused on improving specific aspects of the model responses.This proved not an easy task, since even slight ameliorations of one criterion had asymmetrically negative impacts on others, due to the high sensitivity of some parameters.Often, it was necessary to accept non-optimal transitions, regarding the overall value of the objective function, to ensure the improvement of particular aspects of the hydrographs.In that phase, we focused on the predictive capacity of the model, accounted on the basis of efficiency values in validation, and the consistency of parameters.The latter was evaluated according to our experience, recognising the fact Hydrol.Earth Syst.Sci., 12, 989-1006Sci., 12, 989- , 2008 www.hydrol-earth-syst-sci.net/12/989/2008/    that the model is unavoidably vulnerable to the existing multiple sources of uncertainty, which is further amplified by the high complexity of the system under study.

Results and discussion
The optimized statistical measures against the simulated runoffs are summarized in Table 3, whereas Figs.7-13 compare the observed and the computed hydrographs at the seven control sites.Regarding the runoff at the outlet (Fig. 7), a very good fitting is ensured for both the calibration and validation periods, with efficiency values 87.0% and 75.6%, re-   spectively.The model preserves the important aspects of the hydrograph, namely the high flows and the artificial flow interruption during the summer due to upstream abstractions.Moreover, it reproduces the sequence of high and low flow periods, which is more prominent during the validation period.Analysis with extended historical data further validated the model capacity relating to the prediction of the basin runoff (Koutsoyiannis et al., 2007).
For the Lilaia-Kefalovrysso springs, located in the upper course of the basin, the model provides a very satisfactory performance (efficiency 80.6% in calibration, 60.7% in    validation), apart from a deviation during the winter of 1992, which was characterized by large amounts of snow in contrast to low liquid precipitation depths.
For the Mavroneri springs, located on the middle course and very close to the water supply wells, the fitting was also very satisfactory (efficiency 69.3% in calibration, 60.1% in validation).Indeed, the model represents the two characteristic periods of flow intermittency, where the first (May-December 1990) is entirely due to the persistent drought of the late 1980's, whereas the second one, lasted more than a year (end of 1992 to start of 1994), resulted as combination of unfavourable hydrological conditions and intensive use of the newly constructed boreholes.Between these extremely dry periods, an impressive increase of discharge was observed, well-represented by the model.In general, the overall fitting on this particularly important hydrograph was a major guarantee of the model performance, especially when taking into account the high uncertainty of such a karst system.
Regarding the other springs, the model fitting was less satisfactory, although acceptable.The efficiency values achieved vary from 72.4% for the Agia Paraskevi springs (the less important of the whole system) to 26.5% for the Melas springs.The latter contribute significantly to the total basin runoff and their mechanisms are extremely complex.Previous modelling attempts failed, even when detailed tools were used.For example, a simulation based on the MODFLOW achieved an efficiency value of 10% (Nalbantis et al., 2002), whereas the lumped approach of Rozos et al. (2004)   a value of 19.4%, regarding the combined representation of Melas and Polygyra springs (the corresponding value for the validation period was slightly negative).Hence, the specific approach, based on a rough discretization of the groundwater field and the use of linear (i.e.Darcian) equations for describing water interchanges between groundwater tanks, worked better than the generic MODFLOW, mainly regarding the preservation of the observed water balance (as indicated by the negligible bias values).Yet, the lack of extended spatial information and the uncertainties of natural mechanisms did not allow a more thorough schematization and parameterization, which would probably (but not definitely) improve the model performance.
The physical interpretation of parameters is generally difficult, especially for those of the groundwater model, because of the complexity of the karst system.Table 4 shows the optimal values of the six parameters of the soil moisture model, assigned to the corresponding HRUs.It is not surprising that the direct runoff coefficients, c, and the retention rates for percolation, µ, are mainly affected by the soil permeability, whereas the soil capacities, S max , are more related to the terrain slope.Thus, the plain areas of the basin have almost twice the capacity of the mountainous ones, for the same category of permeability.On the other hand, the percolation rates through the high-permeable soils are significant, which explains the limited contribution of surface flows to the total water potential of the basin.Regarding the river infiltration coefficients, their optimal values are 26.4,8.5 and 3.1%, along the upper, middle and lower course of Boeoticos Kephisos, respectively.The significant percentage of water losses in the upstream segment is consistent with the flow measurements, as mentioned in Sect.5.1.
The semi-distributed formulation of the model provides a much clearer view of the water balance of the basin and the spatial distribution of its water resources (Table 2).The overall water balance components, for the 10-year control period, are summarized in Table 5. Significant part of precipitation is lost due to evapotranspiration (62.3%) and underground runoff, conducted to the sea (10.4%).The results are close to the first calibration scenario reported by Rozos Hydrol.Earth Syst.Sci., 12, 989-1006Sci., 12, 989- , 2008 www.hydrol-earth-syst-sci.net/12/989/2008/ et al. (2004).Percolation reaches 29.6% of precipitation, whereas surface runoff is only 6.9%, due to the dominance of highly permeable areas.Due to increased demands, less than half of the available resources reach the basin outlet, thus indicating a significantly modified hydrosystem.From the 233.0 hm 3 of the estimated annual demand for the specific period (including non-systematic demands for the water supply of Athens), 28.9 hm 3 (12.4%)are externally provided through pumping from Hylike (this was the only value known a priori), 132.6 hm 3 (56.9%)are withdrawn from inbasin boreholes and 71.5 hm 3 (30.7%)are fulfilled from surface abstractions.

Summary and conclusions
HYDROGEIOS integrates a conjunctive surfacegroundwater simulator within a water management scheme, to describe the hydrological processes and the impacts of human interventions.It is suitable for simulating the water balance across modified hydrosystems with conjunctive water uses even with limited data.It aims to treat the issues mentioned in the introduction, such as the faithful representation of the decision-related interactions, through establishing computationally efficient modelling schemes.Other significant issues are the use of GIS for generating various levels of spatial information, the physical consistency related to the modelling components, and the parsimonious use of parameters.The software provides also tools for the calibration of parameters, on the basis of multiple fitting criteria and advanced algorithms for singleand multi-objective optimization.
The key points of our approach were demonstrated through a case study, involving a real management problem in a hydrosystem of many peculiarities.These refer to both the physics of the basin (karst subsurface, high contribution of baseflow to total runoff, significant lateral outflows to the sea) and the water management regime (combined supply from surface and groundwater resources, multiple water uses, negative impacts of pumping on the downstream water availability, lack of real abstraction data).The calibration was based on a combined strategy, where the hydrological experience had the key role and optimization, carried out through the evolutionary annealing-simplex method, was used as an auxiliary tool, which provided fast solutions.Having a leaky basin with ill-posed boundaries and thus unknown allocation of hydrological losses due to evapotranspiration and underground outflows, the hundred parameters to optimize and seven hydrographs to fit, the equifinality problem emerged.However, we attempted to increase the information contained in calibration data, by assigning additional criteria to control specific aspects of the measured outputs (e.g., the reproduction of flow interruptions) and criteria to prohibit the generation of unrealistic trends, regarding the internal variables of the groundwater module, which is another innovation in our approach.Moreover, by emphasising the model performance in validation and the physical interpretation of parameters, we managed to guide the search towards a best-compromise set, which is essential for an engineering application.
Further analysis is now implemented regarding three key issues of the proposed modelling framework.The first involves a sensitivity analysis of parameters accompanied by the examination of different levels of delineation in the study area, in the direction of reducing the number of free variables, thus gaining on model parsimony, without loosing on performance.The second focuses on the practical use of Pareto-based approaches, in combination with classical global optimization, to analyze possible structural and data errors and detect the most promising areas of search spaces in problems involving many parameters and criteria.A last research issue under way is the model implementation on finer time scales, by means of incorporating routing procedures within simulation.The task is not straightforward, given the co-operation of three modules interchanging inputs and outputs (surface hydrological model, groundwater model and water management model); an iterative procedure is required within each time step, which is, however, inconsistent with the condition of successive time periods assumed in all known numerical routing schemes.The results of these investigations will be reported in due course.

Figure 2 :
Figure 1: Schematic layout of the surface water balance model.

Fig. 1 .
Fig. 1.Schematic layout of the surface water balance model.

Figure 1 :Figure 2 :
Figure 1: Schematic layout of the surface water balance model. 38

Figure 3 :
Figure 3: The Boeoticos Kephisos river basin and the main hydrosystem components.

Figure 4 :
Figure 4: Characteristic layers of geographical data for the schematization of the surface system: (a) top left, terrain slope; (b) bottom left, permeability; (c) top right HRUs, derived as the product of slope and permeability; (d) bottom right, river network and sub-basins.

Figure 5 :
Figure 5: Characteristic layers of geographical data for the schematization of the groundwater system: (a) left, formulation of cells, based on permeability; (b) right, product of cells, subbasins, HRUs, springs and boreholes.

Fig. 4 .
Fig. 4. Characteristic layers of geographical data for the schematization of the surface system: (a) top left, terrain slope; (b) bottom left, permeability; (c) top right HRUs, derived as the product of slope and permeability; (d) bottom right, river network and sub-basins.

Figure 5 :
Figure 5: Characteristic layers of geographical data for the schematization of the groundwater system: (a) left, formulation of cells, based on permeability; (b) right, product of cells, subbasins, HRUs, springs and boreholes.

Fig. 5 .
Fig. 5. Characteristic layers of geographical data for the schematization of the groundwater system: (a) left, formulation of cells, based on permeability; (b) right, product of cells, sub-basins, HRUs, springs and boreholes.

Figure 6 :Figure 7 :
Figure 6: Simplified sketch of the water management network, representing abstractions from the surface and groundwater resources.Nodes 1-5 are river points, whereas nodes 7-10 denote agricultural areas across the basin.

Fig. 6 .
Fig.6.Simplified sketch of the water management network, representing abstractions from the surface and groundwater resources.Nodes 1-5 are river points, whereas nodes 7-10 denote agricultural areas across the basin.

Figure 6 :Figure 7 :
Figure 6: Simplified sketch of the water management network, representing abstractions from the surface and groundwater resources.Nodes 1-5 are river points, whereas nodes 7-10 denote agricultural areas across the basin.

Figure 13 :
Figure 12: Computed and observed discharge series at Melas springs.

Table 1 .
List of model parameters.

Table 2 .
Sub-basin properties.Computed discharge at the downstream node, comprising surface (basin) and groundwater (spring) runoff minus river abstractions.

Table 3 .
Optimal values of efficiency and relative bias for the calibration and validation periods. attained