Transboundary Basins Need More Attention: Anthropogenic Impacts on Land Cover Changes in Aras River Basin, Monitoring and Prediction

Changes in land cover (LC) can alter the basin hydrology by affecting the evaporation, infiltration, and surface and subsurface flow processes, and ultimately affect river water quantity and quality. This study aimed to monitor and predict the LC composition of a major, transboundary basin contributing to the Caspian Sea, the Aras River Basin (ARB). To this end, four LC maps of ARB corresponding to the years 1984, 2000, 2010, and 2017 were generated using Landsat satellite imagery from Armenia and the Nakhchivan Autonomous Republic. The LC gains and losses, net changes, exchanges, and the spatial trend of changes over 33 years (1984–2017) were investigated. The most important drivers of these changes and the most accurate LC transformation scenarios were identified, and a land change modeler (LCM) was applied to predict the LC change for the years 2027 and 2037. Validation results showed that LCM, with a Kappa index higher than 81%, is appropriate for predicting LC changes in the study area. The LC changes observed in the past indicate significant anthropogenic impacts on the basin, mainly by constructing new reservoir dams and expanding agriculture and urban areas, which are the major water-consuming sectors. Results show that over the past 33 years, agricultural areas have grown by more than 57% from 1984 to 2017 in the study area. Results also indicate that the given similar anthropogenic activities will keep on continuing in the ARB, and agricultural areas will increase by 2% from 2017 to 2027, and by another 1% from 2027 to 2037. Results of this study can support transboundary decision-making processes to analyze potential adverse impacts following past policies with neighboring countries that share the same water resources.

. Schematic map of the study area.

Overview of Data Collection and Research Methods
The main research questions that this study is trying to address are: (1) Whether land class modeling can accurately simulate LC changes in ARB, (2) given that similar anthropogenic activities occur in the future, how would the LC change over the next two decades be. To answer these questions, the following five-step research is conducted: (1) Preparing satellite images, (2) producing LC maps, (3) determining transition area matrix and simulating LC maps, (4) predicting LC map in the future using the Markov Chain model, and (5) investigating future water consumption. The flowchart of the analysis process explained in Figure 2, and these steps are outlined in the following sections. Along with national projects in ARB, a few joint international projects established, mostly between Turkey and Armenia, Iran and Azerbaijan, and Iran and Armenia. For instance, a run-off-the-river hydropower project is under investigation and negotiation between Iran and Armenia [36]. Accordingly, the river water quality and quantity expected to be altered from natural flows and therefore cause changes in the LC.
Remote Sens. 2020, 12, 3329 4 of 22 This study aims to investigate LC change of the Aras River in the territory of Armenia and the Nakhchivan Autonomous Republic. In the political divisions of Armenia, the study area included the provinces of Ararat, Syunik, Vayost Dzor, Armavir, Kotayl, Aragatsotn, and Shirak, and in the political divisions of Azerbaijan, the whole Nakhchivan Autonomous Republic also studied (see Figure 1).

Overview of Data Collection and Research Methods
The main research questions that this study is trying to address are: (1) Whether land class modeling can accurately simulate LC changes in ARB, (2) given that similar anthropogenic activities occur in the future, how would the LC change over the next two decades be. To answer these questions, the following five-step research is conducted: (1) Preparing satellite images, (2) producing LC maps, (3) determining transition area matrix and simulating LC maps, (4) predicting LC map in the future using the Markov Chain model, and (5) investigating future water consumption. The flowchart of the analysis process explained in Figure 2, and these steps are outlined in the following sections.

