Overview

RUBEM stands for “Rainfall rUnoff Balance Enhanced Model”. It is a spatially distributed hydrological model that integrates classical hydrological processes, i.e.: the rainfall-runoff processes, evapotranspiration processes and aquifer recharge processes. Its formulation is based on the models SPHY [TERINK2015], WEAP [YATES2005], and WetSpass-M [ABDOLLAHI2017].

The model is written in the Python and uses the PCRaster framework [KARSSENBERG2010]. A graphical user interface, called RUBEM Hydrological, that aims to facilitate configuration, execution and visualization of the results of this model is available as a plug-in for QGIS.

The following figure shows the main relevant hydrological processes present in the model.

Model's inputs, water balance calculation and outputs.

Model’s a) inputs, b) water balance calculation and c) outputs. In b) all the hydrological processes are written in black: rainfall, interception, evapotranspiration, runoff, lateral flow, baseflow and recharge. In italic and coloured are the 4 types of soil used for land cover representation: impervious area fraction (ai), vegetated area fraction (av), bare soil area fraction (as) and open water area fraction (ao).

The following sections summarize the model input data as shown in the figure above.

Main Features

The model was developed based on classical concepts of hydrological processes and equations based mainly on the formulations of SPHY models [TERINK2015], WEAP [YATES2005], and WetSpass-M [ABDOLLAHI2017]. The main features of the developed model are:

  • Distributed monthly step model;

  • Hydrological process based on soil water balance in each pixel, and flow total calculated after composition of the resulting accumulated flow, according to Direction drainage network flow established by the digital elevation model (DEM);

  • Calculations for two zones: rootzone and saturated;

  • Evapotranspiration and interception process based on vegetation index: Leaf Area Index (LAI), Photosynthetically Active Radiation Fraction (FPAR) and Normalized Difference Vegetation Index (NDVI); and

  • Sub-pixel level coverage classification, represented by four fractions that represent percentage of total pixel area covered exclusively by: area vegetated, bare soil area, water area and impervious area.

Model Formulation

RUBEM is a spatially distributed hydrological model that integrates classical rainfall-runoff processes. The water balance in the soil layer considers the interception, evapotranspiration, lateral flows, and aquifer recharge. The model integrates classical rainfall-runoff processes, i.e.: the interception evapotranspiration, surface runoff, lateral flows, baseflow, groundwater recharge, and water balance in the soil. Input data are the Normalized Difference Vegetation Index (NDVI), Land Use Land Cover (LULC), weather, rainfall, and soil, including the digital elevation model (DEM).

Interception

Interception is the fraction of precipitation which is retained by the vegetated area canopy. Therefore, its calculation is based on the vegetation cover characteristics, incorporated into the equations through the indexes Fraction of Absorbed Photosynthetically Active Radiation (FPAR), Normalized Difference Vegetation Index (NDVI), and Leaf Area (LAI). Equation (1) shows the interception calculation.

\[I = \alpha_V \cdot I_V\]

where:

  • \(I\) - Interception (mm);

  • \(I_V\)- Interception at the vegetated area (mm);

  • \(\alpha_V\)- Vegetated area fraction (%).

\(I_V\) is calculated through the \(I_R\) (2) and \(I_D\) (4) indexes, based on FPAR (6), NDVI and LAI. The values of LAI (5) are estimated through the FPAR value, which are in turn calculated from NDVI values, obtained through satellite images, according to [PENG2012], and [ABDOLLAHI2017].

\[I_V = P_m \cdot I_R\]
\[I_R = 1 - \exp{\left( -\frac{I_D \cdot d_p}{P_m} \right)}\]
\[I_D = \alpha \cdot LAI \cdot \exp{\left( 1 - \frac{1}{1 + \frac{P_m \cdot [1 - \exp{(-0.463 \cdot LAI)}]}{\alpha \cdot LAI}} \right)}\]
\[LAI = LAI_{max} \cdot \frac{\log{(1-FPAR)}}{\log{(1-FPAR_{max})}} \]
\[FPAR = \min{\left[ \frac{(SR-SR_{min})(FPAR_{max}-FPAR_{min})}{SR_{max}-SR_{min}} + FPAR_{min}, \: 0.95\right]} \]
\[SR = \frac{1 + NDVI}{1 - NDVI} \]

