# A Physically Based Runoff Model Analysis of the Querétaro River Basin

Authors: Carlos Javier Villa Alvarado; Eladio Delgadillo-Ruiz; Carlos Alberto Mastachi-Loza; Enrique González-Sosa; Ramos Salinas Norma Maricela
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Today the knowledge of physical parameters of a basin is essential to know adequately the rainfall-runoff process; it is well known that the specific characteristics of each basin such as temperature, geographical location, and elevation above sea level affect the maximum discharge and the basin time response. In this paper a physically based model has been applied, to analyze water balance by evaluating the volume rainfall-runoff using SHETRAN and hydrometric data measurements in 2003. The results have been compared with five ETp different methodologies in the Querétaro river basin in central Mexico. With these results the main effort of the authorities should be directed to better control of land-use changes and to working permanently in the analysis of the related parameters, which will have a similar behavior to changes currently being introduced and presented in observed values in this basin. This methodology can be a strong base for sustainable water management in a basin, the prognosis and effect of land-use changes, and availability of water and also can be used to determine application of known basin parameters, basically depending on land-use, land-use changes, and climatological database to determine the water balance in a basin.


---

## Body

## 1. Introduction

The main problems in the world are associated with the quantity and quality of existing data [1], which represent more time of process and uncertainty in the obtained results. Moreover the government policies are limited in making growing plans for urban centers and do not consider a scheme of penalties when variables as land-use are not respected because there is no knowledge of the effect of changes in the basin’s hydrology. As mentioned in Amini et al. [2], in which a dendritic watershed system is used to determine the conditions of a watershed affected by land-use changes, they do not consider the full quantity of precipitation in one observed year or the effect of antecedent soil moisture.In Mitsuda and Ito [3], different factors that affect the response of the basin are listed such as the land-use changes, the landowners activities, and the poor political strategies that increase the climate change effects [4], biodiversity alteration, and changes in biogeochemical and hydrological cycle; however they assume that natural and environmental factors are based only in the slope that is one of many variables that should be studied in a hydrologic basin, where the main interest is related to environmental topics as mentioned in Su and Christensen, 2013 [5].In Mendoza et al. [6], the implications for land-use and land cover changes are soil erosion and degradation that disturbs the ecosystem capability to provide the natural benefits such as superficial and underground water retention; nevertheless it is considered a distributed model that considers the runoff as a function of biophysical characteristics and not based on the physical basin variables. With the product of this research, as mentioned in Li and Huang, 2012 [7], a linkage can be applied among environmental regulations and economic implications with the use of informatics tools; finally through numerical modeling and applied research an environmental linkage can be carried out that involves the law and the economic implications [7].In many places there are problems associated with the quantity and quality of existing dataset, representing time and uncertainty of the data. Moreover, government policies are limited in making decisions growing plans for urban centers and do not present a scheme of penalties when the laws on land-use changes are not respected.As mentioned in Ewen [8], the physically based models are used to determine the water flow in the surface, underground, and environment. Surface and underground processes are widely explained in Birkinshaw and Bathurst [9] wherein SHETRAN is used as a distributed hydrologic model.The importance of this model in water balance is related to the places in the world where quantity and quality of measurement stations are poor with lack of data in most of cases. The insufficient quantity of climate datasets is the main problem to solve to make decisions for annual water use management in this basin and its use.The results suggest that the use of SHETRAN in water balance in the Querétaro river basin has the advantage that it depends on the knowledge of physic parameters of the basin and can be applied in unmeasured basins when the stations are damaged or do not have resources to implement a system of hydrometric monitoring.


## 2. Material and Methods

### 2.1. Location

The Querétaro river basin area is 2,707.74 square kilometers, located in central Mexico, 200 km to the northwest of Mexico City (Figure1). Its coordinates are UTM-WGS84_14N x
(
E
) = 358,380; y
(
N
) = 2′282,700; 1,800 meters above sea level (masl). The basin drains to the Ameche hydrometric station (1,787 masl), 20 kilometers west of Querétaro city. The annual rainfall rate is 55 centimeters per year, mean annual temperature is 18.7°C, and mean annual potential evaporation is 2,050–2,200 mm. This basin is close to the Continental Divide of the Americas and drains to Lake Chapala in the state of Jalisco, 450 km to the west of Querétaro city.Figure 1
Central Mexico: Querétaro river basin. In the middle of Guanajuato and Querétaro State. SHETRAN terrain model values were aggregated to a 2,100 m mesh. At north (2,710) and south (2,610) are the major elevations above the sea level.The total rainfall registered in climate stations in 2003 (Table1) represents twice the average annual volume. As mentioned in Nayak and Mandal [13] and Mahmood et al., [14], these effects can be attributed to the phenomenon of global warming so that the high rainfall events should be analyzed in shorter time periods to determine the response of basin to these changes [15].Table 1
Climate stations: the 2003 total rainfall measured is shown. The maximum difference found in altitude is 85 m, where the minimum is Carrillo (1,800 masl) and the highest is Juriquilla (1,885 masl).
Name
ID
Locality
Municipality
Latitudea
Longitudea
Altitude (masl)
Total 2003 rainfallb (mm)
El Batan
22004
El Batan
Corregidora
20°30′17.33′′
100°25′27.08′′
1881
1,119.50
El Pueblito
22006
El Pueblito
Corregidora
20°32′17.92′′
100°26′25.55′′
1826
1,079.00
Carrillo
22027
Querétaro
Querétaro
20°36′12.60′′
100°25′55.49′′
1800
994.80
Juriquilla
22045
Juriquilla
Querétaro
20°43′04.99′′
100°27′34.00′′
1885
904.20
Querétaro
22063
Querétaro
Querétaro
20°33′48.42′′
100°22′09.88′′
1813
934.20
Mean
1,006.34c
aGeographical coordinate system (WGS-84. 14 north).
bTotal annual measured values in the climate station.
cAnnual mean precipitation in Querétaro is 550 mm-year−1.


### 2.2. Basin