Preparing Satellite Image and Preprocessing
Landsat 5 (TM), 7 (ETM), and 8 (OLI) satellite imageries (Collection 1, Level-1) in 1984, 2000, 2010, and 2017 were acquired from Earth Explorer website (https://earthexplorer.usgs.gov/). It should be noted that the coordinate system used in the entire image data in this study is UTM, and all of these images are projected to WGS_1984_UTM_Zone_38N. As shown in Table 1, four windows of Landsat satellite imagery applied at each studied period (with 168-33, 169-32, 169-33, and 170-32 patch-row codes) and mosaicked for covering the study area. To reduce the impacts of seasonal and phenological variation, the image pairs were selected in the same season or nearly the same season timing for all the studied years (Table 1). Therefore, temporal ranges for most scenes were restricted to within about one month. Land Cloud Cover (less than 5%) was also a critical criterion in selecting

Preparing Satellite Image and Preprocessing
Landsat 5 (TM), 7 (ETM), and 8 (OLI) satellite imageries (Collection 1, Level-1) in 1984, 2000, 2010, and 2017 were acquired from Earth Explorer website (https://earthexplorer.usgs.gov/). It should be noted that the coordinate system used in the entire image data in this study is UTM, and all of these images are projected to WGS_1984_UTM_Zone_38N. As shown in Table 1, four windows of Landsat satellite imagery applied at each studied period (with 168-33, 169-32, 169-33, and 170-32 patch-row codes) and mosaicked for covering the study area. To reduce the impacts of seasonal and phenological variation, the image pairs were selected in the same season or nearly the same season timing for all the studied years (Table 1). Therefore, temporal ranges for most scenes were restricted to within about one month. Land Cloud Cover (less than 5%) was also a critical criterion in selecting the images. To satisfy the preprocessing requirement for change detection, multitemporal image registration, radiometric, and atmospheric corrections were accomplished first. For all scenes, preprocessing included checking geometric accuracy (minor modifications were necessary for them), as well as data stacking to place all scenes in a single composite dataset. Radiometric normalization, as well as atmospheric correction, was conducted by converting digital number to at the top of atmospheric reflectance. Finally, all scenes were stacked and used as inputs to the classifier for generating LC maps. All preprocessors were performed in TerrSet software environment.

Producing LC Maps and Accuracy Assessment
A multi-step process involving (1) visual interpretation with true and false color composites of the images and Google Earth images [40,41], (2) images transformation maps i.e., PCA, TEXTURE, TASSCAP, NDVI, and normalized different water indices (NDWI) [42][43][44][45][46][47], and (3) supervised classification (i.e., maximum likelihood, as one of the most popular supervised classification methods used with remote sensing image data [45]) were used in IDRISI-TerrSet (version 18) software for the preparation of LC maps. It should be noted that supervised classification requires the introduction of training samples [48], so we used the visual interpretation and images transformation analysis to generate training samples. In this regard, the composite images, Google Earth images and the image transformation maps were overlapped together for identifying certain training sample in shape of small polygons. To this end, first all images were geometrically corrected based on high-resolution images of Google Earth using an image-to-image method (with an RMS error < 0.5 pixel). Then, for each of the studied LC, 180 training samples were extracted from Google Earth images. The training samples were considered the same for all the studied years, and in their selection, it was considered that the LC displayed by that point/polygon did not change over the studied period. For instance, we selected 180 training samples for the agricultural class, during which the 32-year period (1984-2016) persisted as agricultural and had not changed to other LCs. It should be noted that PCA, TEXTURE, TASSCAP, and NDVI transformation accompanied by the composites images were also used to produce these training samples for forest, agricultural, rangeland, and bare land classes [44,45].
NDWI values derived from satellite images are commonly and successfully utilized in surface water body detection and mapping (Equation (1); [49]). It is calculated as follows: where Green is the 3rd band and NIR is the 5th band of Landsat-8 OLI. NDWI values (−1 ≤ NDWI ≤ +1) can be considered as a basis for classifying water and non-water LCs [50,51]. Also, the shallow water, which represents rivers well, can be discriminated from deep water based on the NDWI values. It is because the shallow water has higher reflectance in the green and red bands due to the bottom effect and can be classified as a separate water class [51]. To this end, NDWI threshold values can be designated both manually and automatically [51][52][53]. Therefore, we experimentally obtained the thresholds for separating non-water areas from shallow and deepwater classes as follows: Non-water: −1 ≤ NDWI < 0, shallow water: 0 ≤ NDWI < 0.3, and deep water 0.3 ≤ NDWI ≤ 1.
Manual digitization was also carried out to obtain higher accuracy in the classification of residential areas [54]. Finally, seven LC types i.e., (1) forest including primary forest, planted forest, woodland, and other naturally regenerated forest, (2) agricultural lands, (3) rangelands (i.e., areas that are covered by shrubland and grassland), (4) bare lands, corresponding to bare soil (i.e., areas that have no vegetation cover except water), barren areas of rock, sand, or snow with NDVI close to 0 [44], (5) deep water, (6) shallow water (river), and (7) residential areas were classified in the study area.
Ground points were used to assess the accuracy of the LC maps extracted from Google Earth high-resolution images based on simple random sampling [55]. To this end, 173, 179, 177, and 166 ground control points selected to assess the accuracy of the classification for 1984, 2000, 2010, and 2017, respectively. Finally, Kappa coefficient was used to assess the accuracy of the produced LC maps [56,57]. The Kappa coefficient can range from −1 to +1, where 1 is perfect agreement, 0 is exactly what would be expected by chance, and negative values indicate agreement less than chance [58,59].

Land Cover Change Prediction
The land change modeler (LCM) as a part of the TerrSet software package applied to predict future LC maps for years 2027 and 2037. LCM is an integrated vertical tool to analyze LC changes and predict future LC trends [60]. Three main steps of LCM were (1) analyzing LC changes, (2) creating transition potential maps, and (3) predicting future LC maps. In the first step, LC changes were analyzed between the generated LC maps, comparing 1984 versus 2000, 2000 versus 2010, and 2010 versus 2017. In this step, LC transition maps created to analyze and determine spatial trends in LC changes as well as LC gains and losses. These maps and analyses were used to define future scenarios of LC change.
In the next step, transition potential maps were generated based on one or a group of sub-model(s) offered by LCM. In LCM, change analysis and prediction were performed by a series of empirically evaluated transition sub-models. Each sub-model simulates the conversion from one LC to another, e.g., forest to agriculture.
The difference between the sub-models is their transition potential map were defined based on scenarios of changes identified in the historical records [60]. In this study, transitions are modeled using a multi-layer perceptron (MLP) artificial neural network (ANN) [61,62]. The driving variables for each sub-model must be identified and introduced as its inputs. We identified a list of environmental variables recommended in the literature as potential drivers of LC changes [63,64]. These variables include elevation, slope, aspect, distance from roads, distance from rivers, distance from rangelands, distance from bare lands, and distance from agricultural areas. The explanatory power of these variables was tested by Cramer's V that should have a value higher than 0.15 for influential variables [65]. Finally, in the third step, the amount of LC change that expected to occur in the future was predicted using the Markov Chain prediction process and the best-calibrated sub-model. The Markov Chain determines the amount of changes using the earlier and later LC maps along with the date specified [66]. The predictions made for LC changes in 2027 and 2037 were based on multi-objective land allocation.

Details of Sub-Models
In order to predict the LC changes, each LC transition must be modeled empirically, called transition potential maps [67]. In other words, this is a past trend-based model such that the rate of change is obtained from the analysis of change, which occurred during a previous period (called calibration period). Also, a transition sub-model can consist of a single land cover transition or a group of transitions that are thought to have the same underlying driver variables [60]. We used two calibration periods: 1984-2000 and 2000-2010 and attempted to identify the most important transition sub-model maps for these periods. The large number of sub-models causes data processing to be slow and the low number of sub-models also leads to data loss. The choice of the optimal number of sub-models can be based on the size of the changes that have been occurred in a region. The LCM has made it possible to choose a threshold for choosing optimal number of sub-models and users can ignore the small changes [67]. Therefore, considering the extent of the study area, we ignored LC transitions that occurred in less than 700 km 2 (less than 3% of the study area). Based on this precondition, we obtained five sub-models for the first calibration period (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and six sub-models for the second calibration period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). In the first calibration period, the main transition sub-models include: (1) Forest to agriculture, (2) agriculture to rangeland, (3) bare land to rangeland, (4) rangeland to agriculture, and (5) all transmission sub-models. Furthermore, over the second calibration period, sub-models include: (1) Agriculture to forest, (2) agriculture to rangeland, (3) forest to agriculture, (4) rangeland to agriculture, (5) rangeland to bare land, and (6) all transmission sub-models. Tables 2 and 3 show the details of sub-models of the first and second calibration periods, respectively. We then used the identified sub-models to generate transmission potential maps, the so-called soft map that assigns a number from [0, 1] to each pixel to demonstrate its potential for the transition from one LC to another one. In the transition potential modeling process, sub-models can be defined by enabling or disabling identified transitions to determine modeling scenarios. These sub-models have been combined in different scenarios to obtain the best modeling results compared to historically observed LC changes. The sub-model (scenario) with the highest prediction accuracy based on the Kappa index for the LC map of 2017 has been selected to predict LC maps for 2027 and 2037.

LC Maps
The  Table 4 gives more details of classification accuracies for different LC classes. The most difficult LC class to distinguish over the years was agriculture that was very similar to rangeland and bare land areas. Residential areas scored a good classification result in all four periods with more than 86% classification accuracy because they distinguished by both supervised and on-screen digitize methods. The overall Kappa index was 87%, 87%, 85%, and 88% for 1984, 2000, 2010, and 2017, respectively.

Changes in Water Bodies
The LC analysis shows that the major changes in water bodies of ARB from 1984 to 2017 include changes in the surface elevation of lakes, changes in the depth and width of rivers, and development of new lakes created by new reservoir dams. More specifically, the satellite image processing of this study identified six important changes in water bodies of ARB over the past 33 years. The geographic location of these changes shown in Figure 4, and the type and extent of them shown in Figure 5. As shown in Figure 4 Table 5 for full details about these six major changes in water bodies of the study area.
the biggest construction of an aquaculture production with an area of 1403 hectares. Dams No. 5 and 6 are located in the eastern part of the basin and aim to supply water for agriculture and increase infiltration. See Table 5 for full details about these six major changes in water bodies of the study area.

LC Map Processing and Validation
The variables listed in Table 6 considered as the most important drivers of LC changes in the study area. To achieve the best-performing LCM, various sub-models developed using different combinations of these variables. LC changes observed from 1984 to 2000 were used to calibrate the LCM model. The calibrated model showed high accuracy in the validation phase by scoring a Kappa index of 89% and 86% to predict the LC map of the study area in 2010 and 2017, respectively ( Figure  8A,B). When LC changes from 2000 to 2010 are used to calibrate the model, its accuracy for predicting the LC map of 2017 dropped to 84%, (see Figure 8C). It is worth noting that rivers were classified with the lowest prediction accuracy in terms of the Kappa index. The reason for this suboptimal accuracy is attributed to the dynamic nature of rivers that changes them over a relatively short period.

LC Map Processing and Validation
The variables listed in Table 6 considered as the most important drivers of LC changes in the study area. To achieve the best-performing LCM, various sub-models developed using different combinations of these variables. LC changes observed from 1984 to 2000 were used to calibrate the LCM model. The calibrated model showed high accuracy in the validation phase by scoring a Kappa index of 89% and 86% to predict the LC map of the study area in 2010 and 2017, respectively ( Figure 8A,B). When LC changes from 2000 to 2010 are used to calibrate the model, its accuracy for predicting the LC map of 2017 dropped to 84%, (see Figure 8C). It is worth noting that rivers were classified with the lowest prediction accuracy in terms of the Kappa index. The reason for this suboptimal accuracy is attributed to the dynamic nature of rivers that changes them over a relatively short period.

Predicting Future LCs
The most accurately calibrated LCM has been used to predict the LC maps for the years 2027 and 2037. The areas of these six classes of LC classes for all years from the past to the future were calculated and are shown in Figure 9. Moreover, the predicted LC maps for 2027 and 2037 are shown in Figure 10.
According to Figure 9, the area of rangeland was expected to decrease by 957 km 2 and 1207 km 2 in 2027 and 2037, respectively. On the other side, the areas of forest, bare land, agriculture, and residential area are expected to increase by 935 km 2 and 949 km 2 , 1216 km 2 and 1463 km 2 , 230 km 2 and 305 km 2 , and 58 km 2 and 125 km 2 in 2027 and 2037, respectively.

Predicting Future LCs
The most accurately calibrated LCM has been used to predict the LC maps for the years 2027 and 2037. The areas of these six classes of LC classes for all years from the past to the future were calculated and are shown in Figure 9. Moreover, the predicted LC maps for 2027 and 2037 are shown in Figure 10.

Predicting Future LCs
The most accurately calibrated LCM has been used to predict the LC maps for the years 2027 and 2037. The areas of these six classes of LC classes for all years from the past to the future were calculated and are shown in Figure 9. Moreover, the predicted LC maps for 2027 and 2037 are shown in Figure 10.
According to Figure 9, the area of rangeland was expected to decrease by 957 km 2 and 1207 km 2 in 2027 and 2037, respectively. On the other side, the areas of forest, bare land, agriculture, and residential area are expected to increase by 935 km 2 and 949 km 2 , 1216 km 2 and 1463 km 2 , 230 km 2 and 305 km 2 , and 58 km 2 and 125 km 2 in 2027 and 2037, respectively.  According to Figure 9, the area of rangeland was expected to decrease by 957 km 2 and 1207 km 2 in 2027 and 2037, respectively. On the other side, the areas of forest, bare land, agriculture, and residential area are expected to increase by 935 km 2 and 949 km 2 , 1216 km 2 and 1463 km 2 , 230 km 2 and 305 km 2 , and 58 km 2 and 125 km 2 in 2027 and 2037, respectively.

Discussion
This research combined LCM and Markov chain prediction models to enable analyzing LC changes in both temporal and spatial dimensions. The model parameters were created based on the historical LC data of 1984, 2000, 2010, and 2017. Environmental drivers and different scenarios analyzed in order to understand the pattern of changes in LC in the area. The models were run several times in order to realize the optimal model and then LCM was used to predict changes in LC for the future. The results of the evaluation of transition potential modeling using ANN-MLP showed more than 80% accuracy in most scenarios. In the study by Pérez-Vega et al. [68], the transition potential modeling with ANN modeling was carried out. They processed three sub-model series, i.e., recovering sub-model (agricultural land to tropical deciduous forests and agricultural land to secondary tropical deciduous forests), deforestation sub-model (tropical deciduous forests to agricultural land, tropical deciduous forests to rangeland, secondary tropical deciduous forests to agricultural land, and secondary tropical forests to rangeland), and perturbation sub-model (tropical deciduous forests to secondary tropical deciduous forests). Finally, they achieved an accuracy of 59.2%, 35.2%, and 59.6%, respectively. As opposed to that study, here we achieved a higher accuracy. Our obtained high accuracy of transition potential modeling using ANN compared to the study by Pérez-Vega et al. [68] is due to some powerful significant qualitative variables in LC changes or the choice of more appropriate input variables in each sub-model in the present study. In another study, Khoi and Murayama [69] achieved similar results in modeling transition potential with an ANN (about 80% accuracy) in Tam Dao National Park. Furthermore, in a study by Oñate-Valdivieso and Bosque Sendra [70], more than 70% accuracy in all transitions was achieved using an ANN, which is similar to the results of the present study.
Since the ARB is a populated basin, the LC analysis was concentrated on anthropogenic impacts mainly on construction, residential areas, forest degradation, and agriculture. Land use change occur at the local scale, but its environmental influences spread across regional and global scales [71]. Expansion of agricultural land is the most pervasive anthropogenic land conversion process in Asia [72]. Based on the obtained results, the agriculture class has increased by 399,930 hectares from 1984 to 2017. Most agricultural lands currently in use occupied former primary forests, rangeland, and wetlands [73,74]. Agricultural activities have caused serious degradation of the soil in South and Southeast Asia [75].
In this study, substantial changes occurred in water bodies during the years 1984 to 2017, such as the decrease and increase in the water level of lakes and changes in depth and width of rivers and depositing and building new constructions. However, processing images showed six important changes in water bodies over the past 33 years, in particular the dam site. The rapid construction of dams in Asia has created wide-spreading loss or fragmentation of freshwater habitats, particularly those of riparian floodplains and wetlands, although dams have executed a significant role in water supply, flood control, irrigation, and hydroelectric power production [72].
According to the captured satellite images for these years 1984, 2000, 2010, and 2017, as well as extracting LC for these years, water bodies have some fluctuations in this period. Most changes took place in 2000 with a decrease in the area from 13,983 to 7120 hectares from 1984 to 2000. This decrease continued from 2010 to 2017. Overall, changes in water bodies in comparison to other LCs during the past 33 years are limited and on a small scale. These findings are consistent with a study by Hajihosseini et al. [76]. They quantified the causes of the decrease of streamflow in the transboundary Hirmand (Helmand) located between in in Iran and Afghanistan. Th same as the current study, they applied remote sensing to understand watershed response to environmental changes and land use change. Their finding showed that agricultural developments in the Helmand basin decreased the discharge with about 346 million cubic meters/year accompanying an increase of 64,000 ha in an irrigated area in the river. Their results also revealed that the impact of land use change increased significantly for more recent periods and causes a reduction of 810 million cubic meters in annual streamflow. In another similar study by Al-saady et al. [77], they investigated impact of land use expansion on the natural environment of the Little Zab River basin in Iraq. Using three indices, NDVI, NDWI, and difference built-up index (NDBI), the change detection of the river has been performed considering six main classes of land use land cover. Their findings showed that continuous population growth with increasing individual demand for water used for irrigation, urban, and industrial uses resulted in intense competition in the distribution of the watershed's scarce water resources.
At the same time, it is common knowledge that LC practices such as deforestation and regeneration of wetlands also affect the quality and wellbeing of freshwater bodies. Deforestation has contributed to river siltation and failing to restore the wetlands has reduced the ability of river basins to control flows, leading to significant flooding. Both of these actions modify freshwater environments causing species to decline or vanish in them [78,79]. Moreover, forests as one of the most important LC classes to the sustainability of a basin underwent a big decrease from 1984 to 2000. Accordingly, the area of forest class decreased from 495,370 hectares in the first period to 164,511 hectares in the second period, and the change affects~66.7 percent of the whole area of forest class in the study area. These results are in line with the results of Moreno-Sanchez et al. [80], who showed that the total area of forests in Armenia decreased from 294,135 in the period 1987-1989 to 246,098 in the period 2000-2001. Furthermore, regarding the transmission sub-models, almost all of these changes have converted to agriculture. This is in line with a general trend for the wider region. According to Houghton and Hackler [81], in most parts of Asia, forest biomass and soil carbon have been reduced by 58% and 18%, respectively. Deforestation and forest fragmentation also created significant threats to the flora and fauna that live inside [72]. Ogutu et al. [82] reported a decline in some wild animal species (e.g., giraffe, waterbuck, impala) in the pastoral ranches bordering the Maasai Mara National Reserve, largely due to gradual degradation of forests and land use change.
In this study, the residential areas have grown during the years 1984, 2000, 2010, and 2017 by 81,290, 85,370, 95,971, and 101,696 hectares, respectively. Residential development in the basin has grown, and considering the direct link between residential area development and water consumption, it implies that much more water consumed in upstream parts of the basin. LC practices such as deforestation and regeneration of wetlands also influence the nature and protection of freshwater bodies. Kiragu [83] reported high sedimentation in Mara River, which is linked to watershed deforestation and intensive agriculture currently practiced in east Africa. According to Kiragu [83] study, agricultural deforestation and intensification are likely to cause an increase in surface runoff due to watershed degradation, which decreases its ability to absorb rainwater (reduced infiltration).
In line with the recent trend, continuous expansion of the agricultural LC class is expected. Consequently, the growth process in agriculture will lead to increased water consumption and water demand. Meantime, the decreasing trend in rangeland class will decrease infiltration and resilience in the basin, and the volume of sedimentation will increase. Since the Aras River Basin is a transboundary basin, more water consumption could lead to environmental risk in the commonwealth, particularly in downstream areas like Iran. According to the data that were collected from the FAO [84], the amount of water consumption in three different types of LCs, including agricultural, residential, and industrial areas, over time was calculated. The population of Armenia was also collected from 1993 to 2017 based on the FAO [84] database. As shown in Figure 11, water consumption in agriculture has a significant increasing trend from 1993 to 2017, which is in good agreement by the extracted agriculture in Landsat images and its expansion; again, this trend can be concluded as industrial consumption. Also, in Armenia, land area devoted to cereals declined over the past decade, with an 18% decrease from 2001 to 2010, where land area devoted to fruits increased by 25%, and land area devoted to vegetables, including melons, increased by 21% from 2001 to 2010 [85], which indicates a change in cultivation pattern and the use of species with higher water requirements, which can be one of the reasons for increasing water harvest and consumption in agriculture. Although water consumption in agriculture and industrial areas has increased over time, based on the results, the amount of water consumption in the residential areas has decreased. This decrease can be due to a decreasing trend in population and changes encouraged by the Armenia government, such as modification of water consumption patterns. These modifications of patterns can include reducing bath time, using water-reducing appliances, identifying, repairing, or replacing damaged appliances, and so on. Also, this dramatic drop is attributed to the introduction of water metering and a volumetric billing system. from the FAO [84], the amount of water consumption in three different types of LCs, including agricultural, residential, and industrial areas, over time was calculated. The population of Armenia was also collected from 1993 to 2017 based on the FAO [84] database. As shown in Figure 11, water consumption in agriculture has a significant increasing trend from 1993 to 2017, which is in good agreement by the extracted agriculture in Landsat images and its expansion; again, this trend can be concluded as industrial consumption. Also, in Armenia, land area devoted to cereals declined over the past decade, with an 18% decrease from 2001 to 2010, where land area devoted to fruits increased by 25%, and land area devoted to vegetables, including melons, increased by 21% from 2001 to 2010 [85], which indicates a change in cultivation pattern and the use of species with higher water requirements, which can be one of the reasons for increasing water harvest and consumption in agriculture. Although water consumption in agriculture and industrial areas has increased over time, based on the results, the amount of water consumption in the residential areas has decreased. This decrease can be due to a decreasing trend in population and changes encouraged by the Armenia government, such as modification of water consumption patterns. These modifications of patterns can include reducing bath time, using water-reducing appliances, identifying, repairing, or replacing damaged appliances, and so on. Also, this dramatic drop is attributed to the introduction of water metering and a volumetric billing system. Finally, water availability in the Aras river basin is of shared interest among Armenia, Turkey, and Iran [86]. The downstream riparian countries including Armenia, Iran, Georgia, and Azerbaijan have identified the major transboundary issues in the Kura-Aras basin as being: Hydrological flow variation and reduction, water quality deterioration, river basin ecosystem degradation, and increased flooding and bank erosion [86].
Construction of dams and the expansion of irrigated farms in upstream areas would undoubtedly affect flow patterns, water availability, and freshwater ecosystems, such as wetlands and alluvial plains in the downstream states, mainly Iran and Azerbaijan [86,87].
The present study is indicative of the large influences that LC change has on water resources in the outlet of the case study as domestic end-users. Based on the data that were recorded by the Ministry of Energy in Iran, approaching 2017, the trend of the observed streamflow showed a remarkable reduction (see Figure 12), which can be attributed to increased water withdrawal because Finally, water availability in the Aras river basin is of shared interest among Armenia, Turkey, and Iran [86]. The downstream riparian countries including Armenia, Iran, Georgia, and Azerbaijan have identified the major transboundary issues in the Kura-Aras basin as being: Hydrological flow variation and reduction, water quality deterioration, river basin ecosystem degradation, and increased flooding and bank erosion [86].
Construction of dams and the expansion of irrigated farms in upstream areas would undoubtedly affect flow patterns, water availability, and freshwater ecosystems, such as wetlands and alluvial plains in the downstream states, mainly Iran and Azerbaijan [86,87].
The present study is indicative of the large influences that LC change has on water resources in the outlet of the case study as domestic end-users. Based on the data that were recorded by the Ministry of Energy in Iran, approaching 2017, the trend of the observed streamflow showed a remarkable reduction (see Figure 12), which can be attributed to increased water withdrawal because of LC change, especially agricultural activities. This finding is inconsistent with a study by Mwangi et al. [88], as they estimated that land use change in the last 50 years contributed to about 97% of the observed increase in mean streamflow of Nyangores River (a headwater tributary of the Mara River in east Africa).
Based on prediction maps for the years 2027 and 2037, as well as the results for the past, the governments of Armenia and the Nakhchivan Autonomous Republic are developing agriculture and residential areas and constructing dams. Consequently, according to the articles 5, 6, and 7 of the Convention on the Law of Non-Navigational Uses of International Watercourses Adopted by the General Assembly of the United Nations on 21 May 1997, "These activities should not violate the rights of stakeholder downstream countries".

Limitations of the Study
Finally, we also encountered some limitations when conducting the analysis. Different seasonal timing of Landsat images can dramatically influence the LC classification, especially in LC such as agriculture, rangeland, and bare land, which are strongly related to NDVI. In this respect, the most important limitations we faced were dealing with the incompleteness of the USGS-Archive as well as the extreme Land Cloud Cover of the study area during our chosen dates. In order to overcome this limitation in future studies, and so to better understand the conversion of rangeland to bare land, it is suggested to work with a more complete version of USGS-Archive. of LC change, especially agricultural activities. This finding is inconsistent with a study by Mwangi et al. [88], as they estimated that land use change in the last 50 years contributed to about 97% of the observed increase in mean streamflow of Nyangores River (a headwater tributary of the Mara River in east Africa). Based on prediction maps for the years 2027 and 2037, as well as the results for the past, the governments of Armenia and the Nakhchivan Autonomous Republic are developing agriculture and residential areas and constructing dams. Consequently, according to the articles 5, 6, and 7 of the Convention on the Law of Non-Navigational Uses of International Watercourses Adopted by the General Assembly of the United Nations on 21 May 1997, "These activities should not violate the rights of stakeholder downstream countries".

Limitations of the Study
Finally, we also encountered some limitations when conducting the analysis. Different seasonal timing of Landsat images can dramatically influence the LC classification, especially in LC such as agriculture, rangeland, and bare land, which are strongly related to NDVI. In this respect, the most important limitations we faced were dealing with the incompleteness of the USGS-Archive as well as the extreme Land Cloud Cover of the study area during our chosen dates. In order to overcome this limitation in future studies, and so to better understand the conversion of rangeland to bare land, it is suggested to work with a more complete version of USGS-Archive.

Conclusions and Outlook
LC changes and modifications of the Aras River in the territory of Armenia and the Nakhchivan Autonomous Republic over 33 years from 1984 to 2017 investigated to predict LC maps in the years 2027 and 2037. Environmental driver variables and different scenarios defined to understand the pattern of changes in LC in the Aras river basin, and the models were operated several times to realize the optimal model and then, using LCM, predict changes in LC for the future. The outcomes showed that the models are able to predict the LC changes in the study area with accuracy from 81% to 85% using the Kappa index. This study focused on anthropogenic impacts mainly on construction and agriculture in the basin, and based on the outcomes, six important dams and aquacultures built on the upstream of the basin over the past 33 years. On the other hand, agriculture and residential classes as the most important water consumers have increased around 399,930 and 20,406 hectares from 1984 to 2017, respectively. In addition, according to prediction maps, the increase in this trend continues

Conclusions and Outlook
LC changes and modifications of the Aras River in the territory of Armenia and the Nakhchivan Autonomous Republic over 33 years from 1984 to 2017 investigated to predict LC maps in the years 2027 and 2037. Environmental driver variables and different scenarios defined to understand the pattern of changes in LC in the Aras river basin, and the models were operated several times to realize the optimal model and then, using LCM, predict changes in LC for the future. The outcomes showed that the models are able to predict the LC changes in the study area with accuracy from 81% to 85% using the Kappa index. This study focused on anthropogenic impacts mainly on construction and agriculture in the basin, and based on the outcomes, six important dams and aquacultures built on the upstream of the basin over the past 33 years. On the other hand, agriculture and residential classes as the most important water consumers have increased around 399,930 and 20,406 hectares from 1984 to 2017, respectively. In addition, according to prediction maps, the increase in this trend continues in 2027 and 2037 for these classes. It should be regarded that water supply for farming and guarantying the security of agriculture will undoubtedly be a concern for the governments of Armenia and the Nakhchivan Autonomous Republic in the Aras River Basin. For these reasons, the governments have some simple options such as constructing a dam and seal and water diversion from the main channel for more irrigation and consumption. Overall, since the basin is a transboundary basin between Iran, Turkey, Armenia, Azerbaijan, and Georgia, and considering the importance of water in political relationships, a more balanced strategy for land covers and harvesting rivers should be applied. The international development experience of transboundary territories and the results of our study demonstrate that sustainable water resource management can be realized in specific territorial trends. Under poor water resource management, one frontier territory may, in general, have a negative impact on the management of a neighboring border territory. Comparative assessments of water resource management trends that are built on the neighboring countries' frontier territories are therefore required. At the same time, studying the long-term shifts in land use land cover change in some part of the world is one of the most critical tasks of development studies.
The study of all the complex problems affecting land use efficiency and ecological conditions in transboundary river basins, as well as the factors undermining their systemic structure, is therefore one of the primary objectives policy makers, which should be designed for all basins. The methods applied in this study can be used for long-term hydrological studies of other transboundary river basins in data-scarce regions like Iran.
Further research on the possible driving forces of deforestation and agricultural expansion found in this study is required to establish effective management strategies. For example, it may be necessary to determine the role of forest management institutions and socio-political drivers, and how and to what extent deforestation affects the livelihood of domestic end-users. It may also be necessary to examine how the impact of the political environment, and the income flow from Aras River Basin protection, influences the conversion of forest to agricultural land, and vice versa.