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

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´etaro 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.


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; 2 Journal of Applied Mathematics 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.

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 ( Figure 1). Its coordinates are UTM-WGS84 14N ( ) = 358,380; ( ) = 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.
The total rainfall registered in climate stations in 2003 (Table 1) 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].

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.

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 where ( ) is the value of the elevation in the position , ( ) is the function that describes the dataset , ( ) is the stochastic local variation in the dataset after the ( ) variation, and , is the no tendencies variability dataset (noise). When the dataset has no tendencies ( ) = 0, then where ℎ is the distance among the evaluated points and the variogram (ℎ) as follows: Var (3) Figure 1 shows the resulting digital terrain model.

Runoff and Rivers.
In this basin there are three main rivers (Figures 2 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]. 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).

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 areabased interception percentages are smaller than reported for natural forests, grass, and possibly agroforests [12,[23][24][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.

Soil.
In Querétaro river basin the following different soils exist: litosol, vertisol, fluvisol, phaeozem, chernozem, castañozem, and yermosol ( Figure 3 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.

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 Table 1

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]: where is potential evapotranspiration, is net radiation, Δ is rate of increase with temperature of the saturation vapor pressure of water at air temperature, is density of air, is vapor pressure deficit of the air, 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, is atmospheric pressure, and is ratio of density of water vapor to density of air ( ≈ 0.622).
is specific heat of air at constant pressure, LE is latent heat flux (Wm −2 ), is net radiation (Wm −2 ), is soil heat flux (Wm −2 ), is air density (kgm −3 ), is specific heat of dry air jkg −1 ( ∘ C −1 ), is the saturation vapor pressure (kPa), is actual vapor pressure of the air (kPa), and ℎ is aerodynamic resistance to turbulent heat (sm −1 ). Canopy resistance 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, , and leaf area index are is height above the ground surface for the wind speed measurement (m), displacement (m), roughness coefficient transference bending (m), ℎ assumed roughness length governing the transfer of sensible heat from the surface, von Karman constant diffusive turbulence, and wind velocity at in ms −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: where is the total precipitation, 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. , , 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:  The Papadakys method is based on saturated vapor deficit ( 0 − ), 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: where the vapor pressure at saturation, 0 , can be calculated with Bossen's relationship using average temperature ( med ) in ∘ C as follows: The Hargreaves method [32] is often used to evaluate potential ETp using temperature and solar radiation. The relationship is known as follows: where 0 (also known as ) 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: where = 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 An average of 10 hours ( value) is applied; saturated vapor is thus calculated:

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 vaportranspiration 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. Runoff is solved by Saint-Venant equations which describe the one-dimensional flux of the wide wave and steady regime as follows Continuity equation (mass water conservation) is as follows: Momentum conservation is as follows: where is time, is measured distance along the channel (m), is discharge (m 3 s −1 ), is hydraulic area (m 2 ), is tributary outflow (m 3 s −1 ), ℎ is channel depth (m), is Chezy coefficient (m 0.5 s −1 ), is ratio hydraulic (m), and is correction factor.

Model Calibration.
The results of calibration are shown in Figure 6, 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: where is the observed discharge, is modeled discharge, and 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 Figure 7. The result shows variability of the response of the basin applying the five different ETp methods.

Results Analysis
Model calibration of 2003 (Figure 6: model calibration) has been carried out by using the Nash-Sutcliffe efficiency factor ( = 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 Table 3 (soil parameters) parameters used for the Querétaro river basin for different coverage and calibration for AE/PE ratio value are listed.  (Table 4). 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.
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 (Figures 4(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.   of the basin allows a better understanding of the processes that are involved in the water balance including dry-wet transitions.

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 landuse, canopy storage capacity, leaf area index, routing depth, actual-potential evaporation relationship (Table 2), 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 Table 2.
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  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.

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.

Conflict of Interests
In this research the results using SHETRAN are shown, applied in the Querétaro river basin, in order to propose a physically based methodology for water balance analysis in this catchment. With the obtained results five methodologies used in basins have been evaluated where data are scarce. The results suggest that all methodologies have differences between calculated and observed data. In Mexico the government agency of regulations on water issues is Comisión Nacional del Agua (CONAGUA), which can use the scientific contribution developed by research centers and universities. The results are intended to show the shortcomings of traditional methods in the calculation, but this does not affect the point of view of the authorities because all the methods are compared on the same analysis. Therefore the authors declare that there is no conflict of interests regarding the publication of this paper.