The last years there has been an increase in urban land-use changes, principally in east, west, and northwest of the basin, due to topographic natural conditions and industrial-commercial activity increasing. Because of this recent increase in land-use changes, the basin is subject to floods more frequently.In this basin there is no methodology to determine the water balance based on the characteristics of the region, so that different formulas are applied based on existing regulations [16], which can be implemented in a better way if the physical parameters of the basin are considered, due to the poor level of hydrometric data quality; nevertheless actually these data are used to determine the water balance. A physically based method can be applied calculating the runoff by using the precipitation and land-use conditions. At the same time the margin of uncertainty that can be very sensitive in low availability of water resources in the basin must be continuously evaluated.


### 2.3. Digital Terrain Model

A digital terrain model (DTM) was obtained from INEGI (Instituto Nacional de Estadística y Geografía, México) topographic maps (INEGI, 2010), at a scale of 1 : 50,000. The data was processed with ArcMap 3D_Analyst to define the catchment [17], using
(1)
Z
(
x
)
=
m
(
x
)
+
,
′
(
x
)
+
,
′
′
,
where Z
(
x
) is the value of the elevation in the position x, m
(
x
) is the function that describes the dataset ,
′
(
x
) is the stochastic local variation in the dataset after the m
(
x
) variation, and ,
′
′ is the no tendencies variability dataset (noise). When the dataset has no tendencies m
(
x
)
=
0, then
(2)
E
[
Z
(
x
)
-
Z
(
x
+
h
)
]
=
0
,
where h is the distance among the evaluated points and the variogram γ
(
h
) as follows:
(3)
Var
⁡
[
Z
(
x
)
-
Z
(
x
+
h
)
]

=

E
[
{
Z
(
x
)
-
Z
(
x
+
h
)
}
2
]

=

2
γ
(
h
)
.
Figure 1 shows the resulting digital terrain model.


### 2.4. Runoff and Rivers

In this basin there are three main rivers (Figures2 and 2(a)): Arenal river (from north to south), Pueblito river (from south to north), and Querétaro river at the center of the basin (from east to west), which all converge close to the Querétaro-Guanajuato state border. All of them flow across the Querétaro city urban area, where in recent years repetitive flooding has been observed. On the one hand the risk of flooding damage to population affects public services, economic activities, and social activities. On the other hand, the water deficit occurs for most of the year and in dry years the total rainfall is less than 450 mm [11].Figure 2
(a) Climate stations and land-use in Querétaro river catchment. Land-use: 65% correspond to arable land (grass, arable land and shrub), urbanized and paved surfaces 14% (mainly at the center of the basin; 1,870 masl) and forest and water 21% (at north and south, in the higher elevations). In (a), Querétaro (32.73 km), Pueblito (16.29 km) and Arenal (11.76 km) main rivers and flow direction are shown.


### 2.5. Satellite Images

A Landsat TM Satellite geotiff30
×
30 meters image (September 2003) was used to generate the grid used in the hydrological model. This image was classified with ArcMap 5_4_3 (red, green and blue) bands a reflectance of R = C
2, G = C
3, B = C
1. The resulting image was aggregated to obtain a 2.1
×
2.1-kilometer grid model. A classification method supervised was used with in situ interpretation and land-use. Satellite data was obtained from the United States Geological Survey (USGS) and was used to reduce the landscape and land-use cover uncertainty [18] in this basin in natural and urban areas. As a result of the satellite image analysis, the forest represents 20%, arable land 65%, water bodies 1%, and urban paved surfaces 14% (Figure 2).


### 2.6. Land-Use and Interception

In semiarid areas of Mexico, matorral intercepts from 1.9% to 7.7% of gross rainfall [19], deciduous trees from 9% to 20%, and coniferous trees from 20% to 48% [12]. In urban areas, rainfall is intercepted by buildings and vegetation. Buildings intercept rainfall on rooftops and laterally on wall surfaces, modifying the hydrological response peak time [20]. Canopy shrub and grass saturation was estimated at 1.5 mm using the linear relationship between throughfall and steamflow, compared to forest saturation, which was observed at 5.0 mm [12].The edge effects present in isolated tree canopies might increase evaporation from wet canopies [21, 22]. Urban vegetation occurs as individuals or parklands, and therefore area-based interception percentages are smaller than reported for natural forests, grass, and possibly agroforests [12, 23–25]. Klaassen et al. [26] affirm that proximity to the forest edge affects both the interception storage capacity and the rate of evaporation of intercepted water, which cancel each other.


### 2.7. Soil

In Querétaro river basin the following different soils exist: litosol, vertisol, fluvisol, phaeozem, chernozem, castañozem, and yermosol (Figure3 and Table 3), according to World Reference Base for Soil Resources of the Food and Agriculture Organization (FAO), where clay, loam, and sand are the principal elements found in the basin. The maximum depth layer is found at 5 m with a rock basement detected, located at the north of the basin. Vertisol is located in the center of the basin, litosol at the south, north, and east, and phaeozem in the upper areas coinciding with the watershed mainly at the north and south; other soils are small area of fluvisol (north), chernozem (south) in Huimilpan municipality, and yermosol and castañozem.Figure 3
SHETRAN soils in the Querétaro river basin (mesh resolution: 2,100 meters). The existing soil composition in this basin isvertisol (52%), phaeozem (23%), litosol (18%), and fluvisol, castañozem, chernozem, and yermosol (7%). Vertisol principally at the center of the basin, atphaeozem north and south, and litosol found at the center, northwest, and southwest.


### 2.8. Precipitation

In this region the main rainfall volume occurs in summer (June to September), but the weather stations usually have breakdowns, resulting in incomplete databases that usually are filled by calculating missing values; to determine the mean precipitation 80 rainfall years have been analyzed. As a result dry years are considered between 257 and 550 mm/year.The maximum rainfall volume occurs between 5,850 and 6,800 hours, but May 2003 (3,700) and June (4,500) showed an increase in volume annual rainfall. The maximum event corresponds to 68 mm (24-hour rainfall register). The 100 cm registered precipitation is considered in this paper to evaluate the basin model calibration and runoff total volume. The spatial rainfall analysis was made with daily precipitation records, using for this 5 climate stations which were visited to evaluate if they comply with the World Meteorological Organization (WMO) principally related to obstructions at surroundings [27].The rainfall applied in the basin was the daily data obtained from the government climate stations in data station: 22027 (Carrillo), 22045 (Juriquilla), El Batan (22004), 22006 (Pueblito), and 22063 (Querétaro). The climate stations details and location are shown in Table1 and Figure 2.