where:

  • \(I_R\)- Interception rate, dependant on land cover characteristics, represented by the Leaf Area Index (LAI) (mm);

  • \(P_m\) - Total monthly precipitation (mm);

  • \(d_P\) - Number of rainy days in the month (days);

  • \(I_D\) - Minimum threshold for daily interception depends on the canopy storage capacity. Its calculation is associated with the LAI (mm);

  • \(LAI\) - Leaf Area Index (-);

  • \(\alpha\) - Interception calibration parameter.

  • \(FPAR\) - Fraction Photosynthetically Active Radiation (-);

  • \(FPAR_{min}\), \(FPAR_{max}\) - minimum and maximum values for FPAR (0.001 and 0.95, respectively), corresponding to the minimum and maximum values for LAI for a particular vegetation class.

Note

If \(P_m = 0\) then \(I_R = 0\).

Evapotranspiration

Evapotranspiration refers to the transfer of water from the soil-plant system to the atmosphere. Actual evapotranspiration is calculated based on the sum of evapotranspiration values for each sub-pixel level coverage classification fraction (vegetation, bare soil, water, and impervious soil). Each evapotranspiration fraction is estimated by the concept of potential evapotranspiration [ABDOLLAHI2017]. The potential evapotranspiration values are an input calculated in this study using the Penman-Montheit method [ALLEN2018]. The following equations presents the calculation process.

\[ET_{REAL} = \alpha_V \cdot ET_{R,V} + \alpha_{SN} \cdot ET_{R,S} + \alpha_A \cdot ET_{R,A} + \alpha_I \cdot ET_{R,I}\]

where:

  • \(ET_{REAL}\) – Total real evapotranspiration (mm);

  • \(ET_{R,V}\) – Real evapotranspiration at the vegetated area (mm);

  • \(ET_{R,S}\) – Real evapotranspiration at the bare soil area (mm);

  • \(ET_{R,A}\) – Real evapotranspiration at the water area (mm);

  • \(ET_{R,I}\) – Real evapotranspiration at the impervious area (mm);

  • \(\alpha_V\) – Vegetated area fraction (%);

  • \(\alpha_S\) – Bare soil area fraction (%);

  • \(\alpha_A\) – Water area fraction (%);

  • \(\alpha_I\) – Impervious area fraction (%).

Note

If \(\alpha_A = 1\) and \(ET_{R,A} > P_m\) then \(ET_{R,A} = P_m\).

Vegetated Area Fraction

\[ET_{R,V} = ET_p \cdot kc \cdot ks\]
\[kc = kc_{min} + (kc_{max} - kc_{min}) \cdot \left( \frac{NDVI-NDVI_{min}}{NDVI_{max}-NDVI_{min}} \right)\]

where:

  • \(ET_p\) – Potential evapotranspiration (mm);

  • \(kc\) – Crop coefficient (-);

  • \(ks\) – Soil moisture reduction coefficient (-).

  • \(NDVI\) - Normalized Difference Vegetation Index (-).

Note

If \(NDVI \leq 1.1 \cdot NDVI_{min}\) then \(kc = kc_{min}\).

Bare Soil Area Fraction

\[ET_{R,S} = ET_p \cdot kc_{min} \cdot ks\]
\[ks = \frac{\ln{(TU_R - TU_{PM} + 1)}}{\ln{(TU_{CC} - TU_{PM} + 1)}} \]

where:

  • \(ET_p\) – Potential evapotranspiration (mm);

  • \(kc\) – Crop coefficient (-);

  • \(ks\) – Soil moisture reduction coefficient (-);

  • \(TU_{PM}\) – Moisture content of the soil at the wilting point (mm);

  • \(TU_{CC}\) – Moisture content of the soil at field capacity (mm);

  • \(TU_R\) – Actual moisture content of the soil (mm).

Note

