Water Balance Analysis Based on a Quantitative Evapotranspiration Inversion in the Nukus Irrigation Area, Lower Amu River Basin

: Human activities are mainly responsible for the Aral Sea crisis, and excessive farmland expansion and unreasonable irrigation regimes are the main manifestations. The conﬂicting needs of agricultural water consumption and ecological water demand of the Aral Sea are increasingly prominent. However, the quantitative relationship among the water balance elements in the oasis located in the lower reaches of the Amu Darya River Basin and their impact on the retreat of the Aral Sea remain unclear. Therefore, this study focused on the water consumption of the Nukus irrigation area in the delta of the Amu Darya River and analyzed the water balance variations and their impacts on the Aral Sea. The surface energy balance algorithm for land (SEBAL) was employed to retrieve daily and seasonal evapotranspiration (ET) levels from 1992 to 2018, and a water balance equation was established based on the results of a remote sensing evapotranspiration inversion. The results indicated that the actual evapotranspiration (ET a ) simulated by the SEBAL model matched the crop evapotranspiration (ET c ) calculated by the Penman–Monteith method well, and the correlation coe ﬃ cients between the two ET a sources were greater than 0.8. The total ET a levels in the growing seasons decreased from 1992 to 2005 and increased from 2005 to 2015, which is consistent with the changes in the cultivated land area and inﬂows from the Amu Darya River. In 2000, 2005 and 2010, the groundwater recharge volumes into the Aral Sea during the growing season were 6.74 × 10 9 m 3 , 1.56 × 10 9 m 3 and 8.40 × 10 9 m 3 ; respectively; in the dry year of 2012, regional ET exceeded the river inﬂow, and 2.36 × 10 9 m 3 of groundwater was extracted to supplement the shortage of irrigation water. There is a signiﬁcant two-year lag correlation between the groundwater level and the area of the southern Aral Sea. This study can provide useful information for water resources management in the Aral Sea region.


Introduction
The Aral Sea crisis has become a global ecological and social problem and is also one of the hotspots in water resources research for Central Asia. As an inland lake, the variations in the water volume of the Aral Sea are mainly determined by surface runoff, evaporation, rainfall and groundwater recharge. Since the 1960s, surface runoff has been greatly reduced, and the Aral Sea has experienced the rapid shrinkage, salinization and deterioration of its ecological environment. According to previous studies, the shrinking of the Aral Sea has resulted from climate change and human activities; however, human activities were the main factor, with natural climate variability being responsible for approximately 25% of the observed level drop and the remaining 75% due to human impacts [1][2][3]. Human activities mainly include unsustainable irrigation methods and the expansion of cultivated land areas around the Amu Darya and Syr Darya River Basins [4][5][6][7], which have led to the overuse and waste of surface water resources and thus have indirectly affected the Aral Sea. For example, nearly all drainage of the Amu Darya River, which is the largest river in central Asia, is presently controlled and used for irrigation [5]. Due to the construction of water conservation facilities and irrigation systems in the upper reaches, a number of immigrants were attracted to the middle and lower reaches to reclaim land, which caused the irrigation area to continue expanding, and agricultural water conditions were aggravated [8], which also influenced the water balance of the Aral Sea by perturbing surface runoff and groundwater recharge [9][10][11][12]. Therefore, it is important to study the water balance dynamics of the irrigation areas in the Aral Sea Basin when agricultural activities prevailed.
Several studies have examined the water balance elements in the irrigation areas of the Aral Sea Basin but have mostly focused on single elements such as evapotranspiration (ET) or runoff. For example, Singh et al. [13] used Gravity Recovery and Climate Experiment (GRACE) and MODIS Global Evapotranspiration Project (MOD16) data with a coarse resolution to analyze changes in water mass and ET in the lower Amu Darya River irrigation area. Nezlin et al. [11] and Asarin et al. [4] analyzed runoff from the Amu Darya River into the Aral Sea over decades. Conrad et al. [14] studied crop water consumption by using coarse-resolution MODIS data. Ochege et al. [15] used remote sensing technology to analyze the annual variations in ET. However, the relationship among river inflow, water utilization and transformation in irrigation areas was not analyzed from the perspective of water resources.
As an important part of water balance, surface runoff, precipitation, groundwater recharge and evapotranspiration have a great influence on regional agricultural production and water resource management. The characteristics of the surface runoff of the Amu Darya River determine the water inflows to the Aral Sea and the irrigation area in the Amu Darya delta [4,16]. Accurate measurement of the region's inflow is essential and can be obtained from the several hydrological stations which located along the main stream of Amu Darya River, and variations of the runoff at the inlet station of the Amu Darya delta showed a significant decreasing trend from 1960 to 2006 [4]. Precipitation is the basic part of the water cycle and one of the main parameters in the water balance equation. However, Uzbekistan is largely an arid region, where evaporation exceeds rainfall and annual precipitation is below 200 mm; this means that precipitation is not the main component of irrigation water in the Amu Darya delta [16,17]. In the lower Amu Darya River Basin, groundwater is a life-sustaining resource that supplies water to local people, plays a central part in irrigated agriculture and influences the health of many ecosystems [18]. Schettler et al. [19] found some isotope evidence that groundwater from cretaceous aquifers in the Amu Darya delta recharged the Aral Sea. But there is little research on the groundwater reserves in the Amu Darya delta and its influence on water balance in the irrigated area. Evapotranspiration includes vegetation transpiration and soil and water evaporation and is an important part of the surface energy balance and water balance as well as a key parameter for studying and simulating land surface processes [20,21]. Accurate estimations of ET are essential for analyzing the water balance in the Aral Sea region because ET is the main process for water dissipation in arid regions [22,23]. Traditional site and field ET measurements include the Eddy correlation method [24], Bowen ratio energy balance method [25], photosynthesis instruments [26], weighing lysimeters [27] Remote Sens. 2020, 12, 2317 3 of 25 and other methods, which can be used to accurately obtain ET at the point or field scales, but these methods are limited by their small spatial scales and higher costs and would be difficult to use in the inland basins of Central Asia, which have sparse observation networks [28]. The water balance can be used to calculate ET values under conditions in which the other hydrologic factors are clear [29]. When integrated with the Penman-Monteith method and the Kc coefficients, the water demand of irrigated crops can be estimated, but a fine and accurate classification of planting structures is required. Increasingly mature remote sensing technology and evapotranspiration research open a new avenue for measuring and modeling ET. However, due to the uncertainty of MOD16 in sparsely vegetated areas, ET is underestimated [13]. Although Priestley-Taylor Jet Propulsion Laboratory Model (PT-JPL) actual evapotranspiration (ET a ) and Global Land Data Assimilation System (GLDAS) data are more accurate than MOD16 data [30], the rough resolution is not sufficient to carefully characterize the ET performance of irrigation areas. In general, existing global remote sensing data products do not meet the requirements for water resources calculations in the Aral Sea region. The increasingly mature evapotranspiration inversion model based on remote sensing can solve this problem, and the surface energy balance algorithm for land (SEBAL) model has been widely used especially.
Therefore, the objective of this study is to estimate the water consumption and analyze water balance in the Nukus irrigation area, which is in the lower reaches of the Amu River Basin. The SEBAL model was applied to model ET a levels over the study area from 1992 to 2018, and the simulated ET a were statistically compared with ground-based ET a values for different land cover types in the study area. The spatial and temporal variations in ET a were analyzed based on years of results. Additionally, combined with rainfall, river inflow and groundwater observation data, a water balance model was developed for analyzing the relationships among the main elements of the regional water resource system. The results can provide useful data support for water resources management and ecological protection in the Aral Sea region.