## 2.1. Location

The Querétaro river basin area is 2,707.74 square kilometers, located in central Mexico, 200 km to the northwest of Mexico City (Figure1). Its coordinates are UTM-WGS84_14N x
(
E
) = 358,380; y
(
N
) = 2′282,700; 1,800 meters above sea level (masl). The basin drains to the Ameche hydrometric station (1,787 masl), 20 kilometers west of Querétaro city. The annual rainfall rate is 55 centimeters per year, mean annual temperature is 18.7°C, and mean annual potential evaporation is 2,050–2,200 mm. This basin is close to the Continental Divide of the Americas and drains to Lake Chapala in the state of Jalisco, 450 km to the west of Querétaro city.Figure 1
Central Mexico: Querétaro river basin. In the middle of Guanajuato and Querétaro State. SHETRAN terrain model values were aggregated to a 2,100 m mesh. At north (2,710) and south (2,610) are the major elevations above the sea level.The total rainfall registered in climate stations in 2003 (Table1) represents twice the average annual volume. As mentioned in Nayak and Mandal [13] and Mahmood et al., [14], these effects can be attributed to the phenomenon of global warming so that the high rainfall events should be analyzed in shorter time periods to determine the response of basin to these changes [15].Table 1
Climate stations: the 2003 total rainfall measured is shown. The maximum difference found in altitude is 85 m, where the minimum is Carrillo (1,800 masl) and the highest is Juriquilla (1,885 masl).
Name
ID
Locality
Municipality
Latitudea
Longitudea
Altitude (masl)
Total 2003 rainfallb (mm)
El Batan
22004
El Batan
Corregidora
20°30′17.33′′
100°25′27.08′′
1881
1,119.50
El Pueblito
22006
El Pueblito
Corregidora
20°32′17.92′′
100°26′25.55′′
1826
1,079.00
Carrillo
22027
Querétaro
Querétaro
20°36′12.60′′
100°25′55.49′′
1800
994.80
Juriquilla
22045
Juriquilla
Querétaro
20°43′04.99′′
100°27′34.00′′
1885
904.20
Querétaro
22063
Querétaro
Querétaro
20°33′48.42′′
100°22′09.88′′
1813
934.20
Mean
1,006.34c
aGeographical coordinate system (WGS-84. 14 north).
bTotal annual measured values in the climate station.
cAnnual mean precipitation in Querétaro is 550 mm-year−1.


## 2.2. Basin

The last years there has been an increase in urban land-use changes, principally in east, west, and northwest of the basin, due to topographic natural conditions and industrial-commercial activity increasing. Because of this recent increase in land-use changes, the basin is subject to floods more frequently.In this basin there is no methodology to determine the water balance based on the characteristics of the region, so that different formulas are applied based on existing regulations [16], which can be implemented in a better way if the physical parameters of the basin are considered, due to the poor level of hydrometric data quality; nevertheless actually these data are used to determine the water balance. A physically based method can be applied calculating the runoff by using the precipitation and land-use conditions. At the same time the margin of uncertainty that can be very sensitive in low availability of water resources in the basin must be continuously evaluated.


## 2.3. Digital Terrain Model

A digital terrain model (DTM) was obtained from INEGI (Instituto Nacional de Estadística y Geografía, México) topographic maps (INEGI, 2010), at a scale of 1 : 50,000. The data was processed with ArcMap 3D_Analyst to define the catchment [17], using
(1)
Z
(
x
)
=
m
(
x
)
+
,
′
(
x
)
+
,
′
′
,
where Z
(
x
) is the value of the elevation in the position x, m
(
x
) is the function that describes the dataset ,
′
(
x
) is the stochastic local variation in the dataset after the m
(
x
) variation, and ,
′
′ is the no tendencies variability dataset (noise). When the dataset has no tendencies m
(
x
)
=
0, then
(2)
E
[
Z
(
x
)
-
Z
(
x
+
h
)
]
=
0
,
where h is the distance among the evaluated points and the variogram γ
(
h
) as follows:
(3)
Var
⁡
[
Z
(
x
)
-
Z
(
x
+
h
)
]

=

E
[
{
Z
(
x
)
-
Z
(
x
+
h
)
}
2
]

=

2
γ
(
h
)
.
Figure 1 shows the resulting digital terrain model.


## 2.4. Runoff and Rivers

In this basin there are three main rivers (Figures2 and 2(a)): Arenal river (from north to south), Pueblito river (from south to north), and Querétaro river at the center of the basin (from east to west), which all converge close to the Querétaro-Guanajuato state border. All of them flow across the Querétaro city urban area, where in recent years repetitive flooding has been observed. On the one hand the risk of flooding damage to population affects public services, economic activities, and social activities. On the other hand, the water deficit occurs for most of the year and in dry years the total rainfall is less than 450 mm [11].Figure 2
(a) Climate stations and land-use in Querétaro river catchment. Land-use: 65% correspond to arable land (grass, arable land and shrub), urbanized and paved surfaces 14% (mainly at the center of the basin; 1,870 masl) and forest and water 21% (at north and south, in the higher elevations). In (a), Querétaro (32.73 km), Pueblito (16.29 km) and Arenal (11.76 km) main rivers and flow direction are shown.


## 2.5. Satellite Images

A Landsat TM Satellite geotiff30
×
30 meters image (September 2003) was used to generate the grid used in the hydrological model. This image was classified with ArcMap 5_4_3 (red, green and blue) bands a reflectance of R = C
2, G = C
3, B = C
1. The resulting image was aggregated to obtain a 2.1
×
2.1-kilometer grid model. A classification method supervised was used with in situ interpretation and land-use. Satellite data was obtained from the United States Geological Survey (USGS) and was used to reduce the landscape and land-use cover uncertainty [18] in this basin in natural and urban areas. As a result of the satellite image analysis, the forest represents 20%, arable land 65%, water bodies 1%, and urban paved surfaces 14% (Figure 2).


## 2.6. Land-Use and Interception