If \(TU_R < TU_{PM}\) then \(ks = 0\).

Water Area Fraction

\[ET_{R,A} = \frac{ET_p}{kp}\]
\[kp = 0.482 + 0.024 \cdot \ln{(B)} - 0.000376 \cdot U_2 + 0.0045 \cdot UR \]

where:

  • \(ET_p\) – Potential evapotranspiration (mm);

  • \(kp\) – Pan Coefficient (-).

Impervious Area Fraction

\[ET_{R,I} = I_I\]

where:

  • \(I_I\) – Interception for impervious areas (3 -5 mm);

  • \(ET_{R,I}\) – Real evapotranspiration at the impervious area (mm).

Surface Runoff

Surface runoff refers to the water that flows on the cell surface. The monthly calculation of surface runoff is based on the rational method, considering two flow coefficients, as shown in equation (16) [ABDOLLAHI2017].

\[SR = C_{SR} \cdot C_h \cdot (P_m - I) \]

where:

  • \(SR\) – Surface runoff (mm);

  • \(C_{SR}\) – Real runoff coefficient (-);

  • \(C_h\) – Coefficient related to soil moisture (-);

  • \(P_m\) – Total monthly precipitation (mm);

  • \(I\)– Total interception (mm).

\(C_{SR}\) is a real flow coefficient, obtained through a flow coefficient which considers pervious and impervious areas and adjusted by the mean monthly precipitation. \(C_h\) is related to soil moisture conditions. When the soil is saturated at the current time step, surface runoff is calculated as the difference between precipitation and interception.

\[C_{SR} = \frac{C_{wp} \cdot \overline{P}_{MD}}{C_{wp} \cdot \overline{P}_{24} - RCD \times C_{wp} + RCD} \]
\[C_{wp} = (1 - A_{imp}) \cdot C_{per} + A_{imp} \cdot C_{imp} \]
\[C_{imp} = 0.09 \cdot \exp{(2.4 \cdot A_{imp})} \]
\[A_{imp} = \alpha_O + \alpha_I \]
\[C_{per} = w_1 \left(\frac{0.02}{n} \right) + w_2\left(\frac{\theta_{PM}}{1-\theta_{PM}} \right) + w_3\left(\frac{S}{10+S} \right) \]
\[C_h = \left(\frac{\theta_{TUR}}{\theta_{POR}} \right)^b \]
\[\theta_{TUR} = \frac{TU_R}{d_g \cdot Z_r \cdot 10}\]

where:

  • \(SR\) – surface runoff in the cell (mm);

  • \(C_{SR}\) – actual flow coefficient (-);

  • \(C_h\) – coefficient related to soil moisture (-);

  • \(P_m\) – total monthly precipitation (mm);

  • \(I\) – total interception total (mm);

  • \(\overline{P}_{MD}\) – average daily rain on rainy days (mm day-1month-1);

  • \(RCD\) – regional consecutive dryness factor [parameter to be calibrated] (mm);

  • \(C_{wp}\) – potential runoff coefficient (-);

  • \(C_{imp}\) – potential runoff coefficient of impermeable areas, empirical relationship (-);

  • \(C_{per}\) – potential runoff coefficient of permeable areas (-);

  • \(A_{imp}\) – total fraction of impermeable area per cell (-);

  • \(w1\), \(w2\) e \(w3\) – Weights for the three components related to runoff: cover, soil and slope characteristics [parameter to be calibrated] (-);

  • \(n\) – Manning roughness (-);

  • \(\theta_{PM}\) – root layer wilting point volumetric content (-);

  • \(S\) – terrain slope in each cell (%);

  • \(\theta_{TUR}\) – volumetric moisture content of the root layer (-);

  • \(\theta_{POR}\) – volumetric porosity content of the root layer (-);

  • \(b\) – calibration exponent, ranges between 0 and 1 [parameter to be calibrated] (-);

  • \(TU_R\) – moisture content for the root zone (mm);

  • \(d_g\) – overall soil density in the root layer (g cm3);

  • \(Z_r\) – root layer thickness (cm).