Study Area
The Nukus region is a typical agricultural irrigation area in central Asia, is located in the northwestern part of Uzbekistan and is part of the Khorezm region; the Nukus region is connected with Turkmenistan to the south and the Aral Sea to the north (Lat. 42 • to 44 • N and Long. 58 • to 61 • E). The elevation of the study area is between 5 and 10 meters with a total area of 14,247 km 2 . The Amu Darya River crosses the Nukus irrigation area and empties into the Aral Sea, while its high sediment discharge has formed very dynamic deltas. Since the 1980s, only a small amount of river inflow to the Aral Sea has occurred, and most of the water is utilized for Nukus agricultural irrigation [4]. Annual precipitation in the region is less than 200 mm, which is concentrated in summer and has shown a decreasing trend in recent years; the seasonal and intraday temperature differences are significant, where the lowest temperatures are between −1.4 • C and 4 • C and the highest temperature can reach 30 • C [31]. The typical crop types are cotton, wheat, rice, corn and alfalfa, and flood irrigation with low efficiency is the main irrigation method. The main wild vegetation types are ammodendron and tamarisk, with Euphrates poplar and reeds near puddles. Farming and animal husbandry are important sources of local income.
As shown in Figure 1, there are three meteorological stations, namely Nukus, Cungrad and Chimboy, which have obtained continuous observation records since the 1960s. A hydrological station called Samambay is located in the inlet of the Nukus irrigation area. In addition, an observation site for measuring pan evaporation was built in Nukus city, which provided validation data for the 2019 growing season.

Data Availability
Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+) and Operational Land Imager (OLI) images for the Nukus irrigation area (path/row is 161,160/030) were downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/). Taking into account the quality of the images, the number of images available each year and the time intervals, a total of 144 images, including 42 TM images, 59 ETM images and 43 OLI images, were selected from the growing seasons (April to October) of 1992,1996,2000,2005,2010,2012,2015,2018 and 2019. The amount of cloud cover above the irrigation areas in all images was below 5%, and the preprocessing for radiometric calibration, Fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) atmospheric correction and geographic correction using control points was carried out on ENVI (The Environment for Visualizing Images). Figure 2 shows the dates of the 144 remote sensing images.