In semiarid areas of Mexico, matorral intercepts from 1.9% to 7.7% of gross rainfall [19], deciduous trees from 9% to 20%, and coniferous trees from 20% to 48% [12]. In urban areas, rainfall is intercepted by buildings and vegetation. Buildings intercept rainfall on rooftops and laterally on wall surfaces, modifying the hydrological response peak time [20]. Canopy shrub and grass saturation was estimated at 1.5 mm using the linear relationship between throughfall and steamflow, compared to forest saturation, which was observed at 5.0 mm [12].The edge effects present in isolated tree canopies might increase evaporation from wet canopies [21, 22]. Urban vegetation occurs as individuals or parklands, and therefore area-based interception percentages are smaller than reported for natural forests, grass, and possibly agroforests [12, 23–25]. Klaassen et al. [26] affirm that proximity to the forest edge affects both the interception storage capacity and the rate of evaporation of intercepted water, which cancel each other.


## 2.7. Soil

In Querétaro river basin the following different soils exist: litosol, vertisol, fluvisol, phaeozem, chernozem, castañozem, and yermosol (Figure3 and Table 3), according to World Reference Base for Soil Resources of the Food and Agriculture Organization (FAO), where clay, loam, and sand are the principal elements found in the basin. The maximum depth layer is found at 5 m with a rock basement detected, located at the north of the basin. Vertisol is located in the center of the basin, litosol at the south, north, and east, and phaeozem in the upper areas coinciding with the watershed mainly at the north and south; other soils are small area of fluvisol (north), chernozem (south) in Huimilpan municipality, and yermosol and castañozem.Figure 3
SHETRAN soils in the Querétaro river basin (mesh resolution: 2,100 meters). The existing soil composition in this basin isvertisol (52%), phaeozem (23%), litosol (18%), and fluvisol, castañozem, chernozem, and yermosol (7%). Vertisol principally at the center of the basin, atphaeozem north and south, and litosol found at the center, northwest, and southwest.


## 2.8. Precipitation

In this region the main rainfall volume occurs in summer (June to September), but the weather stations usually have breakdowns, resulting in incomplete databases that usually are filled by calculating missing values; to determine the mean precipitation 80 rainfall years have been analyzed. As a result dry years are considered between 257 and 550 mm/year.The maximum rainfall volume occurs between 5,850 and 6,800 hours, but May 2003 (3,700) and June (4,500) showed an increase in volume annual rainfall. The maximum event corresponds to 68 mm (24-hour rainfall register). The 100 cm registered precipitation is considered in this paper to evaluate the basin model calibration and runoff total volume. The spatial rainfall analysis was made with daily precipitation records, using for this 5 climate stations which were visited to evaluate if they comply with the World Meteorological Organization (WMO) principally related to obstructions at surroundings [27].The rainfall applied in the basin was the daily data obtained from the government climate stations in data station: 22027 (Carrillo), 22045 (Juriquilla), El Batan (22004), 22006 (Pueblito), and 22063 (Querétaro). The climate stations details and location are shown in Table1 and Figure 2.


## 3. Theory

A SHETRAN grid soil column model [15] is used to determine superficial and subsuperficial water interchanges [28]. Evaporation is calculated with the Penman-Monteith equation [29]:
(4)
E
p
=
R
n
Δ
+
(
ρ
c
p
δ
e
/
r
a
)
λ
(
Δ
+
γ
)
,
γ
=
P
c
p
σ
λ
,
LE
=
Δ
(
R
n
-
G
)
+
ρ
C
p
(
(
e
s
-
e
)
/
r
a
)
Δ
+
γ
*
,
γ
*
=
γ
(
1
+
r
c
r
h
)
,
where E
p is potential evapotranspiration, R
n is net radiation, Δ is rate of increase with temperature of the saturation vapor pressure of water at air temperature, ρ is density of air, δ
e is vapor pressure deficit of the air, r
a is aerodynamic resistance to transport of water vapor from the canopy to a plane 2 m above it, λ is latent heat of vaporization of water, γ is psychrometric constant, P is atmospheric pressure, and σ is ratio of density of water vapor to density of air (σ
≈
0.622).C
p is specific heat of air at constant pressure, LE is latent heat flux (Wm−2), R
n is net radiation (Wm−2), G is soil heat flux (Wm−2), ρ is air density (kgm−3), C
p is specific heat of dry air jkg−1 (°C−1), e
s is the saturation vapor pressure (kPa), e is actual vapor pressure of the air (kPa), and r
h is aerodynamic resistance to turbulent heat (sm−1). Canopy resistance r
c is bulk surface resistance that describes the resistance to flow of water vapor from inside the leaf, vegetation canopy, or soil to outside.The surface calculated with stomatal resistance,r
s, and leaf area index are
(5)
r
c
=
2
r
s
L
t
,
r
h
=
ln
⁡
⁡
(
(
z
-
d
)
/
z
M
)
ln
⁡
⁡
(
(
z
-
d
)
/
z
h
)
k
2
U
z
.
z is height above the ground surface for the wind speed measurement (m), d displacement (m), z
M roughness coefficient transference bending (m), z
h assumed roughness length governing the transfer of sensible heat from the surface, k von Karman constant diffusive turbulence, and U
z wind velocity at z in ms−1.


### 3.1. ETp Methods

