Assessing water resources availability - Deductive methods

Deductive assessment of water resources availability is carried out by applying mathematical models which simulate the water flow in the surface river network or underground aquifers. Modeling of surface water is usually carried out through hydraulic models or rainfall-runoff models. Aquifers are modeled by applying groundwater models.

1. Rainfall-runoff modeling

Rainfall-Runoff modeling is one of the most classical applications of hydrology. It has the purpose of simulating a river flow hydrograph induced by an observed or a hypothetical rainfall forcing. Rainfall-runoff models may include other input variables, like temperature, information on the catchment or others. Within the context of this subject, we are studying rainfall-runoff models with the purpose of producing simulation of long term hydrographs. Rainfall-runoff modelling is a cross cutting topic over several of the major issue this subject is focusing on.

Rainfall-runoff models describe a portion of the water cycle (see Figure 1) and therefore the movement of a fluid - water - and therefore they are explicitly or implicitly based on the laws of physics, and in particular on the principles of conservation of mass, conservation of energy and conservation of momentum. Depending on their complexity, models can also simulate the dynamics of water quality, ecosystems, and other dynamical systems related to water, therefore embedding laws of chemistry, ecology, social sciences and so forth.

Models are built by constitutive equations, namely, mathematical formulations of the above laws, whose number depends on the number of variables to be simulated. The latter are the output variables, and the state variables, which one may need to introduce to describe the state of the system. Constitutive equations may include parameters: they are numeric factors in the model equations that can assume different values therefore making the model flexible. In order to apply the model, parameters needs to be estimated (or calibrated, or optimized, and we say that the model is calibrated, parameterized, optimized). Parameters usually assume fixed value, but in some models they may depend on time, or the state of the system.

Figure 1. The hydrologic cycle (from Wikipedia)

Rainfall-runoff models can be classified within several different categories. They can distinguished between event-based and continuous-simulation models, black-box versus conceptual versus process based (or physically based) models, lumped versus distributed models, and several others. It is important to note that the above classifications are not rigid - sometimes a model cannot be unequivocally assigned to one category. We will treat rainfall-runoff models by taking into consideration models of increasing complexity.

1.1. The rational formula

The rational formula is discussed here. We just note in this page that it delivers an estimate of the peak river flow only, and that it is implicitly based on the principle of mass conservation. There is also the implicit use of the principle of conservation of energy in the estimation of the time of concentration. In fact, even if such time is estimated empirically, the mass transfer to the catchment outlet is actually governed by transformation and conservation of energy. Being focused on the estimation of peak river flow, the rational formula is not generally used to assess water resources availability.

1.2. The linear reservoir

The linear reservoir is one of the most used rainfall-runoff models, together with the time area method which we will not discuss here. The linear reservoir model assimilates the catchment to a reservoir, for which the conservation of mass applies. The reservoir is fed by rainfall, and releases the river flow through a bottom discharge, for which a linear dependence applies between the river flow and the volume of water stored in the reservoir, while other losses - including evapotranspiration - are neglected. Therefore, the model is constituted by the following relationships:


where W(t) is the volume of water stored in the catchment at time t, p(t) is rainfall, q(t) is the river flow at time t and k is a constant parameter with the dimension of time (if the parameter was not constant the model would not be linear).

The second equation above assumes a linear relationship between discharge and storage into the catchment. The properties of a linear function are described here. Actually, the relationship between storage in a real tank and bottom discharge is an energy conservation equation that is not linear; in fact, it is given by the well-known Torricelli's law. Therefore the linearity assumption is an approximation, which is equivalent to assuming that the superposition principle applies to runoff generation. Actually, such assumption does not hold in practice, as the catchment response induced by two subsequent rainfall events cannot be considered equivalent to the sum of the individual catchment responses to each single event. However, linearity is a convenient assumption to make the model simpler and analytical integration possible.

A nice feature of the linear reservoir is that the above equations can be integrated analytically, under simplifying assumption. In fact, by substituting the second equation into the first one gets:

Then, by multiplying both sides by et/k and dividing by k one gets:

which can be written as:

By integrating between 0 and t one obtains:

and, by assuming q(0)=0, one finally gets:

which can be easily integrated numerically by using the Euler method.

If one assumes p(t)=constant, an explicit expression is readily obtained for the river flow:

The parameter k is usually calibrated by matching observed and simulated river flows. Its value significantly impacts the catchment response. A high k value implies a large storage into the catchment. Therefore, large values of k are appropriate for catchment with a significant storage capacity. Conversely, a low k value is appropriate for impervious basins. k is also related to the response time of the catchment. A low k implies a quick response, while slowly responding basins are characterised by a large k. In fact, k is related to the response time of the catchment. The linear reservoir is also described here.

1.3. The non-linear reservoir