Data Availability
Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+) and Operational Land Imager (OLI) images for the Nukus irrigation area (path/row is 161,160/030) were downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/). Taking into account the quality of the images, the number of images available each year and the time intervals, a total of 144 images, including 42 TM images, 59 ETM images and 43 OLI images, were selected from the growing seasons (April to October) of 1992,1996,2000,2005,2010,2012,2015,2018 and 2019. The amount of cloud cover above the irrigation areas in all images was below 5%, and the preprocessing for radiometric calibration, Fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) atmospheric correction and geographic correction using control points was carried out on ENVI (The Environment for Visualizing Images). Figure 2 shows the dates of the 144 remote sensing images. Meteorological data from 1990 to 2019 at the Nukus, Cungrad and Chimboy stations were obtained from the National Oceanic and Atmospheric Administration (https://gis.ncdc.noaa.gov/maps/ncei/cdo/daily) and included average temperature, daily maximum temperature, daily minimum temperature, dew point temperature, average wind speed, average air pressure and precipitation. The observation data of pan evaporation in 2019 were obtained from the Karapakstan Branch of the Institute of Water Problem, Uzbekistan. Monthly discharges at the Samambay hydrological station were collected from 1990 to 2015, and monthly groundwater level data from 1999 to 2017 were collected from the Ministry of Agriculture and Water Resources (MAWR) of Uzbekistan.
The sampling points of the agricultural planting structure were taken from field surveys in September 2018 and March 2019 and consisted of a total of 446 sampling points, including 121 rice sample points, 162 wheat sample points, 161 cotton sample points and 2 bare land points, which were located near three weather stations; for the specific locations, the reader is referred to the verification section. The detailed data are listed in Table 1.  Meteorological data from 1990 to 2019 at the Nukus, Cungrad and Chimboy stations were obtained from the National Oceanic and Atmospheric Administration (https://gis.ncdc.noaa.gov/maps/ncei/cdo/ daily) and included average temperature, daily maximum temperature, daily minimum temperature, dew point temperature, average wind speed, average air pressure and precipitation. The observation data of pan evaporation in 2019 were obtained from the Karapakstan Branch of the Institute of Water Problem, Uzbekistan. Monthly discharges at the Samambay hydrological station were collected from 1990 to 2015, and monthly groundwater level data from 1999 to 2017 were collected from the Ministry of Agriculture and Water Resources (MAWR) of Uzbekistan.
The sampling points of the agricultural planting structure were taken from field surveys in September 2018 and March 2019 and consisted of a total of 446 sampling points, including 121 rice sample points, 162 wheat sample points, 161 cotton sample points and 2 bare land points, which were located near three weather stations; for the specific locations, the reader is referred to the verification section. The detailed data are listed in Table 1.

SEBAL Model
The SEBAL model was initially developed by Bastiaanssen et al. [32]; this model is based on the energy balance equation of the land surface, and the distribution of sensible heat flux is calculated by the cyclic recursive method, while the ET a of a region is obtained by the energy residual method [33]. Bhattarai et al. [34] devised a simple decision-tree-based classifier to find "cold" and "hot" pixels. DEM and slope data are used to adjust the parameters of the SEBAL model for mountainous terrain, and good simulation results have been obtained [35]. The SEBAL model does not need large amounts of auxiliary ground-based data and reaches high spatial resolutions in ET simulations; the Monin-Obukhov similarity theory can reduce the uncertainty of iterations. However, the SEBAL model still has some limitations such that the choice of "cold" and "hot" pixels has a great influence on the results, and due to the unfavorable timing of remote sensing images, it is difficult to obtain accurate simulated results for continuous time periods.
The basic principle of the SEBAL model is: where R n is the net radiant flux, G is the soil heat flux, H is the sensible heat flux and LET is the latent heat flux.
(1) The net radiation flux is the net energy obtained after shortwave radiation from the sun reaches the Earth's surface and is reflected by the surface and is exchanged with the longwave radiation of the atmosphere [32]. The calculation formula for net radiation flux is as follows: where R n is the surface net radiation (Wm −2 ), R S↓ is the short-wave solar radiation reaching the surface (Wm −2 ), R L↓ is the downgoing atmospheric longwave radiation (Wm −2 ) and R L↑ is the longwave radiation emitted from the Earth's surface (Wm −2 ). α and ε o represent the surface albedo and surface emissivity, respectively. (2) Soil heat flux refers to the heat exchange energy that enters the soil. According to Bastiaanssen's research, the soil heat flux is calculated by the following method [36]: Note that the above formula applies to vegetation; when the Normalized Difference Vegetation Index (NDVI) is less than 0.157, use G = 0.2R n according to experience, and for water bodies, G = 0.5R n .
(3) Sensible heat flux represents the process of energy transfer from the land surface to the atmosphere and takes place mainly by conduction and convection. Sensible heat flux can be computed as follows: In the formula, ρ a represents the air density, C p is the air specific heat, dT = T 1 − T 2 and represents the temperature difference between two heights (Z 1 and Z 2 ) and r ah represents the aerodynamic resistance to heat transport. However, dT and r ah are unknown variables but can be calculated using an iterative algorithm based on Monin-Obukhov similarity theory: where k, u * , d m , Z 1 , Z 2 , Z x , u x and Z om represent von Karman's constant, the friction velocity, zero plane displacement, reference heights (Z 1 , Z 2 , and Z x ), wind speed at height Z x and momentum transfer roughness, respectively. Z om can be calculated by using the NDVI: The calculation process of r ah and u * involves an iterative algorithm. Since the algorithm is complex, please refer to previous articles for more detail [37]. As mentioned above, dT can be calculated as: where a and b are regression coefficients and the SEBAL model assumes that there are "cold" and "hot" pixels in the image. In this case, a and b can be calculated by the following formula: "Cold" pixels represent cases for which the sensible heat flux of a pixel is 0 and can be found under conditions when the NDVI is greater than 0.3 and the temperature is lowest; "hot" pixels represent cases for which the latent heat flux is 0 and can be found where the NDVI is less than 0.157 and greater than 0 and the pixels have the highest temperature.
(4) Sensible heat flux can be calculated according to the energy balance formula: Remote Sens. 2020, 12, 2317 8 of 25 Only the instantaneous latent heat flux can be obtained above, and one day of ET a can be obtained through time extrapolation. By introducing the concept of the evaporative fraction (EF ins ), which is considered a constant throughout the day, EF ins can be estimated as follows: By using EF ins , the daily ET a (ET 24_SEBAL ) can be determined. λ stands for the latent heat of vaporization, and R n24 and G n24 represent the net radiation flux and the soil heat flux over one day, respectively; the calculation method can be referred to in Allen [38]. ET r−period represents the overall ET r from the observations for the studied period from m to n [15,37]. The seasonal ET a (ET 24−period ) is derived by cumulatively adding the daily ET a grids of the growing season (April to October) [15,39,40]. Additionally, for Landsat TM/ETM/OLI images, the surface albedo (a) is calculated according to different methods [37,41,42]: where α toa is the reflectance at the top of the atmosphere, α path radiance refers to the path radiance, τ sw is the atmospheric transmissivity and α 1~α7 correspond to the surface reflectance for bands 1~7 of the satellite sensor. The operation of the SEBAL model, including the preprocessing of the Landsat images and retrieval of surface parameters, is entirely based on Interactive Data Language (IDL).

The FAO Penman-Monteith Equation
The FAO Penman-Monteith equation is considered a universal standard to estimate ET O because it closely approximates the short-grass ET O [38,43] and is widely used around the world in the absence of measured data [44][45][46][47]. The calculation formula is as follows: where is the slope of the saturated vapor pressure curve, e s is the saturated vapor pressure, e a is the actual vapor pressure, r is the psychrometric constant and u 2 is the wind speed at a height of 2 m. The reference ET can be calculated using ET O Calculate software by using the maximum temperature, minimum temperature, average temperature, wind speed and relative humidity as input sources, and the missing weather data are replaced by default values in the software. The crop evapotranspiration (ET c ) is actual evapotranspiration, which differs distinctly from ET O , as the ground cover, crop canopy properties and aerodynamic resistance of crops are different from the properties of grass. The effects of the characteristics that distinguish field crops from grasses are integrated into Kc. Crop ET c is calculated by multiplying the reference ET by a crop coefficient [38]: where Kc is the crop coefficient at a specific growth stage. In this study, the Kc coefficients for rice, wheat and cotton in this research area are referenced in the relevant research [48], refer to Table 2 for detailed values; the international Food and Agriculture Organization (FAO) considers the K C coefficient of dry, bare land to be between 0.15 and 0.2, and a K C coefficient for bare land of 0.2 was used in this study.

Statistical Evaluation
Validation is one of the most challenging problems in satellite-based ET estimation research. However, no actual ground-measured ET data were available in the study area. This study attempted to assess the ET estimation accuracy via comparisons with ET c [38]. The agreement between the SEBAL model-derived ET a and ET c calculated from ground meteorological datasets was assessed using the Pearson correlation, root mean square error (RMSE), mean absolute error (MAE) and percent bias (PBIAS); these methods are frequently used to measure the differences between the values predicted by a model or an estimator and the values observed. In the following formulas, x represents the true value and y represents the approximate value.

Water Balance Analysis
The water balance explores the relationship between the amount of water available under natural conditions and the demand for water in the socioeconomic environment. According to the characteristics of incoming water and considering that water consumption in the Nukus irrigation area involves most of the inflow water from the Amu Darya River being used for agricultural activities with little or no surface runoff flowing outward, a water balance formula over the growing season is established based on the available data: where P is the precipitation in the Nukus irrigation area, R in f low is the inflow from the Amu Darya River, which is also the main source of irrigation water in the study area, G recharge is the groundwater recharge from the upstream region to the Nukus irrigation area, E is the total ETa, G leaving is the total groundwater recharge to the Aral Sea or elsewhere (it is an unknown variable need to be calculated) and G is the variation in the groundwater volume, which can be calculated from the groundwater level. P, R in f low , E and G are all available from the data sources mentioned above. However, G recharge data cannot be obtained at present. Here, it is assumed that the groundwater recharge of the Nukus irrigation area is 0 to explore the changes in the water balance of the irrigation area under closed groundwater conditions. The water balance formula can be used to calculate G leaving .

Lag Correlation Analysis
In the relationship between two time series (yt and xt), the series yt may be related to past lags of the x-series. The cross-correlation function (CCF) is helpful for identifying lags of the x-variable that might be useful predictors of yt, and τ is the lag time.
There is a correlation between the hydrological elements and the area of the Aral Sea, and the time may be delayed. For example, the groundwater needs the propagation time to recharge. To explore the lag correlation between the hydrological elements and the area of the Aral Sea, cross-correlation analysis is used to extract the lag time and correlation coefficient, which can be implemented by the CCF function in R language. For continuous functions, CCF functions are defined as follows:

Accuracy Assessment of the SEBAL Model
The estimation results are evaluated at 446 sample locations and compared with independently estimated ET c following the FAO Penman-Monteith equation [38].  Figure 4 shows the correlation analysis results of ET c and ET SEBAL for all sampling points. To be more specific, there were a total of 446 crop sampling points, and the simulated ET a at each point was modeled by SEBAL ordered in a chronological sequence, which correlated with ET c . The results show that the correlation coefficients (R) at most sample points range from 0.8 to 1, especially for wheat ( Figure 4). More statistics are summarized in Table 3. In most regions, the correlation is greater than 0.6, except for a minority of rice sampling points near the Chimboy meteorological station; in particular, the correlations of all wheat sampling points in the Cungrad region exceed 0.6, the MAE of each crop is within 1.43 mm and the percentage deviation is within 11.42%. These results demonstrate that SEBAL has good applicability in the study area.   Figure 4 shows the correlation analysis results of ETc and ETSEBAL for all sampling points. To be more specific, there were a total of 446 crop sampling points, and the simulated ETa at each point was modeled by SEBAL ordered in a chronological sequence, which correlated with ETc. The results show that the correlation coefficients (R) at most sample points range from 0.8 to 1, especially for wheat ( Figure 4). More statistics are summarized in Table 3. In most regions, the correlation is greater than 0.6, except for a minority of rice sampling points near the Chimboy meteorological station; in particular, the correlations of all wheat sampling points in the Cungrad region exceed 0.6, the MAE of each crop is within 1.43 mm and the percentage deviation is within 11.42%. These results demonstrate that SEBAL has good applicability in the study area.     The simulation results of the monthly ET SEBAL and ET c values also exhibit good consistency.

Annual Variations of the Simulated Evapotranspiration
Evapotranspiration was spatiotemporally variable over the study area due to differences in the moisture content of the surface and the energy exchange properties of the different land covers. Furthermore, the ETa was also spatiotemporally variable within the same land cover type due to variations within land cover types [49]. Figure 6 shows the SEBAL-estimated monthly ETa for two typical years. Generally, the moist areas over the southern farmland part of the study area had higher ETa than the northern areas. Water bodies, the reed swamp and well-irrigated and mid-stage crops exhibited the maximum ETa. On the other hand, the dry, bare land and sparsely vegetated parts of the plain had limited moisture, which resulted in low ETa in the study area. The monthly ETa clearly changed during the growing season. From April to August, the maximum increased month by month, and after August, the maximum ETa generally decreased. The spatial pattern of the monthly ETa varied and was dependent on land use type and climate conditions at the same time. The monthly ETa for cultivated land and bare land in April and May exhibited low values in the growing season, and the differences between the two were small. Between June and August, the monthly ETa values for cultivated and bare land were significantly different. In July 2012, the monthly average ETa of

Annual Variations of the Simulated Evapotranspiration
Evapotranspiration was spatiotemporally variable over the study area due to differences in the moisture content of the surface and the energy exchange properties of the different land covers. Furthermore, the ET a was also spatiotemporally variable within the same land cover type due to variations within land cover types [49]. Figure 6 shows the SEBAL-estimated monthly ET a for two typical years. Generally, the moist areas over the southern farmland part of the study area had higher ET a than the northern areas. Water bodies, the reed swamp and well-irrigated and mid-stage crops exhibited the maximum ET a . On the other hand, the dry, bare land and sparsely vegetated parts of the plain had limited moisture, which resulted in low ET a in the study area. The monthly ET a clearly changed during the growing season. From April to August, the maximum ET a increased month by month, and after August, the maximum ET a generally decreased. The spatial pattern of the monthly ET a varied and was dependent on land use type and climate conditions at the same time. The monthly ET a for cultivated land and bare land in April and May exhibited low values in the growing season, and the differences between the two were small. Between June and August, the monthly ET a values for cultivated and bare land were significantly different. In July 2012, the monthly average ET a of cultivated land in the southern irrigation area was approximately 209 mm, while the monthly average ET a for bare land in the northern irrigation area was 77.50 mm. The ET a for cultivated land was clearly greater than that for bare land until September and October, at which point the differences diminished again.    7 shows the regional monthly average ET a in the growing season. The results roughly formed a single peak, and the highest ET a values were found in July with a range of 187.79 to 231.46 mm. The ET a in June and August ranked second, and the monthly ET a values were in the ranges from 120.02 to 197.67 mm and from 137.83 to 195.46 mm, respectively. In spring, the monthly ET a in April and May ranged from 114.80 to 131.34 mm and from 121.35 to 186.33 mm, respectively; the main ET a consumption was by water bodies and by cultivated land near water bodies. In autumn, the harvest season lasted from September to October, and the monthly average ET a were 141.59 mm and 119.31 mm respectively; the values became smaller because of the truncated R n and low air temperatures. The total ETa in the growing season were between 900 mm and 1200 mm.

Interannual Variations of the Simulated Evapotranspiration
Figure 8a-g shows the spatial distributions of total ETa in the growing season from 1992 to 2018. Overall, the spatial ETa distributions over the years are generally the same, and local differences are caused by changes in planting structure and land use. The ETa of water bodies (mazarine blue-dark green block in Figure 8a-g) can reach 1100-1500 mm during the growing season, while the ETa of cultivated land are between 700 and 1100 mm. Figure 8h shows the annual change rate of ETa from 1992 to 2018, and it is easy to see that the ETa near most of the water bodies (red-orange block in the northern irrigation area in Figure 8h) showed a decreasing trend, while the ETa of cultivated land in the southeast and southwest parts exhibited an increasing trend. These observations indicate that there were obvious expansions or abandonments of farmland; meanwhile, the area of reservoirs and lakes that met the needs of agricultural production in the Nukus irrigation area was shrinking.

Interannual Variations of the Simulated Evapotranspiration
Figure 8a-g shows the spatial distributions of total ET a in the growing season from 1992 to 2018. Overall, the spatial ET a distributions over the years are generally the same, and local differences are caused by changes in planting structure and land use. The ET a of water bodies (mazarine blue-dark green block in Figure 8a-g) can reach 1100-1500 mm during the growing season, while the ET a of cultivated land are between 700 and 1100 mm. Figure 8h shows the annual change rate of ET a from 1992 to 2018, and it is easy to see that the ET a near most of the water bodies (red-orange block in the northern irrigation area in Figure 8h) showed a decreasing trend, while the ETa of cultivated land in the southeast and southwest parts exhibited an increasing trend. These observations indicate that there were obvious expansions or abandonments of farmland; meanwhile, the area of reservoirs and lakes that met the needs of agricultural production in the Nukus irrigation area was shrinking.
According to the 1992-2015 ESA LUCC data, the average daily ET a trends for cultivated land, forest and bare land were generally consistent, but there were some differences (Figure 9). In general, the ET a of cultivated land and forestland were higher than that of bare soil. Evapotranspiration in cultivated land was higher, and most values were between 5 mm and 8 mm, whereas that in bare land was between 0 mm and 4.5 mm. cultivated land are between 700 and 1100 mm. Figure 8h shows the annual change rate of ETa from 1992 to 2018, and it is easy to see that the ETa near most of the water bodies (red-orange block in the northern irrigation area in Figure 8h) showed a decreasing trend, while the ETa of cultivated land in the southeast and southwest parts exhibited an increasing trend. These observations indicate that there were obvious expansions or abandonments of farmland; meanwhile, the area of reservoirs and lakes that met the needs of agricultural production in the Nukus irrigation area was shrinking.  According to the 1992-2015 ESA LUCC data, the average daily ETa trends for cultivated land, forest and bare land were generally consistent, but there were some differences (Figure 9). In general, the ETa of cultivated land and forestland were higher than that of bare soil. Evapotranspiration in cultivated land was higher, and most values were between 5 mm and 8 mm, whereas that in bare land was between 0 mm and 4.5 mm.  Figure 10 shows that the changes in ETa in the study area mainly exhibited two periods. One period was from 1992 to 2005, in which the ETa showed a decreasing trend, and the total ETa in the growing season (May to September) decreased from 960 mm to 690 mm; the second period was from 2005 to 2015, in which the seasonal ETa in the irrigation areas increased from 690 mm to 780 mm.   Figure 10 shows that the changes in ET a in the study area mainly exhibited two periods. One period was from 1992 to 2005, in which the ET a showed a decreasing trend, and the total ET a in the growing season (May to September) decreased from 960 mm to 690 mm; the second period was from 2005 to 2015, in which the seasonal ET a in the irrigation areas increased from 690 mm to 780 mm. Figure 10 shows that the changes in ETa in the study area mainly exhibited two periods. One period was from 1992 to 2005, in which the ETa showed a decreasing trend, and the total ETa in the growing season (May to September) decreased from 960 mm to 690 mm; the second period was from 2005 to 2015, in which the seasonal ETa in the irrigation areas increased from 690 mm to 780 mm.

Water Balance Analysis
The main factors for the water balance in Nukus irrigation area are the inflow from the AmuDarya River, evapotranspiration, precipitation and groundwater inflow and outflow. Figure 11a shows the changes in the river inflow and underground water quantities from 1999 to 2015. The results indicate that the inflows from the inlet station fluctuated during 1999 and 2017, as is

Water Balance Analysis
The main factors for the water balance in Nukus irrigation area are the inflow from the AmuDarya River, evapotranspiration, precipitation and groundwater inflow and outflow. Figure 11a shows the changes in the river inflow and underground water quantities from 1999 to 2015. The results indicate that the inflows from the inlet station fluctuated during 1999 and 2017, as is observable in the moving average curve. There were obvious dry and wet years. For example, in 2003, 2005, 2010 and 2012, which were typical wet years, the inflow water amounts were greater than 1.00 × 10 10 m 3 per year, while the inflow in dry years, such as 2001, 2008 and 2011, was less than 2.00 × 10 9 m 3 per year. The changes in groundwater also fluctuated with the river inflows. When inflow water was abundant, the groundwater increased, and groundwater decreased with low inflows.
Due to the availability of evapotranspiration data (SEBAL modeled), rainfall data, groundwater data and runoff data, the years 2000, 2005, 2010 and 2012 were selected as typical years for a regional water balance analysis (Figure 11b). The results indicate that ET and surface runoff were the two major components in the water balance of the irrigation area, and the rainfall values were small. In 2005, the river inflow was 1.11 × 10 10 m 3 , ET consumption was 9.87×10 9 m 3 , the groundwater volume increased by 1.20 × 10 8 m 3 , and ET consumption accounted for most of the inflow water. In 2010, river inflow was 1.31 × 10 10 m 3 , ET consumption was 9.79 × 10 9 m 3 and groundwater increased by 2.56 × 10 9 m 3 . In 2012, the inflow water was 7.84 × 10 9 m 3 , ET consumption was 1.31 × 10 10 m 3 , the groundwater decreased by 2.66 × 10 9 m 3 and the ET was greater than the river inflow amount, which resulted in decreased groundwater levels.
2005, the river inflow was 1.11 × 10 10 m 3 , ET consumption was 9.87×10 9 m 3 , the groundwater volume increased by 1.20 × 10 8 m 3 , and ET consumption accounted for most of the inflow water. In 2010, river inflow was 1.31 × 10 10 m 3 , ET consumption was 9.79 × 10 9 m 3 and groundwater increased by 2.56 × 10 9 m 3 . In 2012, the inflow water was 7.84 × 10 9 m 3 , ET consumption was 1.31 × 10 10 m 3 , the groundwater decreased by 2.66 × 10 9 m 3 and the ET was greater than the river inflow amount, which resulted in decreased groundwater levels.  Table 4 shows the calculated groundwater loss volumes for the Nukus irrigation area in the growing seasons of 2000, 2005, 2010 and 2012. Generally, ET was lower than the amount of incoming water, which came from precipitation and river inflow in wet years, so the surface runoff could meet  Table 4 shows the calculated groundwater loss volumes for the Nukus irrigation area in the growing seasons of 2000, 2005, 2010 and 2012. Generally, ET was lower than the amount of incoming water, which came from precipitation and river inflow in wet years, so the surface runoff could meet the water consumption requirements in the study area, and a part of the surface water converted to underground runoff through soil infiltration; the groundwater could recharge to the Aral Sea. In the wet years of 2005 and 2010, the groundwater losses in 2005 and 2010 were 1.56 × 10 9 m 3 and 8.40 × 10 8 m 3 , respectively. However, in the dry years, ET exceeded the amount of incoming water, and groundwater was extracted for irrigation. In 2000, the Nukus irrigation area had high groundwater levels, which not only satisfied the regional crop irrigation requirements but also replenished the Aral Sea by 6.74 × 10 9 m 3 . In 2012, the groundwater levels were low, the surface water and groundwater in the study area could not meet the irrigation requirements and underground runoff from the upstream regions of the Nukus was needed for supplementation.

Evaluation the Accuracy of the ET Products Modeled by SEBAL
According to previous studies, the ET a results obtained by the SEBAL model in highly vegetated areas will be slightly higher than those calculated using the Penman-Monteith formula by 5% and slightly lower for water bodies [50][51][52]. To eliminate this error as much as possible, the following improvements were made to the regular SEBAL model. Considering that the surface and atmospheric temperatures decrease with increasing altitude, it was necessary to adjust the surface temperature in the entire research area to the same reference height [53]. When calculating some parameters of the regular SEBAL model, a single empirical formula is adopted for the whole region, which is insufficient to reflect the complexity of the underlying surface. In this paper, different empirical formulas were used to calculate parameters such as soil heat flux and emissivity for different underlying surfaces [35,54,55]. To reduce the influence of clouds on the accuracy of "cold" pixel selection, higher NDVIs with the lowest surface temperatures were selected as the "cold" pixels. After adjustments, the simulation results agreed well with the ET c results. Ideally, a greater number of measured ET should have been used in the assessment but were not available; only the observed evaporation levels in 2019 from a 20 cm pan were obtained at the Nukus station. However, there was some deviation between the evaporation of the water surface observed with an evaporating pan and the actual water evaporation, and more experiments are needed to eliminate these biases; the verification accuracy from using the evaporation pan was not high. Generally, the ET sebal has a high fitting degree with ET c in this study. However, the fitting degrees were different among the main crop types, and higher RMSE values over 2 mm were found in the evapotranspiration estimation of the cotton. The source of error may contain the random error caused by crop sampling points and the impact from different soil backgrounds. The sensible heat flux and the soil heat flux are affected by soil adjusted vegetation index (SAVI), referring to the Allen et al. [37] study, and the soil impact factor (L) in SAVI was defined as 0.1, which means SEBAL model is insensitive to low vegetation coverage. The rice evapotranspiration is close to water evaporation due to its high water consumption, affected by field water, and the soil background has little disturbance to it. The soil also has less effect on wheat, which has a higher NDVI in early spring. As for the cotton, between April and October, it goes through the process of germination, growth, ripening and picking, and its NDVI shows the change of unimodal shape. The soil background has greater influence on the inversion results of cotton evapotranspiration. From the perspective of regional water balance, the error of crop evapotranspiration inversion would increase the uncertainty in regional ET calculation and water balance analysis. Future studies will improve the inversion accuracy by combining fine crop classification and field-scale evapotranspiration observation data.
With reference to previous studies, the SEBAL-modeled seasonal ET a were over 1200 mm from April to October in 2004 in Khorezm, which is near the Nukus irrigation area [14], the seasonal ET a of another Area Sea basin irrigation area to the north of the Nukus irrigation area were as high as 850 mm and the ET a of water bodies were as high as 1135 mm in 2012; the daily ET a values of the oasis-irrigation area in the growing season were within the range from 4 mm to 8 mm [15], while evaporation from an oasis agriculture evaporation pan could be as high as 1546 mm in Central Asia [56]. Annual evaporation from the Aral Sea as calculated by water balance was approximately 900-1100 mm [7]. In the 1990s, the observed evaporation of the Aral Sea actually reached 1220 mm [28]. Létolle et al. [57] predicted the evaporation of the Aral Sea from 2000 to 2050 with an average annual increase of 1.5 m, which could reach 1500 mm. In conclusion, compared with previous studies, the seasonal ET a results are within an acceptable range, and detailed information can be found in Table 5.

The Impact of LUCC on ET Variations
In this study, the results of seasonal ET a from 1992 to 2015 show a trend of decreasing first and then increasing; this phenomenon can be explained by LUCC. Table 6 presents the area changes of the main LUCC types in the Nukus irrigation area in 1992, 2000, 2010 and 2015. Cultivated land, grassland and bare land were the main land cover types in the Nukus irrigation area, which accounted for approximately 90% of the total area. Due to the characteristic high water consumption and large area of cultivated land, the ET of cultivated land is the main body of the total ET in the research area. Since the 1960s, the land reclamation movement throughout the Soviet Union led to desertification; soil fertility declined, and large areas of natural vegetation were destroyed. After the collapse of the Soviet Union, farmland irrigation technology was not mature, and flood irrigation was the main irrigation method, which caused the surface area of the Aral Sea to shrink significantly; when people realized the seriousness of the rapid shrinkage of the Aral Sea, they began to abandon the cultivated land. The cultivated land area in the Nukus irrigation area decreased from 1992 to 2010 and increased from 2010 to 2015. Other studies showed that the cultivated land in Uzbekistan showed a slight decline from 1993 to 2003 [16,58].
To further explore the relationship between changes in cultivated land area and ET changes, Figure 12 shows the change trend of the total ET and total ET of cultivated land. It is found that the change trend of cultivated land ET was the same as that of the cultivated land area and that total ET and ET of cultivated land accounted for 55% of total ET. All evidence indicates that the changes in the cultivated land area are the main reason for regional ET variations, especially in the irrigation areas where is dominated by highly water-consuming crops such as cotton and rice.
To further explore the relationship between changes in cultivated land area and ET changes, Figure 12 shows the change trend of the total ET and total ET of cultivated land. It is found that the change trend of cultivated land ET was the same as that of the cultivated land area and that total ET and ET of cultivated land accounted for 55% of total ET. All evidence indicates that the changes in the cultivated land area are the main reason for regional ET variations, especially in the irrigation areas where is dominated by highly water-consuming crops such as cotton and rice.

Exploring the Influences of the Nukus Water Balance on the Aral Sea
The hydrological factors of the Amu Darya delta are closely related to the changes in the Aral Sea. Previous studies have shown that the high groundwater levels in the lower Amu Darya irrigation area caused the groundwater to recharge to the Aral Sea [10,19]. Moreover, human activities and climate change also affected the Aral Sea volume. Human activities include the excessive reclamation of cultivated land and the construction of irrigation channels and various hydraulic facilities [16,59,60], which resulted in increased irrigation water consumption and decreased surface inflow to the Aral Sea. Climate change mainly affects temperature and precipitation and thus indirectly affects runoff. Several studies have shown that in the past few decades, rising temperatures and decreasing precipitation have caused the source of Amu River runoff to exhibit a decreasing trend [60][61][62]. Here, lag cross-correlation analysis can more clearly provide the connections between the hydrological factors of the Nukus irrigation area with the Aral Sea area. Figure 13 shows the lag cross-correlation analysis of the inflow volume to the Nukus irrigation area (a), groundwater level (b), annual precipitation (c) and annual average temperature (d) with the area of the southern Aral Sea. Only groundwater levels exhibit dominant cross-correlations with the Aral Sea area, while the other hydrological factors do not exhibit strong lag correlations. This is because Nukus is the nearest irrigation area to the southern Aral Sea with minimal surface runoff flows into the Aral Sea in recent years, while groundwater recharge is the main direct supply source to the southern Aral Sea. Based on the 1999-2017 data, the dominant cross-correlations occur when Lag = 2, the groundwater levels show a significant positive correlation with b and the correlation coefficient is 0.44. These results indicate that rising groundwater levels in the Nukus irrigation area will lead to expansion of the Aral Sea with a two-year delay. Two years may be the time needed for

Exploring the Influences of the Nukus Water Balance on the Aral Sea
The hydrological factors of the Amu Darya delta are closely related to the changes in the Aral Sea. Previous studies have shown that the high groundwater levels in the lower Amu Darya irrigation area caused the groundwater to recharge to the Aral Sea [10,19]. Moreover, human activities and climate change also affected the Aral Sea volume. Human activities include the excessive reclamation of cultivated land and the construction of irrigation channels and various hydraulic facilities [16,59,60], which resulted in increased irrigation water consumption and decreased surface inflow to the Aral Sea. Climate change mainly affects temperature and precipitation and thus indirectly affects runoff. Several studies have shown that in the past few decades, rising temperatures and decreasing precipitation have caused the source of Amu River runoff to exhibit a decreasing trend [60][61][62]. Here, lag cross-correlation analysis can more clearly provide the connections between the hydrological factors of the Nukus irrigation area with the Aral Sea area. Figure 13 shows the lag cross-correlation analysis of the inflow volume to the Nukus irrigation area (a), groundwater level (b), annual precipitation (c) and annual average temperature (d) with the area of the southern Aral Sea. Only groundwater levels exhibit dominant cross-correlations with the Aral Sea area, while the other hydrological factors do not exhibit strong lag correlations. This is because Nukus is the nearest irrigation area to the southern Aral Sea with minimal surface runoff flows into the Aral Sea in recent years, while groundwater recharge is the main direct supply source to the southern Aral Sea. Based on the 1999-2017 data, the dominant cross-correlations occur when Lag = 2, the groundwater levels show a significant positive correlation with b and the correlation coefficient is 0.44. These results indicate that rising groundwater levels in the Nukus irrigation area will lead to expansion of the Aral Sea with a two-year delay. Two years may be the time needed for groundwater to flow from the Nukus irrigation area to the southern Aral Sea, although further research is needed to prove this hypothesis.
Remote Sens. 2020, 12, 2317 20 of 23 groundwater to flow from the Nukus irrigation area to the southern Aral Sea, although further research is needed to prove this hypothesis.

Conclusions
In this study, the SEBAL model was developed to simulate the in the Nukus irrigation area in the lower reaches of the Amu Darya River Basin, and station-based calculated by the FAO Penman-Monteith method were used to evaluate the accuracy of the simulated results. Based on multiple years of results, the spatial and temporal variations in were analyzed. Moreover, a water balance model of the study area was established by coupling the estimated by SEBAL to observed inflow, precipitation and groundwater data. The ETSEBAL values were more consistent with the ETc values at daily and monthly scales, and the correlation coefficients (R) at most sample points were above 0.8. The exhibited spatiotemporal variations for different land cover types and varied within specific land cover types; the daily of cultivated land were approximately 5 mm to 8 mm, and the of bare land varied from 0 mm to 4.5 mm. Due to the uneven vegetation coverage, the in irrigation areas showed spatial heterogeneity such that the southern irrigation area had much higher than the northern area because of the high vegetation coverage of the former. Due to the effects of changes in the cultivated land area and the river inflow from the Amu Darya River, the accumulated in the growing seasons decreased from 1990 to 2010 and increased from 2010 to 2018. The seasonal in the Nukus irrigation area ranged from 900 mm to 1200 mm

Conclusions
In this study, the SEBAL model was developed to simulate the ET a in the Nukus irrigation area in the lower reaches of the Amu Darya River Basin, and station-based ET c calculated by the FAO Penman-Monteith method were used to evaluate the accuracy of the simulated ET a results. Based on multiple years of results, the spatial and temporal variations in ET a were analyzed. Moreover, a water balance model of the study area was established by coupling the ET a estimated by SEBAL to observed inflow, precipitation and groundwater data. The ET SEBAL values were more consistent with the ET c values at daily and monthly scales, and the correlation coefficients (R) at most sample points were above 0.8. The ET a exhibited spatiotemporal variations for different land cover types and varied within specific land cover types; the daily ET a of cultivated land were approximately 5 mm to 8 mm, and the ET a of bare land varied from 0 mm to 4.5 mm. Due to the uneven vegetation coverage, the ET a in irrigation areas showed spatial heterogeneity such that the southern irrigation area had much higher ET a than the northern area because of the high vegetation coverage of the former. Due to the effects of changes in the cultivated land area and the river inflow from the Amu Darya River, the accumulated ET a in the growing seasons decreased from 1990 to 2010 and increased from 2010 to 2018. The seasonal ET a in the Nukus irrigation area ranged from 900 mm to 1200 mm during the past 30 years. In wet years, irrigation water can replenish the groundwater and thus indirectly feed the southern Aral Sea; in dry years, groundwater meets the needs of crop evaporation, but even upstream groundwater needs to recharge. In 2000, 2005 and 2010, the groundwater supply quantities to the southern Aral Sea from the Nukus irrigation area were 6.74 × 10 9 m 3 , 1.56 × 10 9 m 3 and 8.40 × 10 8 m 3 , respectively, and these supplies have declined in recent years. In addition, there is a significant positive correlation between the groundwater level and southern Aral Sea area with a two-year lag. The SEBAL modeled ET product in this study can provide scientific data for estimating regional crop water consumption, and combined with the crop classifications based on remote sensing, irrigation water requirement at different scales with higher accuracy can be obtained. These should be the basis for irrigation water allocation. Moreover, based on the results of water balance analysis, a surface water and groundwater joint regulation strategy to minimize the irrigation water resources waste can be developed.