The Thornthwaite [30] method is one of the most applied methods to determine a basin’s water balance and thus originally was described:
(6)
P
=
I
+
AET
+
OF
+
Δ
SM
+
Δ
GWS
,
where P is the total precipitation, I is the interception, AET is the evapotranspiration, OF is the overland flow, Δ
SM is the soil moisture storage change, and Δ
GWS is the groundwater ponding changes. P, I, and AET are the main water balance variables that consider changes in soil moisture and groundwater storage changes.In SHETRAN, the ET ratio, actual evaporation (AE), and potential evaporation (ETp) are the main parameters that increase or decrease the calculated annual runoff volume in the basin discharge and therefore it is necessary to calculate the ETp considering the database old weather stations, in which the temperature is one of the most important variables.Five different methodologies have been analyzed to determine ET values to calculate and compare the total discharge and runoff volume: Thornthwaite and Papadakys (monthly data, M), Blaney and Criddle, Hargreaves, and Hamon (daily data, D). These methods are briefly explained.Originally, the Thornthwaite method [30] was developed with rainfall-runoff data applied to different basins. Results were an empirical relationship between ETp, air temperature that involves the vapor flux, and heat balance:
(7)
ETp
=
16
(
10
T
I
)
a
,
where
(8)
I
=
∑
1
12
i
,
a
=
(
675
*
10
-
9
I
3
)
-
(
771
*
10
-
7
I
2
)
+
(
1792
*
10
-
5
I
)
+
0.49239
,
i
=
(
T
5
)
1.514
.
The Papadakys method is based on saturated vapor deficit (
e
0
-
e
), relative humidity, and temperature. The potential ET was proposed in cases where the climate database is not complete. In this method, Papadakys [31] proposed that:
(9)
ETp
=
5.625
(
e
0
(
t
max
⁡
)
-
e
0
(
t
min
⁡
-
2
)
)
,
where the vapor pressure at saturation, e
0, can be calculated with Bossen’s relationship using average temperature (t
med) in °C as follows:
(10)
e
0
=
33.8639
[
(
0.00738
*
t
med
+
0.8072
)
8
-
0.000019
g
g
g
g
g
g
*
(
1.8
*
t
med
+
48
)
+
0.001316
(
0.00738
*
t
med
+
0.8072
)
8
]
.
The Hargreaves method [32] is often used to evaluate potential ETp using temperature and solar radiation. The relationship is known as follows:
(11)
ETp
=
0.0023
(
t
med
+
17.78
)
*
R
0
*
(
t
max
⁡
-
t
min
⁡
)
0.5
,
where R
0 (also known as R
a) is the extraterrestrial solar radiation as a function of latitude. In Querétaro (México) this value is north 21°.The Blaney and Criddle method [33] is used in arid and semiarid regions to calculate ETp for periods of one month or greater. It is commonly used to estimate reference crop evapotranspiration in soil with water deficit using the temperature variable as follows:
(12)
ETp
=
p
*
[
(
0.46
*
T
med
)
+
8.13
]
,
where p
=
100 * (daylight hours per day/daylight hours per year).The Hamon method is applied mainly in places where there are no reliable climate databases. Thornthwaite developed an empirical method with mean temperature data and average daylight hours. Different authors found similarities between the empirical Thornthwaite relationship and saturated vapor pressure. The Hamon relationship is(13)
ETp
=
29.8
*
D
*
e
a
(
T
med
)
T
med
+
273.2
.
An average of 10 hours (D value) is applied; saturated vapor is thus calculated:
(14)
e
a
(
T
med
)
=
0.611
exp
⁡
⁡
(
17.3
T
med
T
med
+
237.3
)
.


### 3.2. Hydrological Model

SHETRAN is a physically based distributed hydrological model for water and sediment flow in basins. This model includes the Gash model [34] and Rutter model [35, 36] applying their rainfall interception model to vegetation and to different land-uses as well as to ETp phenomena. This software uses variables as vapor-transpiration components for variable saturated subsurface flow interception, overland flow through the numerical solution of the Saint-Venant equations, and interactions between surface and subsurface waters with the Richards equation; the mass and momentum are solved with differential equations in a three-dimensional finite difference grid [37]. The different applications of this tool have been demonstrated to represent water flow, land-use changes, soils production erosion, transport effect simulation [38], and contaminant transport by runoff in Birkinshaw and Ewen [39].The SHETRAN platform is based on a vertical column model representation dividing the basin into finite difference cells where upper cells are the superficial and subsuperficial water. At the end of column is the confined groundwater. The superficial layers are the vegetation and soil; over this layer are canopy and other interceptors. This tool is physically based because it allows solving the physically based equations of momentum and mass-energy conservation. Hence the parameters have physical meaning and allow evaluating the field measurements across the partial spatially distributed results.In the Querétaro river basin the different vegetation layers have been included with experimental results parameters [12]. The physically based model allows the change in forest cover to be modeled by changing the vegetation parameter values to represent the forest and urban centers condition. In this model a 2,100 × 2,100 m vegetation types mesh is used, composed as in Table 2.Table 2
Land-use SHETRAN parameter values in the Querétaro river basin. AE/PE is the relationship for actual and potential evaporation. Values of canopy storage capacity, leaf area index, and maximum routing depth are based on Guevara-Escobar et al., 2000 [10], Mastachi-Loza et al., 2010 [11], and Guevara-Escobar et al., 2007 [12].
Type
Vegetation type
Canopy storage capacity (mm)
Leaf area index
Max routing depth(m)
Total area(%)
AE/PE at field capacitya
1
Deciduous forest
5.0
6
1.5
21.00
1.10
2
Grass
1.5
6
1
20.00
1.10
3
Arable
1.5
4
0.8
45.00
1.05
4
Urban
0
0.3
0.5
13.50
1.00
5
Urban streets
0
0.3
0.5
0.50
1.00
aFinal calibrated value.Table 3
SHETRAN soil parameters applied to model.
Soilcategory
Soillayer
Soil type
Depth at base of layer (m)
Saturated water content
Residual water content
Saturated conductivity (m/day)
Van Genuchten-Alpha(cm−1)a
Van Genuchten-n
b
1
1
Silt loam(10% sand : 10% clay)
1
0.452
0.093
0.163
5.15
E
-
03
1.681
2
1
Loamy sand(85% sand : 6% clay)
1
0.370
0.075
0.467
1.99
E
-
02
1.793
3
1
Sandy silt loam(35% sand : 10% clay)
1
0.434
0.086
0.317
8.38
E
-
03
1.587
3
2
Silty clay loam(10% sand : 27% clay)
2
0.507
0.144
0.036
7.24
E
-
03
1.608
3
3
Bedrock
5
0.200
0.030
0.100
3.00
E
-
02
1.900
4
1
Clay loam(35% sand : 27% clay)
1
0.489
0.153
0.055
9.23
E
-
03
1.657
6
1
Sandy clay(52% sand : 40% clay)
1
0.499
0.233
0.029
1.07
E
-
02
1.879
7
1
Silty clay(10% sand : 40% clay)
1
0.529
0.212
0.019
6.54
E
-
03
1.531
8
1
Sandy loam(65% sand : 10% clay)
1
0.412
0.098
0.622
1.44
E
-
02
1.736
aVan Genuchten exponent α for soil moisture content-tension curve in cm−1.
bVan Genuchten exponent n for soil moisture content-tension curve.Runoff is solved by Saint-Venant equations which describe the one-dimensional flux of the wide wave and steady regime as followsContinuity equation (mass water conservation) is as follows:(15)
∂
Q
∂
x
+
∂
A
∂
t
=
q
.
Momentum conservation is as follows:
(16)
∂
Q
∂
t
+
∂
(
α
(
Q
2
/
A
)
)
∂
x
+
g
A
∂
h
∂
x
+
g
Q
|
Q
|
C
2
A
R
=
0
,
where t is time, x is measured distance along the channel (m), Q is discharge (m3s−1), A is hydraulic area (m2), q is tributary outflow (m3s−1), h is channel depth (m), C is Chezy coefficient (m0.5s−1), R is ratio hydraulic (m), and α is correction factor.


