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.
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.
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].
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.
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
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
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
where:
\(ET_p\) – Potential evapotranspiration (mm);
\(kp\) – Pan Coefficient (-).
Impervious Area Fraction
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].
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.
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].
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):
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.
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).
where:
\(P_E\) – Effective precipitation (mm);
\(P_m\) – Total Monthly precipitation (mm);
\(I\) – Total interception (mm).
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).
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].
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
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.
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.
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.
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.
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.
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.
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
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
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.
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.
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
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
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
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