Note

If \(\theta_{TUR} > \theta_{POR}\) then \(\theta_{TUR} = \theta_{POR}\).

Note

If \(\alpha_A = 1\) then \(SR = P_m - ET_{R,A}\).

Lateral Flow

Lateral flow occurs in the unsaturated layer under the surface of the soil. The groundwater recharge is the downflow from the surface to the water aquifer by infiltration and deep percolation. Lateral flow and groundwater recharge (equations below) are calculated by values of root zone moisture, hydraulic conductivity, and a calibrated partitioning coefficient, which controls horizontal and vertical flow [YATES2005].

\[LF = f \cdot K_R \cdot \left( \frac{TU_R}{TU_{SAT}} \right)^2 \]

where:

  • \(LF\) – Lateral flow (mm);

  • \(f\) – Vertical and horizontal flow partitioning coefficient (-);

  • \(K_R\) – Root zone hydraulic conductivity coefficient (mm/month);

  • \(TU_R\) – Root zone moisture content (mm);

  • \(TU_{SAT}\) – Root zone moisture content at saturation (mm).

Aquifer Recharge

Groundwater recharge is complementary to lateral flow, and is calculated as follows (25):

\[REC = (1-f) \cdot K_R \cdot \left( \frac{TU_R}{TU_{SAT}} \right)^2 \]

where:

  • \(REC\) – Recharge (mm);

  • \(f\) – Vertical and horizontal flow partitioning coefficient (-);

  • \(K_R\) – Root zone hydraulic conductivity coefficient (mm/month);

  • \(TU_R\) – Root zone moisture content (mm);

  • \(TU_{SAT}\) – Root zone moisture content at saturation (mm).

Baseflow

The basic flow is the one that occurs in the saturated soil layer. The base flow is determined by [TERINK2015] by groundwater recharge and the recession coefficient. It is calculated only if the moisture content in the saturated zone exceeds a specified threshold.