### 3.3. Model Calibration

The results of calibration are shown in Figure6, which shows the observed precipitation in millimeters and the total runoff observed in one year (8,760 hours). The Nash-Sutcliffe model efficiency coefficient is used to assess the predictive power of hydrological models; this method is defined as follows:
(17)
E
=
1
-
∑
t
=
1
T
(
Q
o
t
-
Q
m
t
)
2
∑
t
=
1
T
(
Q
o
t
-
Q
-
o
)
2
,
where Q
o is the observed discharge, Q
m is modeled discharge, and t is time. The Nash-Sutcliffe model was used to compare the results.Once the reference model represents the existing conditions in the basin, the different ETp methodologies have been introduced to the model with the results shown in Figure7. The result shows variability of the response of the basin applying the five different ETp methods.


## 3.1. ETp Methods

The Thornthwaite [30] method is one of the most applied methods to determine a basin’s water balance and thus originally was described:
(6)
P
=
I
+
AET
+
OF
+
Δ
SM
+
Δ
GWS
,
where P is the total precipitation, I is the interception, AET is the evapotranspiration, OF is the overland flow, Δ
SM is the soil moisture storage change, and Δ
GWS is the groundwater ponding changes. P, I, and AET are the main water balance variables that consider changes in soil moisture and groundwater storage changes.In SHETRAN, the ET ratio, actual evaporation (AE), and potential evaporation (ETp) are the main parameters that increase or decrease the calculated annual runoff volume in the basin discharge and therefore it is necessary to calculate the ETp considering the database old weather stations, in which the temperature is one of the most important variables.Five different methodologies have been analyzed to determine ET values to calculate and compare the total discharge and runoff volume: Thornthwaite and Papadakys (monthly data, M), Blaney and Criddle, Hargreaves, and Hamon (daily data, D). These methods are briefly explained.Originally, the Thornthwaite method [30] was developed with rainfall-runoff data applied to different basins. Results were an empirical relationship between ETp, air temperature that involves the vapor flux, and heat balance:
(7)
ETp
=
16
(
10
T
I
)
a
,
where
(8)
I
=
∑
1
12
i
,
a
=
(
675
*
10
-
9
I
3
)
-
(
771
*
10
-
7
I
2
)
+
(
1792
*
10
-
5
I
)
+
0.49239
,
i
=
(
T
5
)
1.514
.
The Papadakys method is based on saturated vapor deficit (
e
0
-
e
), relative humidity, and temperature. The potential ET was proposed in cases where the climate database is not complete. In this method, Papadakys [31] proposed that:
(9)
ETp
=
5.625
(
e
0
(
t
max
⁡
)
-
e
0
(
t
min
⁡
-
2
)
)
,
where the vapor pressure at saturation, e
0, can be calculated with Bossen’s relationship using average temperature (t
med) in °C as follows:
(10)
e
0
=
33.8639
[
(
0.00738
*
t
med
+
0.8072
)
8
-
0.000019
g
g
g
g
g
g
*
(
1.8
*
t
med
+
48
)
+
0.001316
(
0.00738
*
t
med
+
0.8072
)
8
]
.
The Hargreaves method [32] is often used to evaluate potential ETp using temperature and solar radiation. The relationship is known as follows:
(11)
ETp
=
0.0023
(
t
med
+
17.78
)
*
R
0
*
(
t
max
⁡
-
t
min
⁡
)
0.5
,
where R
0 (also known as R
a) is the extraterrestrial solar radiation as a function of latitude. In Querétaro (México) this value is north 21°.The Blaney and Criddle method [33] is used in arid and semiarid regions to calculate ETp for periods of one month or greater. It is commonly used to estimate reference crop evapotranspiration in soil with water deficit using the temperature variable as follows:
(12)
ETp
=
p
*
[
(
0.46
*
T
med
)
+
8.13
]
,
where p
=
100 * (daylight hours per day/daylight hours per year).The Hamon method is applied mainly in places where there are no reliable climate databases. Thornthwaite developed an empirical method with mean temperature data and average daylight hours. Different authors found similarities between the empirical Thornthwaite relationship and saturated vapor pressure. The Hamon relationship is(13)
ETp
=
29.8
*
D
*
e
a
(
T
med
)
T
med
+
273.2
.
An average of 10 hours (D value) is applied; saturated vapor is thus calculated:
(14)
e
a
(
T
med
)
=
0.611
exp
⁡
⁡
(
17.3
T
med
T
med
+
237.3
)
.


## 3.2. Hydrological Model