Several variants of the linear reservoir modeling scheme can be introduced, for instance by adopting a non linear relationship between discharge and storage. Moreover, an upper limit can be fixed for the storage in the catchment, and additional discharges can be introduced, which can be activated for different levels of storage. All of the above modifications make the model non-linear so that an analytical integration is generally not possible. One should also take into account that increasing the number of parameters implies a corresponding increase of estimation variance and therefore simulation uncertainty. Furthermore, it is often observed that introducing additional discharges or thresholds may induce discontinuities in the hydrograph shape. The non-linear reservoir is described here.

Figure 2. A non-linear reservoir (from Wikipedia)

1.4. The Hymod model

The Hymod model is a flexible solution that is increasingly adopted for its capability of providing a good fit in several practical applications. It was originally proposed by Boyle (2000). It is based on the assumption that each point location i in the basin is characterised by a local value of soil water storage Ci, which varies from 0 in the impervious areas up to a maximum value Cmax in the most permeable location of the catchment. Ci is assumed to be randomly varying, so that for an assigned value C* of soil water storage a probability distribution is introduced that gives the probability that a randomly selected location j is characterised by Cj less then, or equal to, C*. Such probability may be interpreted as the fraction F(C*) of the catchment area where CjC*. The above probability distribution is written as

Here, βk is a parameter which quantifies the variability of the soil water storage over the catchment. One can easily verify by numerical simulation that βk=0 implies that the soil water storage is constant over the basin and equal to Cmax; βk=1 implies that the soil water storage is linearly varying from 0 to Cmax; βk→∞ implies that the soil water storage is tending to the null value over the whole catchment, which is therefore impervious.

Let us assume that a storm event occurs over the basin and let us define with the symbol C(t) the time varying water depth stored in the unsaturated locations of the catchment. If we ignore any water losses, like evapotranspiration, C(t) is equal to the rainfall amount from the beginning of the event, which we will indicate with the symbol P(t). If one assumes that the shape of the above probability distribution, now expressed in terms of C(t), is the one reported in Figure 3, it can be easily proved that the water volume stored in the catchment at time t is given by

Figure 3. Distribution of soil water storage, surface runoff and stored volume in Hymod

In fact, the integral at the right hand side of the above equation is the area below the red line in Figure 3. Elementary increment of that area are given by the product of the rainfall at each time step by the fraction F(c(t)) of saturated area at the same time, namely, the product of rainfall by saturated area, which is indeed the surface runoff. Conversely, the area above the curve gives the global storage into the catchment W(t), which can be interpreted as a weighted average value of C(t). The progress of surface runoff and water storage is depicted by the animated picture in Figure 4. Note that after saturation, the contribution of surface runoff is given by the product of the rainfall itself by 1, which is the fraction of saturated area, taking and keeping unit value when the catchment is saturated. After saturation, the storage in the catchment of course does not increase any more.

Figure 4. Progress of surface runoff and stored volume in Hymod

By computing the integral in the latter relationship one obtains:

Therefore the water volume stored in the catchment is equal to

By inverting the above equation one gets

Note that one obtains an estimate of the upper value of the water volume that the catchment can store by imposing C(t)=Cmax, therefore obtaining

The above equations allow an easy application of the Hymod model through a numerical simulation, that is usually carried out by adopting a time step Dt that is equal to observational time step of rainfall and river flow. At a given time step t, one knows the value of C(t) which is equal to the cumulative rainfall depth from the beginning of the event at time t. Therefore, W(t) can be easily computed as well by using the above relationships. At the time t+1, C(t+1)=C(t)+P(t), where P is rainfall, under the condition that C(t+1)=Cmax if C(t)+P(t)>Cmax. Therefore, one can compute a first contribution to surface runoff through the relationship

Finally, one can compute a second contribution to the surface runoff which is given by the water volume that cannot be absorbed by the catchment because part of the catchment area got saturated in the last time step. Such second contribution is given by


By summing ER1+ER2 one obtains the total contribution to the surface runoff during the time step from t to t+1.

At this stage, one may evaluate the water losses given by evapotranspiration at the current time step, and subtract them from W(t+1). Such water losses are computed within Hymod through the relationship

Here, Ep(t) is the potential evapotranspiration at time t. Then, the water storage at time t+1 is given by

One should note that the evapotranspiration is subtracted from the stored water volume after ER1 and ER2 are computed.

The total contribution ER(t)=ER1(t)+ER2(t) to the surface runoff is then divided into 2 components: αER(t) which represent the fast runoff and (1-α)ER(t) which is the slow runoff. αER(t) is propagated through a series of linear reservoirs with the same bottom discharge time constant kq, while (1-α)ER(t) is instead propagated through a single linear reservoir with parameter ks.

The computation moves forward through the sequence of time steps. Hymod counts 5 parameters, namely: Cmax, β, α, kq and ks. These parameters need to be calibrated by using observed data.

Figure 5. A schematic representation of the Hymod model

2. Groundwater models