\[ BF = \left\{ \begin{array}{ll} 0, & \mbox{if } TU_S \leq BF_{thresh} \\ BF_{T-1} \cdot e^{-\alpha_{GW}} + (1-e^{-\alpha_{GW}}) \cdot REC, & \mbox{if } TU_S > BF_{thresh} \\ \end{array} \right. \]

where:

  • \(BF\) – Baseflow (mm);

  • \(BF_{T-1}\) – Baseflow at the previous time step (mm);

  • \(\alpha_{GW}\) – Baseflow decay coefficient (-) [parameter to be calibrated];

  • \(REC\) – Recharge (mm);

  • \(TU_S\) – Saturated zone moisture content (mm);

  • \(BF_{thresh}\) – Threshold baseflow, attributed for each watershed (mm).

Water Balance

RUBEM is governed by the mass balance for the water in the soil layer. Soil moisture at the root and saturated zones are calculated through a water balance equation on each respective layer according to equations (27), (28) and (29).

\[P_E = P_m - I\]

where:

  • \(P_E\) – Effective precipitation (mm);

  • \(P_m\) – Total Monthly precipitation (mm);

  • \(I\) – Total interception (mm).

\[TU_S = TU_{S,T-1} - BF + REC \]

where:

  • \(TU_S\) – Saturated zone moisture content (mm);

  • \(TU_{S,T-1}\) – Saturated zone moisture content at the previous time step (mm);

  • \(BF\) – Baseflow (mm).

\[TU_R = TU_{R,T-1} + P_E - SR - LF - REC - ET_{REAL} \]

where:

  • \(TU_R\) – Root zone moisture content (mm);

  • \(TU_{R,T-1}\) – Root zone moisture content at the previous time step (mm);

  • \(SR\) – Surface runoff (mm);

  • \(LF\) – Lateral flow (mm);

  • \(REC\) – Recharge (mm);

  • \(ET_{REAL}\) – Total real evapotranspiration (mm).

Note

If \(\alpha_A = 1\) then \(TU_R = TU_{SAT}\).

In the root zone, the surface runoff relies on the soil moisture conditions, based on the lateral flow, recharge, and evapotranspiration. The saturated zone of the soil affected the calculation of the base flow and recharge parameters. The water balance equation adopted in RUBEM aims to calculate the total superficial flow.

Total Discharge

Total Discharge corresponds to the contribution of the flows produced in each cell. It is computed at the grid level of the model, and the total flow is computed using the Digital Elevation Model to generate the Local Drainage Direction (LDD). The superficial flow is computed at any grid cell, Equation (31), and the total flow is computed using an aggregate result according to the LDD [KARSSENBERG2010]; [TERINK2015].

Total superficial discharge, Equation (30), is the sum of flow components: surface runoff (SR), lateral flow (LF), and baseflow (BF) [ABDOLLAHI2017].

\[Q_{Tot} = SR + LF + BF \]
\[Q_t = x \cdot Q_{t-1} + \frac{0.001 \cdot (1 - x) \cdot A \cdot Q_{Tot}}{days \cdot 24 \cdot 3600} \]

where:

  • \(Q_{Tot}\) – Total pixel runoff (mm);

  • \(SR\) – Surface runoff (mm);

  • \(LF\) – Lateral flow (mm);

  • \(BF\) – Baseflow (mm);

  • \(Q_t\) – Total pixel flow (m3s-1);

  • \(Q_{t-1}\) – Total pixel flow at the previous time step (m3s-1);

  • \(A\) – Pixel area (m2);

  • \(x\) – Damping coefficient (-).

The model operates on a monthly basis and calculates the water balance cell-by-cell. Local Drainage Direction (LDD) is then calculated based on the Digital Elevation Model (DEM), which is used for obtaining the total discharge.

Calibration and Validation

RUBEM calibration is based on 9 parameters as described in the table below. Their range is defined according to the watershed physical features (soil and rainfall) and literature review [TERINK2015]; [ABDOLLAHI2017] and [WANG2018].

Parameter

Description

Restriction

Interception

Parameter (\(\alpha\))

Interception parameter value. Represents the

daily interception threshold that depends on

land use.

\(0.01 ≤ \alpha ≤ 10\)

Rainfall Intensity

Coefficient (\(b\))

Exponent value representing the effect of

rainfall intensity in the runoff.

\(0.01 ≤ b ≤ 1\)

Land Use Factor

Weight (\(w_1\))

Weight of land use factor. It measures the

effect of the land use in the potential

runoff produced.

\(w_1 + w_2 + w_3 = 1\)

Soil Factor

Weight (\(w_2\))

Weight of soil factor. It measures the

effect of the soil classes in the potential

runoff produced.

\(w_1 + w_2 + w_3 = 1\)

Slope Factor

Weight (\(w_3\))

Weight of slope factor. It measures the

effect of the slope in the potential

runoff produced.

\(w_1 + w_2 + w_3 = 1\)

Regional Consecutive

Dryness Level (\(RCD\))

Regional Consecutive Dryness level,

incorporates the intensity of rain and the

number of consecutive days in runoff

calculation.

\(1 ≤ RCD ≤ 10\)

Flow Direction

Factor (\(f\))

Used to partition the flow out of the root

zone between interflow and flow to the

saturated zone.

\(0.01 ≤ f ≤ 1\)

Baseflow Recession

Coefficient (\(\alpha_{GW}\))

Decimal value referring to the recession

coefficient of the baseflow.

\(0.01 ≤ \alpha_{GW} ≤ 1\)

Flow Recession

Coefficient (\(x\))

Flow recession coefficient value, it

incorporates a flow delay in the accumulated

amount of water that flows out of the cell

into its neighboring downstream cell.

\(0 ≤ x ≤ 1\)

For calibration purposes, the Differential Evolution algorithm [STORN1997] can be used. The method executes three operations for candidate choice, mutation, crossing and selection. For each basin, spatialized distributed stream gauges have to be considered. The choice criteria must include series with few flaws, location in the basin, LULC, and rainfall regimes.

Traditional indexes and statistics can be used to evaluate the accuracy of the results. Root Mean Square Error – RMSE and Nash-Sutcliffe Efficiency – NSE are examples.

References

[ABDOLLAHI2017] (1,2,3,4,5,6,7)

Abdollahi, K., Bashir, I., Harouna, M., Griensven, A., Huysmans, M., Batelaan, O., Verbeiren, B., A distributed monthly water balance model: formulation and application on Black Volta Basin, Environ Earth Sci, 76:198, 2017.

[ALLEN2018]

Allen, M., Dube, O. P., Solecki, W., Aragón-Durand, F., Cramer, W., Humphreys, S., … & Mulugetta, Y. (2018). Global warming of 1.5° C. An IPCC Special Report on the impacts of global warming of 1.5° C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty. Sustainable Development, and Efforts to Eradicate Poverty.

[KARSSENBERG2010] (1,2)

Karssenberg, D., Schmitz, O., Salamon, P., de Jong, K., and Bierkens, M. F. P.: A software framework for construction of process-based stochastic spatio-temporal models and data assimilation, Environmental Modelling & Software 25(4), 489-502, 2010. doi: 10.1016/j.envsoft.2009.10.004.

[MACHADO2017]

Machado, A.R. (2017). Alternativas de restauração de florestas ripárias para o fornecimento de serviços ecossistêmicos. 149 p. Tese (Doutorado em Engenharia de Recursos Hídricos) – Escola politécnica, Universidade de São Paulo, São Paulo, 2017.

[MCCUEN2006]

McCuen, R.H.; Knight, Z; Cutter, A.G. (2006). “Evaluation of the Nash–Sutcliffe efficiency index”. Journal of Hydrologic Engineering. v. 11, n. 6, 2006.

[OLIVOS2017]

Olivos, L.M. (2017). Sustentabilidade do uso de recursos hídricos superficiais e subterrâneos no município de São Carlos, SP. 2017. Dissertação (Mestrado em Recursos Hídricos) – Escola Politécnica, Universidade de São Paulo, São Paulo, Brazil.

[PENG2012]

Peng, D.; Zhang, B.; Liu, L. Comparing spatiotemporal patterns in Eurasian FPAR derived from two NDVI based methods. International Journal of Digital Earth, v. 5, n. 4, p. 283-298, 2012

[RITTER2013]

Ritter, A.; Muñoz-Carpena, R. (2013). “Performance evaluation of hydrological models: statistical significance for reducing subjectivity in goodness-of-fit assessments”. Journal of Hydrology. 480 (1): 33–45

[SINGH2004]

Singh, J., Knapp, H.V., Demissie, M. (2004). Hydrologic modeling of the Iroquois River watershed using HSPF and SWAT. 2004. Illinois State Water Survey Contract Report 2004-08. Illinois Department of Natural Resources: Illinois State Geological Survey.

[STORN1997]

Storn, R., Price, K. (1997). Differential Evolution – A simple and efficient heuristic for global optimization over Continuous spaces. Journal of Global Optimization, v. 11 p. 341–359.

[TENA2019]

Tena, T.M.; Mwaanga, P.; Nguvulu, A. (2019). Hydrological modelling and water resources assessment of Chongwe River Catchment using WEAP model. Water, v. 11, p. 839-856

[TERINK2015] (1,2,3,4,5)

Terink, W., Lutz, A. F., Simons, G. W. H., Immerzeel, W. W., & Droogers, P. (2015). SPHY v2.0: Spatial Processes in HYdrology. Geoscientific Model Development, 8(7), 2009–2034. https://doi.org/10.5194/gmd-8-2009-2015

[YATES2005] (1,2,3)

Yates, D., Sieber, J., Purkey, D., & Huber-Lee, A. (2005). WEAP21—A Demand-, Priority-, and Preference-Driven Water Planning Model. Water International, 30(4), 487–500. https://doi.org/10.1080/02508060508691893

[WANG2018]

Wang, L., Wang, Z., Liu, C., Bai, P., & Liu, X. (2018). A Flexible Framework HydroInformatic Modeling System—HIMS. Water, 10(7), 962. https://doi.org/10.3390/w10070962