SHETRAN is a physically based distributed hydrological model for water and sediment flow in basins. This model includes the Gash model [34] and Rutter model [35, 36] applying their rainfall interception model to vegetation and to different land-uses as well as to ETp phenomena. This software uses variables as vapor-transpiration components for variable saturated subsurface flow interception, overland flow through the numerical solution of the Saint-Venant equations, and interactions between surface and subsurface waters with the Richards equation; the mass and momentum are solved with differential equations in a three-dimensional finite difference grid [37]. The different applications of this tool have been demonstrated to represent water flow, land-use changes, soils production erosion, transport effect simulation [38], and contaminant transport by runoff in Birkinshaw and Ewen [39].The SHETRAN platform is based on a vertical column model representation dividing the basin into finite difference cells where upper cells are the superficial and subsuperficial water. At the end of column is the confined groundwater. The superficial layers are the vegetation and soil; over this layer are canopy and other interceptors. This tool is physically based because it allows solving the physically based equations of momentum and mass-energy conservation. Hence the parameters have physical meaning and allow evaluating the field measurements across the partial spatially distributed results.In the Querétaro river basin the different vegetation layers have been included with experimental results parameters [12]. The physically based model allows the change in forest cover to be modeled by changing the vegetation parameter values to represent the forest and urban centers condition. In this model a 2,100 × 2,100 m vegetation types mesh is used, composed as in Table 2.Table 2
Land-use SHETRAN parameter values in the Querétaro river basin. AE/PE is the relationship for actual and potential evaporation. Values of canopy storage capacity, leaf area index, and maximum routing depth are based on Guevara-Escobar et al., 2000 [10], Mastachi-Loza et al., 2010 [11], and Guevara-Escobar et al., 2007 [12].
Type
Vegetation type
Canopy storage capacity (mm)
Leaf area index
Max routing depth(m)
Total area(%)
AE/PE at field capacitya
1
Deciduous forest
5.0
6
1.5
21.00
1.10
2
Grass
1.5
6
1
20.00
1.10
3
Arable
1.5
4
0.8
45.00
1.05
4
Urban
0
0.3
0.5
13.50
1.00
5
Urban streets
0
0.3
0.5
0.50
1.00
aFinal calibrated value.Table 3
SHETRAN soil parameters applied to model.
Soilcategory
Soillayer
Soil type
Depth at base of layer (m)
Saturated water content
Residual water content
Saturated conductivity (m/day)
Van Genuchten-Alpha(cm−1)a
Van Genuchten-n
b
1
1
Silt loam(10% sand : 10% clay)
1
0.452
0.093
0.163
5.15
E
-
03
1.681
2
1
Loamy sand(85% sand : 6% clay)
1
0.370
0.075
0.467
1.99
E
-
02
1.793
3
1
Sandy silt loam(35% sand : 10% clay)
1
0.434
0.086
0.317
8.38
E
-
03
1.587
3
2
Silty clay loam(10% sand : 27% clay)
2
0.507
0.144
0.036
7.24
E
-
03
1.608
3
3
Bedrock
5
0.200
0.030
0.100
3.00
E
-
02
1.900
4
1
Clay loam(35% sand : 27% clay)
1
0.489
0.153
0.055
9.23
E
-
03
1.657
6
1
Sandy clay(52% sand : 40% clay)
1
0.499
0.233
0.029
1.07
E
-
02
1.879
7
1
Silty clay(10% sand : 40% clay)
1
0.529
0.212
0.019
6.54
E
-
03
1.531
8
1
Sandy loam(65% sand : 10% clay)
1
0.412
0.098
0.622
1.44
E
-
02
1.736
aVan Genuchten exponent α for soil moisture content-tension curve in cm−1.
bVan Genuchten exponent n for soil moisture content-tension curve.Runoff is solved by Saint-Venant equations which describe the one-dimensional flux of the wide wave and steady regime as followsContinuity equation (mass water conservation) is as follows:(15)
∂
Q
∂
x
+
∂
A
∂
t
=
q
.
Momentum conservation is as follows:
(16)
∂
Q
∂
t
+
∂
(
α
(
Q
2
/
A
)
)
∂
x
+
g
A
∂
h
∂
x
+
g
Q
|
Q
|
C
2
A
R
=
0
,
where t is time, x is measured distance along the channel (m), Q is discharge (m3s−1), A is hydraulic area (m2), q is tributary outflow (m3s−1), h is channel depth (m), C is Chezy coefficient (m0.5s−1), R is ratio hydraulic (m), and α is correction factor.


## 3.3. Model Calibration

The results of calibration are shown in Figure6, which shows the observed precipitation in millimeters and the total runoff observed in one year (8,760 hours). The Nash-Sutcliffe model efficiency coefficient is used to assess the predictive power of hydrological models; this method is defined as follows:
(17)
E
=
1
-
∑
t
=
1
T
(
Q
o
t
-
Q
m
t
)
2
∑
t
=
1
T
(
Q
o
t
-
Q
-
o
)
2
,
where Q
o is the observed discharge, Q
m is modeled discharge, and t is time. The Nash-Sutcliffe model was used to compare the results.Once the reference model represents the existing conditions in the basin, the different ETp methodologies have been introduced to the model with the results shown in Figure7. The result shows variability of the response of the basin applying the five different ETp methods.


## 4. Results Analysis

