The Effect of Agricultural Development on Streamow

In recent decades, the amount of water in Urmia Lake has faced a considerable decrease. Several studies evaluated agricultural development as one of the main governing factors. The current study investigates the effect of agricultural development on the amount of streamow in one of the main Urmia Lake subbasins. To focus on agricultural development, a novel method using DEM and NDVI information was utilized to distinguish the agriculture class for selected years. Outputs of SWAT model for different years and scenarios revealed that agricultural development from the year 1977 to 2015, increased the annual evapotranspiration from 295 to 308 mm, and decreased the annual streamow from 317 mcm to 300 mcm for the irrigation scenario. In addition, the groundwater level uctuations captured by the SWAT, showed 1.9 meters drawdown in the Sarab plain from the year 1992 to 2003, which was in accordance with the observed trend during this period.


Introduction
In recent decades, the amount of water in the Urmia Lake faced a considerable decrease. Several studies evaluated different parameters responsible for this issue, including climate change, constructing dams on increase over the surface runoff. As another example of application of remote sensing in water-related studies, Thanapakpawin et al. (2007) applied DHSVM hydrological model to evaluate the effect of landuse change on the hydrological regime for irrigation and no-irrigation scenarios in the Mae Chaem basin located in the North-West of Thailand. It was concluded that in the irrigation scenario, changes from forest land into agricultural land decreased basin discharge, whereas in the no-irrigation scenario, changes from forest land into agricultural land increased basin discharge.
For Satellite image classi cation, several digital change detection techniques, namely univariate image differencing, image regression, image rationing, vegetation index differencing, principal components analysis, post-classi cation comparison, direct multi-data classi cation, change vector analysis, and background subtraction, have been developed (Singh, 1989). On the other hand, the standard deviation threshold, supervised and unsupervised analysis techniques, as shown by Nelson (1983), are among the most well-known techniques that have been frequently used to classify images. Among the aforementioned digital change detection techniques, univariate image differencing techniques was determined as the most accurate method (Singh, 1989).
The main satellite imagery projects frequently used in remote sensing studies are Spot, NOAA, MODIS, IKONOS, QuickBird, and Landsat. The basic advantage of Landsat imagery for change detection studies, among the others, is its long temporal coverage, which makes it the most suitable imagery product for the change detection studies (Xie et al., 2008). Besides, multispectral remote sensing is capable of capturing re ected radiation in the visible and near-infrared regions in which the most vegetation radiative changes occur (Muttitanon and Tripathi, 2005). The Landsat imagery covers a broad historical timeline since 1972. In this time period, Landsat imagery has been continuously providing scientists with accurate multispectral images that cover all parts of Earth. Furthermore, the capabilities of Landsat imagery have been enhancing in time by adopting more sophisticated sensors. The Landsat 1 to 3, covered the duration of 1972-1983 years, were equipped with Multispectral Scanner (MSS). The more advanced sensor was loaded on satellite 4 and 5 named Thematic Mapper (TM), which provides images with a spatial resolution of 30 meters. Landsat 6 was equipped with the next-generation sensor named Enhanced Thematic Mapper (ETM). However, it did not achieve orbit and went out of service. The Landsat 7 was launched in 1999, loading Enhanced Thematic Mapper Plus (ETM+) sensor armed with a panchromatic band with 15m spatial resolution. However, the project was unsuccessful due to the Scan Line Corrector (SLC) problem resulted in non-accurate imagery. The most recent Landsat product is Landsat 8, which has been in operation since 2013 and provides the most accurate Landsat images by taking advantage of its two instruments, namely the Operational Land Imager (OLI) and the Thermal Infrared Sensors (TIRS), which are able to provide seasonal images with the spatial resolution of 30 meters that can be improved to 15 meters using its panchromatic sensor (https://landat.gsfc.nasa.gov/). Numerous research have recently employed different Landsat imagery products to detect land-use and land cover changes in different parts of the world (Coulter et al., 2016;Devolent et al., 2017;Islam et al., 2018). Nevertheless, it is worth bearing in mind that the digital change detection techniques classify images solely based on the re ected radiation, and as a result, it cannot distinguish between different classes with the same spectral characteristics. This disadvantage becomes noticeable, more speci cally, in studies with the focus on the irrigation practices, in which agricultural elds need to be differentiated from rangeland elds. In this paper, a novel post-classi cation method using the NDVI (Normalized Difference Vegetation Index) and DEM (Digital Elevation Model) has been applied to images derived from Landsat 5 to reclassify vegetation class into agriculture and rangeland classes.
In order to quantify the land-use change effects on hydrology, different hydrological models have been utilized by different researchers, including PRMS (Leavesley et al., 1984), DHSVM (Wigmosta et al., 1994), and SWAT (Arnald et al., 1998). SWAT is a physically-based semi-distributed hydrological model developed to carry out basin-scale simulations. SWAT model can be used in watersheds without su cient or well-organized data (Arnald et al., 1998). The model is executable on a daily, monthly, or annual basis. For simulation, the whole watershed is divided into several Hydrologic Response Units (HRUs) based on the soil classes, land-use classes, and slope characteristics.
The model e cacy in studying the effects of land-use changes on stream ow has been widely investigated and proven (e.g., Du et al., 2013;Fohrer et al., 2001;Ghaffari et al., 2010). However, the effect of land-use changes, which is likely to be followed by changes in management practices, more speci cally in semi-arid regions, on both groundwater and surface ow has been rarely analyzed. Since the growth in agricultural elds and subsequent irrigation practices have been highlighted as the main factors that resulted in a decrease in both Urmia Lake area and level in several studies (Dariane and Eamen, 2017; Eamen and Alireza B. Dariane, 2014;Eimanifar and Mohebbi, 2007;Hassanzadeh et al., 2012), the effect of land-use change on irrigation practices from both surface and ground reservoirs must be taken into account. Different pieces of research have made efforts to solve groundwater equations in the SWAT model. The most famous method in this regard has been the use of MODFLOW to solve the groundwater ow equations and the SWAT for handling the surface ow equations (Dowlatabadi and Zomorodian, 2016;Kim et al., 2008). This approach was later developed as a coupled SWAT-MODFLOW model (Bailey et al., 2016). Although the coupled SWAT-MODFLOW model has some advantages, including geographically located HRUs and effective HRU mapping procedure between two models, but it also has a substantial drawback of not calibrating both models simultaneously. Accordingly, relying solely on the integrated SWAT hydrological model for both groundwater and surface water uctuations has been suggested (Ahmadzadeh et al., 2016). Furthermore, SWAT model has shown its capability of assessing the management procedures on the hydrologic responses in terms of groundwater and surface water alterations (Baker and Miller, 2013).
The main objective of this paper is to study the land-use change effects on the amount of stream ow in Ajichai basin. In order to ful ll this objective, different land-use maps were generated using satellite data. Afterward, SWAT model was calibrated and validated for one of the land-use maps and then executed for the other generated land-use maps. In order to accurately simulate hydrological processes within the basin, we applied the irrigation practices to the basin with the help of the study area's aquifer map and water usage reports provided by the energy ministry of Iran. In addition, we aimed to compare groundwater drawdown simulated by SWAT with the observed values to show SWAT capability in simulating groundwater drawdown if irrigation data is accurately imported to the model.

Study Area
The case study is Ajichai basin above Merkid hydrometric station (Fig. 1). Ajichai located in East Azerbaijan, Iran, is one of the biggest sub-basins in Urmia Lake and has a pivotal role in feeding the Lake, mainly through early spring oods and snowmelts. The basin is situated between 46º 33´ and 47º 54é astern longitudes and 37º 36´ and 38º 24´ northern latitudes. The basin area and length of main river branch above Merkid are 5460 km 2 and 111 km, respectively.

Data
The climatic variables needed for executing the SWAT model are daily precipitation, daily minimum and maximum temperature, solar radiation, wind speed, and relative humidity. Nevertheless, the SWAT model is executable, relying solely on daily minimum and maximum temperature as well as the daily precipitation. In this study, six precipitation stations and two temperature stations were chosen as data stations. These stations were selected according to their data availability and proper distribution within the case study area. Figure 1 and Table 1 respectively show the dispersion of the stations and information regarding each station.
The elevation of the basin ranges between 1503 to 3724 meters above sea level ( Fig. 1) and is captured using DEM data obtained from the Advanced Spaceborne Thermal Emission and Re ection Radiometer  . It should be noted that three satellite images (three different paths and rows according to the USGS website) were obtained for each year due to extensive coverage of the case study. Afterward, the satellite images were processed to generate land-use maps for the selected years. The land-use map generation procedure using the obtained satellite images is discussed in the following sections.

Images classi cation
In order to classify images, Environment software for Visualizing Images (ENVI) was used. It has been concluded that the pre-classi cation process is necessary for a classifying scheme (Coppin et al., 2004). In this study, image enhancement techniques that are available in ENVI software were utilized to enhance the visual quality of images. Pre-classi cation steps, such as radiometric and atmospheric correction, were applied in advance to prepare images for classi cation. Overall, classi cation methods are categorized into two main groups, namely unsupervised and supervised. In unsupervised method, pixels are divided into separated clusters based on their re ectance characteristics. K-means and ISO (Iterative Self Organizing) are examples of widely used algorithms in this regard. On the other hand, in the supervised classi cation method, classi cation is carried out based on samples that are attributed to each class by the user. Mostly used supervised classi cation algorithms are Maximum Likelihood, Minimum Distance, and Arti cial Neural Network (ANN). In this study, the ANN algorithm was employed to classify the images into four classes, namely residential, water, soil-pasture, and vegetation. To do so, samples were selected by drawing polygons representing each class. Representative polygons for each class are de ned as ROI (Region of Interest) in the ENVI software. It should be noted that for classi cation results to be bene cial for change detection studies, classi cation accuracy should be evaluated initially (Owojori and Xie, 2005). Accordingly, to validate the accuracy of the classi cation method, 200 well-distributed points in each image were taken into account as the ground truth data. The confusion error matrix in the ENVI software was employed to make a comparison between classi cation results and ground truth data. As said before, we had to use three images for each year due to the huge coverage of our study area. Thus, the abovementioned steps proceeded for every three images related to each year. Afterward, mosaic tools in ArcGIS software were used to merge images for each year.

Splitting vegetation class into agriculture and rangeland
As discussed before, the increase in agricultural elds and subsequent irrigation practices have been considered as one of the main factors resulting in a decrease in the Urmia Lake area and level. Therefore, we needed to separate the agriculture class from the vegetation class to detect where irrigation practices should be applied. To ful ll this purpose, we employed a novel technique for splitting the vegetation class into agriculture and rangeland classes. In this technique, the DEM image and NDVI index are used simultaneously to distinguish agriculture class from rangeland class.
NDVI is a well-known index for detecting vegetation (Lunetta et al., 2006). Besides, the NDVI index is also able to differentiate between the agricultural elds and rangeland elds in a region as it is related to the vegetation density. Generally, NDVI values between 0 to 0.1, 0.2 to 0.3, and 0.3 to 1 are in accordance with barren, rangeland, and denser vegetation elds, respectively (Gross, 2005). Additionally, DEM images can be used for splitting agricultural elds from rangeland elds since DEM images represent topographic information that can be related to vegetation characteristics (Moore et al., 1991). Meanwhile, we expect rangeland elds to be located at higher elevations than agricultural elds. Therefore, the DEM image was used to specify an elevation threshold to stand for various elevations that are expected for agricultural and rangeland elds. In order to use both NDVI and DEM information to ner classify the vegetation class, a decision tree process was developed using ENVI software. In this regard, it was necessary to rstly produce NDVI images. To do so, we initially rescaled pixels values of our multispectral satellite images to lie between 0 and 1. Afterward, NDVI values were computed for each image by applying the NDVI equation using the Band math tool in the ENVI software. In Fig. 3 the resulted NDVI image is shown.
In order to apply our decision tree, two-band combined NDVI-DEM images were required. As a result, the DEM image was merged with the generated NDVI image. Having two-band combined NDVI-DEM images, the threshold value of 0.3, as suggested by (Gross, 2005), was applied on the NDVI band to separate the agriculture elds from the other classes. Subsequently, to further exclude the agriculture elds from the rangeland elds, the elevation threshold was used on the DEM image. An intuitive scheme of the decision tree discussed above is demonstrated in Fig. 4.
After splitting the vegetation class into agriculture and rangeland classes, ArcMap software was used again to eventually produce land-use maps for each year containing ve classes, namely residential, water, soil-pasture, agriculture, and rangeland.

Hydrologic Modelling using SWAT model
As mentioned earlier, in order to evaluate how different land-use practices affect river discharge, we used SWAT hydrologic model. SWAT is a physically-based hydrological model that is able to simulate different hydrological processes, including stream ow. The model uses general hydrological equations to simulate stream ow and the other hydrologic outputs. SWAT is a continuous-time model and is capable of simulating hydrological processes on a daily, monthly, or annual time scale. SWAT divides each basin to ner sub-basins, and each sub-basin to several HRUs based on similar land-use, soil, and slope characteristics. For simulation, the amount of water stored in the soil, stream ow, nutrient cycle, sediment, plant growth, and management practices are initially calculated for each HRU, and subsequently computed for each sub-basin using a weighted average. SWAT requires different inputs to proceed, including DEM, location of hydrometric stations, rainfall and temperature data, soil and land-use map, and management data (Arnald et al., 1998).

Model setup
In this study, ArcSWAT interface for SWAT2012 was utilized for the simulation processes. To set up the model for our case study, SWAT requires different inputs, ranging from land-use map, soil map, DEM, climate inputs, and management information. Depending on the study's purpose, the management information can be imported into the model or can be neglected. SWAT has shown its robustness in taking into account the management practices (Neitsch et al., 2011). Furthermore, as noted before, since this study's main focus is to consider the probable effect of the management operations on the amount of stream ow, management information needs to be involved.
In practice, the basin was delineated using ArcSWAT interface. The interface carries out this process automatically by importing the DEM, drainage area, and the location of the hydrometric station at the outlet of the basin. Afterward, the generated land-use map and soil map derived from the FAO were imported to stand for the HRUs characteristics. In addition, the agriculture class is further divided into the main crops of the study area, which are spring wheat, spring barley, green beans, alfalfa, potato, tomato, and apple. Consequently, a total of 44 sub-basins and 556 HRUs were delineated. Precipitation data for six rainfall gauging stations were obtained from the Iran Water Resources Management Organization, and temperature (minimum and maximum) and two synoptic stations were obtained from the Iran Meteorological Organization.

Irrigation
Since one of the main objectives of the current study is to tackle land-use change and the resulted irrigation practice changes on the amount of stream ow, data and statistics have been assessed to provide insights into the cultivation patterns and the irrigation sources.
We used data from the National Water Document of Iran (NWDI) to gain information concerning the cultivation pattern for our case study. This document, for a various number of basins in Iran, represents planting and harvesting date for the main crops.
Another required irrigation data is the source of irrigation that can be speci ed in SWAT model as reach, shallow aquifer, deep aquifer, reservoir, or unlimited source outside basin (Neitsch et al., 2011). To determine the irrigation source for each sub-basin, we used the report on water master plan studies for Urmia Lake developed by the ministry of energy of Iran. Based on this report, reach and aquifer were detected as two primary sources of irrigation for the study area. In addition, on a sub-divided map, the utilization percentage of each of these primary irrigation sources has been determined. We used this map and its utilization percentage information to determine irrigation sources for our sub-basins. Moreover, in order to better assign irrigation sources, we matched the basin's aquifer map with the modeled basin to detect which sub-basins are more likely to be fed from the aquifer (Fig. 5). By analyzing data and maps mentioned above, we allocated the irrigation source codes to each sub-basin.

Sensitive parameters
In this study, we used SWAT-CUP software ( In this software, the calibration and validation processes need to start with selecting the sensitive parameters to the stream ow. Furthermore, the changing range of each parameter, and the type of change to be applied to the parameters should be determined. SWAT is equipped with three different types of sampling methods to replace parameter values with the new ones within each parameter's changing range. These are relative change, replace, and absolute methods. In the relative change method, the parameter value of each HRU is replaced by a fraction of its initial value after each iteration. In the replace method, the new values of each parameter replace the former values for all HRUs. Finally, in the absolute method, a given value is added to the former parameter value, and the new value replaces the parameter value for all HRUs. It should be noted that for spatial parameters such as curve number (CN2), the relative change method is suggested as these parameters should be set independently from one HRU to another (Abbaspour et al., 2007a). Moreover, it is also important to know that the changing range of some parameters is limited to a predetermined range and cannot exceed that range. For instance, CN2 value is closely correlated with soil and land-use maps that are inputs to the model. As a result, by changing the CN2 values by more than 20 percent, the inputs, in this case, soil and land-use maps information will also change. According to the above-mentioned explanations, the sensitive parameters to the stream ow, their sampling method and their changing range were chosen as shown in Table 2.

Different parameterization schemes
In order to optimize the calibration and validation processes, two schemes were considered in the parameterization process, which is described below.
Modifying or not modifying the changing range after each set of iterations As explained in the SWAT-CUP manual, the changing range of parameters is updated after each set of iterations. Nonetheless, in many cases, the suggested range for the next set of iterations exceeds the model's acceptable range for the parameter. In these cases, it was suggested to manually modify the suggested range to the acceptable range in order to optimize the calibration process (Abbaspour et al., 2007b). However, in this study, we also explored not modifying the suggested range for the next set of iterations. The reason is that in case of unacceptable values, SWAT replaces the unacceptable values by either the maximum or minimum value of that parameter. By not limiting the proposed range, we explored new combinations of parameters values. Since the parameterization is a collaborative process, this exploring procedure might produce more favorable outcomes.

Choosing the objective function
Changing range proposed by the SWAT-CUP for the next set of iterations also depends on the selected objective function for the optimization. We explored and Nash-Sutcliffe E ciency Index (NS) as two different objective functions in this study.
By applying different parameterization schemes mentioned above, the changing range modi cation, and NS as the objective function generated better results.

Generated land-use maps
Based on the explanations brought in Sect. 3.1, the land-use maps were generated for the years 1977, 1993, 2005, and 2005. Each land-use map was classi ed into 5 classes, namely water, residential, soilpasture, agriculture, and rangeland (Fig. 6).

Changing trend of agriculture and residential areas
The agriculture and residential areas of each image were computed using the aforementioned maps for the selected years (Fig. 7a). After calculating the agricultural and residential land areas for each year, the change percentage of each class was computed as shown in Fig. 7b.
According to Figs. 7, although agricultural land expansion is viewed in all the selected years, in some periods, the expansion trend is steeper. According to Fig. 7a, agricultural land area for years 1977, 1993,   2005, and 2015 is 218, 313, 799, and 896 km 2 , respectively. Moreover, according to Fig. 7b, the R 2 agricultural land area faced 44%, 267%, and 312% increase from the year 1977 to the years 1993, 2005, and 2015 respectively. As it can be seen from Fig. 7b, the agricultural land did not grow remarkably from the year 1977 to the year 1993 due to the Iran-Iraq war that began in 1980 and ended in 1988. The war resulted in a slight increase and in some cases a decrease in agricultural land area in different parts of Iran. As opposed to this period, the agricultural land area surged up from the year 1993 to 2005. Generally speaking, a part of this increase can be related to the population growth in this time period, and accordingly, increase in need for agriculture, however the government policy toward self-su ciency in important goods such as wheat is the main reason behind this sharp increase. The same policy has caused devastating impacts on water scarcity in the area and Urmia Lake drying up crises.
Furthermore, other factors, such as desirable weather conditions, lack of proper careers in different elds, and the availability of agricultural facilities in this time period also caused an increase in agricultural land. During this period, the governmental policies in Iran toward self-su ciency caused a 155% expansion in agricultural land area. As the drought period began between the years 2005 and 2015, the agricultural land growing trend, once more, faced a decrease as a consequence of precipitation decrease, the decline in groundwater level, and other obstacles that farmers faced at that time.
In the study area, as seen in Fig. 7a, the residential area has always increased as analogous to the other regions of Iran. According to Fig. 7a, the residential area was around 3 km 2 in the year 1977 and increased to 9, 21, and 44 km 2 for the years 1993, 2005, and 2015 respectively. Additionally, as shown by

Calibration and Validation
After generating the land-use maps for the proposed years (i.e. 1977, 1993, 2005, and 2015), the 1993 land-use map was taken as the primary map for proceeding the calibration and validation steps. Therefore, the model was initially run for the 1993 land-use map on monthly time scale and then calibrated against monthly discharge observations of the Merkid hydrometric station from the year 1992 to 1999, and later, validated for the year 2000 to 2003. Figure 8 shows the monthly average of the simulated and observed stream ow for both the calibration and validation periods. Furthermore, the R 2 and NS values for the calibration and validation periods were found to be 0.74, 0.73 and 0.74, 0.64 respectively.
Next, to evaluate the effect of land-use change on the amount of stream ow for other years, the 1993 land-use map was replaced by 1977, 2005, and 2015 land-use maps and the calibrated model was executed again. Advantage of this approach to study the effect of land-use change on the amount of stream ow is that the land-use map is the only variable in different model runs. In contrast, in some other land-use change effect studies, the land-use update option in the SWAT model is used to stand for the land-use changes in time. However, in those approaches, the variation of other climatic variables, and more speci cally rainfall, also take parts and are responsible for the stream ow changes. Therefore, the changes in the simulated stream ow cannot be solely related to the changes in land-use.

Comparing the simulated and observed groundwater level
As explained in Sect. 3.2.2, the aquifer map of the study area was utilized to more accurately assign the groundwater irrigation source to each sub-basin. In order to show the capability of the SWAT model in simulating the amount of water derived from groundwater source in each sub-basin, accumulated groundwater drawdown simulated by SWAT for the Sarab plain, one of the central plains of the study area, was compared with the observation values from 1992 to 2003 (Fig. 9).
According to Fig. 9, SWAT proves its ability to capture the groundwater drawdown trend. However, in some years, a noticeable gap is evident between the simulated and observed values. This gap is mainly as a result of not calibrating the simulated groundwater drawdown against the observed values, noting that the stream ow calibration was the priority of this research. Nevertheless, the calibration of the simulated groundwater drawdown against the observation values is suggested for future studies. In order to also address the effect of irrigation practices on these changes, the model was run, once more According to Fig. 11a, for both scenarios, ET increased in all simulated years as a result of an increase in the agricultural land. However, the ET values are higher for the irrigation scenario than the no-irrigation scenario due to the fact that irrigation practices make more water accessible for evapotranspiration. Additionally, as a re ection on the increase in ET, the stream ow values declined in time for both scenarios. However, this decrease is more evident for the irrigation scenario due to more water consumption, both from surface and groundwater sources (Fig. 11b).

Conclusions
The main objective of this study was to investigate the effect of land-use change (i.e., agricultural development) on the amount of stream ow considering surface-groundwater interactions. To do so, we initially produced land-use maps for the years 1977, 1993, 2005, and 2015, including water, residential, soil-pasture, and vegetation classes. Afterward, using a novel approach we further classi ed the vegetation class into agriculture and rangeland classes to apply the management practices only to the agriculture class. We concluded that the agricultural land areas have continuously increased during the years considered in a way that the agricultural land areas faced a 312% increase from the year 1977 to the year 2015.
In the second part of the current study, in order to see how this increase in agricultural land affected the ET and stream ow values, we, rstly, set up and calibrated the SWAT model for the 1993 land-use map, and eventually, ran the model for the other land-use maps individually using two scenarios of irrigation and no-irrigation for agricultural areas. Results indicated that the amount of ET increased and the amount stream ow decreased in both scenarios from the year 1977 to the year 2015 due to agricultural developments in the basin. Nonetheless, the increase in the ET values and the decrease in the stream ow values were more evident in the irrigation scenarios due to more water consumption, both from surface and groundwater sources. According to the results obtained in this research, the agricultural development could be partially responsible for the Lake Urmia desiccation.

Declarations
Declaration of competing interest Figure 1 Dispersion of precipitation, temperature, and hydrometric stations Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 2
Soil map of the case study Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. NDVI image of the case study. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Decision tree to divide the vegetation class into agriculture and rangeland classes Figure 5 The basin's aquifer map pictured on the modeled basin Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.       Comparison of (a) annual average ET and (b) annual average stream ow for irrigation and no-irrigation scenarios