Groundwater models are computer models of groundwater flow systems. As for rainfall-runoff models, there are several different variants of groundwater models, but the physical principles are the same as for surface models: the conservation of mass and momentum. Groundwater models can be coupled with chemical and/or ecological models to describe water quality therefore allowing to design conservation and remediation strategies.

Groundwater models are based on the application of the groundwater flow equation. A distinguishing behaviour of groundwater flow is its reduced velocity, which makes the assumption of steady flow a reasonable working hypothesis. In steady state conditions groundwater flow is described by Laplace equation which describes potential flow as we will see below.

The groundwater flow equation is applied to a control volume, which in technical applications is the region of space occupied by an aquifer. Let us assume that:

  • Fluid is incompressible;
  • There are no sources of water mass in the control volume,

and let us denotes with W(x,y,z,t), Qin(x,y,z,t) and Qout(x,y,z,t) the water storage, total inflow into, and total outflow from, respectively, the elementary volume with coordinates x,y,z at time t (see Figure 6). Then, the continuity equation (see Figure 6) states that:

To estimate Qin(x,y,z,t), one may write:

  • Volume per unit time into face x = qx(x,y,z,t)Ax = qx(x,y,z,t) dy dz;
  • Volume per unit time into face y = qy(x,y,z,t)Ay = qy(x,y,z,t) dx dz;
  • Volume per unit time into face z = qz(x,y,z,t)Az = qx(x,y,z,t) dx dy.

where qx(x,y,z,t), qy(x,y,z,t) and qz(x,y,z,t) are volumes per unit time and area of the control volume through the faces x, y and z, respectively.

Therefore the total inflow volume per unit time may be written as:

Qin(,y,z,t) = qx(x,y,z,t) dy dz + qy(x,y,z,t) dx dz + qz(x,y,z,t) dx dy .

The total outflow volume Qout(t) per unit time may be written as:


Therefore, one may write that:


We now define Ss as the volumetric specific storage, which is often simply called specific storage: it is the volume of water that an aquifer releases from storage, per volume of aquifer, per unit decline in hydraulic head h(x,y,z,t) in the aquifer. It has the dimension of 1/length. Therefore, we can write that


Therefore, one obtains:


or, in compact form:


Figure 6. A schematic representation of the control volume for a groundwater model

The above equation has both the flux q(x,y,z,t) and the head h(x,y,z,t) unknown. Therefore, it is necessary to couple the continuity equation with a second equation, which is the momentum equation that is derived from the Navier-Stokes equations. Under the assumption that the flow is slow and viscous, the momentum equation reduces to the Darcy's law, which reads as:

where K is hydraulic conductivity, which has the dimension of a velocity and is here assumed to be an isotropic property. By substituting in the above equation one obtains:


which is the desired groundwater flow equation. In steady-state conditions ∂ h(x,y,z,t)/∂ t = 0 and therefore the above equation reduces to Laplace's equation:

2h(x,y,z) = 0

which describes a potential flow. The latter is a flow field where the velocity field can be expressed as the gradient of a scalar function, namely, the velocity potential, which in the above case is the hydraulic head h(x,y,z).

To solve the groundwater flow equation for the sake of describing the hydraulic head, one needs to:

  • specify the geometry of the control volume;
  • specify the material properties (including hydraulic conductivity);
  • specify the initial conditions;
  • specify the boundary conditions;
  • specify the method of solution.

Several groundwater models have been proposed by the literature which provide the above specifications by following different approaches and simplifications.

2.1. Application of groundwater models for estimating water resources availability

Application of groundwater models is subjected to relevant uncertainty. Besides the usual approximation affecting environmental models, when dealing with groundwater one needs to face uncertainty originated by the limited knowledge of the geometry of the control volume and the physical and hydraulic properties of the porous media. Geological prospections and field measurements are needed in order to parametrize models and define initial and boundary conditions. Notwithstanding the presence of uncertainty, models deliver useful information that may allow one to integrate the assessment provided by inductive methods. However, uncertainty needs to be communicated to end users and policy makers, through a quantitative assessment.

Uncertainty may be quantitatively estimated by applying the theory of probability. Usual approaches are based on including an error term to model performances, which may be a multiplicative or additive factors. The error can be defined in probabilist terms, through model validation, which is essentially based on comparing model performances with observed data. Model calibration and validation is described here.

Another approach to quantify uncertainty is to run several model simulation for different geometric configuration of the control volume, different initial and boundary conditions and different parameters, therefore obtaining a range for the model output, which allows to define a range for the desired design variables.

Uncertainty assessment for hydrological models is still an active research field. Applications in the professional practice are still rare.


Boyle, D. P. (2000), Multicriteria calibration of hydrological models, Ph.D. dissertation, Dep. of Hydrol. and Water Resour., Univ. of Ariz., Tucson.

Download the powerpoint presentation of this lecture

Last modified on March 26, 2017