Model calibration of 2003 (Figure6: model calibration) has been carried out by using the Nash-Sutcliffe efficiency factor (E = 90.25%) that represents the approximation of the method of calibration basin. The hydrometric data compared with those observed have differences attributable to the size of the basin, detention, or extractions.The obtained values of calibrated model were used to implement ETp results obtained from each method used. In Table3 (soil parameters) parameters used for the Querétaro river basin for different coverage and calibration for AE/PE ratio value are listed.From the calculated daily (D) and monthly (M) methods it can be appreciated that Hamon (D) and Thornthwaite (M) methods have water volume of 5.15 and 3.82 times the SHETRAN calculated volume whereas in the Blaney (D), Hargreaves (D), and Papadakys (M) methods the calculated volume is 31.10%, 37.92%, and 46.83% (68.97%, 59.12%, and 56.98% less) of the total, respectively (Table4). According to these results obtained the Papadakys method has the best similarity with the reference model but does not have significant differences with Blaney and Hargreaves results.Table 4
SHETRAN peak discharge and volume calculated data. Significant differences are observed from daily (D) and monthly (M) methods where Thornthwaite and Hamon overestimate volume; on the other hand Blaney, Hargreaves, and Papadakys underestimate the calculated total volume (31.10%, 37.92%, and 46.83%, resp.).
ETp method
SHETRAN
Blaney(D)
Hamon(D)
Hargreaves(D)
Papadakys(M)
Thornthwaite(M)
Volume × 1000 m3
76,273.78
23,667.24
393,299.12
31,176.94
32,809.94
291,799.56
Peak discharge (m3s−1)
60.58
18.84
191.24
22.97
28.37
132.35
Peak discharge percent
—
−68.90%
315.68%
−62.08%
−53.17%
218.47%
Total volume percent
—
−68.97%
515.64%
−59.12%
−56.98%
382.57%According to the results presented, methods that have an overestimation of runoff volumes (compared with SHETRAN model results) are Hamon and Thornthwaite (315.66% and 218.46%); Blaney, Hargreaves, and Papadakys (31.09%, 37.92%, and 46.82%) represent a smaller runoff volume. None of the shown methods has values close to the observed values in the hydrometric station.The Blaney method presents greater discharge values compared to the other methods in the representation of the ETp, where it is observed that small rainfall (10 mm accumulated in 24 hours) is intercepted by vegetation; the total volume in the year (Figures4(b) and 5(b)) has values of 147.14 mm (Blaney) and 191.53 (Hamon). Rainfall higher than 25 mm accumulated in 24 hours represents the runoff observed in the basin outlet.Cumulated and calculated values for canopy storage, transpiration, evaporation, and soil evaporation.
(a)
(b)
(c)
(d)The different total calculated values are shown to the main variables of the five ET methods shown.
(a)
(b)
(c)
(d)Figure 6
Model calibration: hydrometrical station discharge data, measured at the catchment outlet, and calculated SHETRAN discharge. In the ordinate axis the 2003 year hours are shown (0 hours is January 1st; 0:00 and 8,760 hours are December 31st, 24:00).Figure 7
SHETRAN model using different ET methods. Observed discharge data in Ameche hydrometrical station are shown.The 8,760 hours of modeling are shown in Figure4 in which the evaporation ranges are among 195.26 and 307.67 mm accumulated (Figures 4(a) and 5(a)); soil evaporation (Figures 4(c) and 5(c)) has values ranging from 110.67 (Hamon) to 120.09 (Blaney) mm per year and the total cumulative transpiration (Figures 4(d) and 5(d)) is 862.12 mm (Thornthwaite) and 973.07 mm (Papadakys). The results represent an approach to values that official sources have not been available in the basin and can be obtained through measurements and instrumentation. This means that disaggregating the results in each of the different parameters of the basin allows a better understanding of the processes that are involved in the water balance including dry-wet transitions.


## 5. Discussion

In this paper a digital terrain model was implemented in which the concentration of runoff is observed in the upstream to the lower parts of the basin, so the ETp phenomenon considers the location and concentration of flows in the lowlands, streams, and surface runways with elevations ranging from 1,787 masl to 2,710 masl.Different variables have been used that are involved in the hydrological cycle. The most important are land-use, canopy storage capacity, leaf area index, routing depth, actual-potential evaporation relationship (Table2), depth at base of layer, saturated water content, residual water content, saturated conductivity, and Van Genuchten parameters (Table 3). The variables have been recovered from past investigations in the basin but must be evaluated in order to have more and better results. The obtained parameters allow having a knowledge of the different variables used and they can help authorities to determine penalties or rewards to users that affect or contribute to its preservation. With this model the runoff volume is permanently determined with the climatological data or if they do not exist, they can be inferred from the model calculation.The rainfall-runoff results from 2003 shows that for the basin conditions the SHETRAN model results represent the hydrological phenomena of the Querétaro river basin. Calibration parameters of the relationship AE/PE are shown in Table2.These values indicate that the losses in the basin are significant due to the detention of the basin, as well as the water volume losses of water extractions for agricultural and cattle uses.Once the model has been calibrated, five different ETp methods analysis is made. The runoff maximum volume was obtained by the Thornthwaite and Hamon methods, while the methods results with smaller runoff volume were Blaney, Papadakys, and Hargreaves compared with the model calibrated.The results for the different methods show that the source of information based on variables such as saturated vapor deficit, relative humidity, temperature, solar radiation, and saturated vapor affects calculated volume, so care must be taken in evaluating and selecting the ETp method because the results may not be reliable. It is recommended to have extensive knowledge of the physical conditions of the basin, so that parameters such as AE/PE and the Strickler coefficient represent the phenomenon.This tool allows an initial scheme for the regulation applied to the variables involved in the hydrologic cycle and carries out government programs of the environment preservation with economic activities, with the results of physically based numerical models.


## 6. Conclusions

Knowledge of the hydrological parameters of a watershed is intended to enable authorities to carry out actions related to the rational use of resources and the preservation of hydrologic historical conditions as rainfall and runoff, and thus it is possible to regulate and restrict land-uses according to the effects on natural resources. In this paper data from the weather stations of the basin for a full year of observations have been used, which has made it possible to determine the conditions of evaporation in dry weather periods and rain season which is a short part of the year.The main limitations in this model are the fact that it has been applied to 2003 database that is the most complete year and is based on a basin of 2,700 square kilometers so that the methodology can be applied with different basins size. An advantage of using this model to calculate flow in ungagged basins is that once a hydrometric system is implemented, the climate data can be introduced according to the time-step observed as minutes or hours of accumulated rainfall.This research in North America represents semiarid conditions; the methodological differences in the evaluated methods determine the total evaporation variation values for ETp methods of 1 (daily) and 30 days (monthly). Climatological information represents the rainfall in this basin; the urban area is located at the center of the basin, so that the results suggest that urban area affects the calculated hydrograph shape and particularly the observed peak discharge.The main objective of this methodology is to infer results based on a physically based model on watersheds where there is missing data or no data and that has parameters obtained from research, so that it can be applied to determine the amount of water available from rainfall-runoff processes in places where the total rainfall volume is low and the continuous changes in the components of vegetation, soil, and superficial drainage system demonstrate the lack of analysis tools to predict the increase in runoff volumes, given the continuous land-use changes.According to the results, there is a big difference in the basin’s runoff volume and peak discharge, calculated by different ETp methods, and therefore the application of each method for purposes of water balance should be considered to avoid wrong calculated values.The obtained data show the differences between the methods used; however it is too important to pay special attention to the existing data in the basin, as there are variables that are very important in obtaining better results, such as the basin’s size, quantity and quality of the weather stations data, the water bodies within the watershed, and operating policies. The proposed method allows collecting the different parameters of a watershed; however for the best results better measurements are needed in continuous periods which can be validated to identify significant basin alterations.



---

*Source: 586872-2014-02-16.xml*