Abstract A WATER INFORMATICS APPROACH TO EXPLORING THE HYDROLOGICAL SYSTEMS OF BASINS WITH LIMITED INFORMATION; THE CASE OF THE BUSTILLOS LAGOON, CHIHUAHUA, MEXICO. BY Hugo Luis Rojas Villalobos, M.E.E. A dissertation submitted to the Graduate School in partial fulfillment of the requirements for the degree Doctor of Philosophy Major: Water Science and Management Concentration: Water Informatics NEW MEXICO STATE UNIVERSITY LAS CRUCES, NEW MEXICO October 2019 ProQuest Number: All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent on the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. Published by ProQuest LLC ( ProQuest ). Copyright of the Dissertation is held by the Author. All Rights Reserved. This work is protected against unauthorized copying under Title 17, United States Code Microform Edition © ProQuest LLC. ProQuest LLC 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, MI 48106 - 1346 27663222 27663222 2020 ii iii ACKNOWLEDGEMENTS Foremost, I want to express my eternal gratitude to Dr. Christopher Brown for his unconditional support and encouragement: "la vida es bella." I want to thank my doctoral committee for valuable support, time, and dedication. Profe Brown, Dr. Samani, Dr. Stringam, and Dr. Alatorre: you are responsible for who I am now :). Dr. Alfredo Granados, thanks for encouraging me, you are also responsible for this. To my children Luis Ricardo and Leonardo, if I did not have you, my life would be a disaster. Every time I see you, I realize that I am facing the living portraits of your mother and me, you are my main motivation, my motor of life. Olaya, thank you very much for your patience, dedication, and love for this family. To my parents and brother Hugo, Licha, and Luis, who always believed in me, especially my mother, who already rests in eternity, I love you. I thank Rojas and Villalobos families as well as the Aragones family for your constant support. To my friends and coworkers, you are the “neta” of the planet! Finally, to all members of the Academic department "Geoinformática Aplicada a Procesos Geo- ambientales" (UACJ-CA-94). Thanks to NMWRRI and the NMSU Department of Geography for your support. Part of this material is supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture, under award number 2015-68007-23130. iv VITA Educational Information December 16, 1973 Born in Mazatlán, Sinaloa, México 1992 Graduated from CBTis #51 High School Mazatlán, Sinaloa, México 1996 Graduated from Durango Institute of Technology Computer Systems Engineer Durango, Durango, México 2003 Specialization from Kanazawa Institute of Technology Systems designer. Microcontrollers. Kanazawa, Japan 2008 Graduated from Autonomous University of Ciudad Juarez Master of Environmental Engineer Ciudad Juárez, Chihuahua, México Employment Information 1996-2001 Operations Manager Principal Financial Group 2001-2009 Geographic Information Center Autonomous University of Ciudad Juarez 2009 - Full time professor - Geoinformatics Autonomous University of Ciudad Juarez Field of Study Major Field: Water Science and Management Concentration: Water Informatics v ABSTRACT A WATER INFORMATICS APPROACH TO EXPLORING THE HYDROLOGICAL SYSTEMS OF BASINS WITH LIMITED INFORMATION; THE CASE OF THE BUSTILLOS LAGOON, CHIHUAHUA, MEXICO By Student Name BY Hugo Luis Rojas Villalobos, M.E.E. Doctor of Philosophy NEW MEXICO STATE UNIVERSITY LAS CRUCES, NEW MEXICO October 2019 Christopher Brown, PhD. The analysis of hydrological basins requires information that is often not available or non-existent when the study areas are far from large urban centers. In the case of Bustillos Lagoon in the Mexican state of Chihuahua, hydrological information is limited, and government agencies do not share data with interested persons and research institutions. Given this barrier, this research contributes to filling information gaps concerning the geometry of the Bustillos Lagoon, evaporation, and morphometric parameters through the use of current technology in remote sensors, geographic information systems, and programming techniques that are used to extract, transform and process information. Chapter 2 deals with a new methodology that vi generates a 3D model of the bottom of the lagoon, which uses high-precision GPS surveys, bathymetry, regional digital terrain models, and satellite image time series. The analysis using the Kappa coefficient demonstrates that the overall performance of the 3D model is more significant than 0.89, which means that the model has a very high level of agreement. The analysis also showed that at greater depth, the agreement between the coverage of the water surface of the model and the images is relatively low (0.89), and this is due to the spatial resolution of the satellite images and strip banding errors of Landsat ETM +. On the other hand, on the upper level, there is an agreement close to 0.99 of the Kappa coefficients. Chapter 3 presents a performance comparison between the Regional Evapotranspiration Estimate Model (REEM) and the Earth Engine Evapotranspiration Flux (EEFlux) model, which are evapotranspiration models based on energy balances. These models can estimate the evaporation of water bodies. After applying statistical analysis, REEM performed better than EEFlux in quantifying the evaporation of the Bustillos Lagoon. Chapter 4 proposes an iterative algorithm to calculate morphometric variables (volume-area-height) using 3D models of water bodies. The implementation of the algorithm in the Python programming language showed that it is not necessary to develop complex equations that interrelate the morphometric variables, which by their nature, lead to more considerable uncertainty from the data source for their construction. This research document highlights the importance of cumulative multi-faceted knowledge to support and respond to regional water issues. Keywords: 3D model, topobathymetry, remote sensing, evaporation, iterative algorithm, lagoon geomorphometry 7 Table of contents Table of contents List of tables ...................................................................................................................................10 List of figures .................................................................................................................................11 Chapter 1 ........................................................................................................................................13 Introduction ................................................................................................................................13 Research location .......................................................................................................................15 Chapter 2 Topobathymetric 3D model reconstruction of shallow water bodies through remote sensing, GPS, and bathymetry .......................................................................................................17 Abstract ......................................................................................................................................17 Introduction ................................................................................................................................18 Materials and methods ...............................................................................................................20 Bathymetry .............................................................................................................................22 Contour extraction from remote sensors ................................................................................24 Topography ............................................................................................................................27 Digital elevation model ..........................................................................................................27 Topobathymetric 3D model and volume estimation ..............................................................27 Statistical Evaluation .............................................................................................................28 Results and discussion ...............................................................................................................31 8 Conclusions ................................................................................................................................38 References ..................................................................................................................................41 Chapter 3 Comparison of evaporation estimates from REEM and EEFlux models in a shallow water body. Case: Bustillos Lagoon, Chihuahua, Mexico. ............................................................44 Abstract ......................................................................................................................................44 Background ................................................................................................................................50 Reference evapotranspiration model ......................................................................................50 Brief description of remote sensing models ...........................................................................50 REEM ....................................................................................................................................51 METRIC ................................................................................................................................53 EEFlux ...................................................................................................................................55 Material and methods .................................................................................................................55 Agro-meteorological data ......................................................................................................58 Landsat 8 OLI selection .........................................................................................................58 REEM and EEFlux raster .......................................................................................................59 Lagoon delineation .................................................................................................................59 Statistical evaluation ..............................................................................................................59 Results ........................................................................................................................................61 Discussion ..................................................................................................................................67 Conclusions ................................................................................................................................70 References ..................................................................................................................................72 9 Chapter 4 Single-input, multiple-output iterative algorithm for the calculation of volume, area, elevation, and shape using 3D topobathymetric models. ...............................................................78 Abstract ......................................................................................................................................78 Introduction ................................................................................................................................79 Study area ...............................................................................................................................81 Material and methods .................................................................................................................82 Results and discussion ...............................................................................................................87 Recommendations ......................................................................................................................90 Conclusions ....................................................................................................................................94 References ......................................................................................................................................95 10 List of tables Table 1. Collection of remote sensing data used in this article. ................................................... 24 Table 2. List of multispectral images used to compare 3D model contours. ................................ 28 Table 3. List of multispectral images used to compare areas between reality and 3D model. Added images are identified with *. ............................................................................................. 30 Table 4. Kappa coefficient values and overall accuracy between imagery (reality) and simulation (3D model). ................................................................................................................................... 34 Table 5. Confidence Interval analysis for the percentage of the matching area between the three- dimensional model and the sample images. .................................................................................. 35 Table 6. Landsat 8 OLI imagery used to estimate ETa through REEM and EEFlux. Source: USGS (2019). ................................................................................................................................ 58 Table 7. Comparative table of errors between the reference evaporation and the models based on remote sensors (REEM and EEFlux). Source: Rojas Villalobos. ................................................. 62 Table 8. Summary of the ranked results of the comparative statistical indicators applied to the REEM and EEFlux versus S-Penman. Source: Rojas Villalobos. ................................................ 65 Table 9. Result of the calculations of the implementation of the algorithm in Python language. Study site: the Bustillos Lagoon, Chihuahua, Mexico. Error threshold = 0.01%. * Input data. Source: Rojas Villalobos............................................................................................................... 87 Table 10. Iterative model processing times with various storage volume input values using two DTMs with different pixel dimensions. Pixel spatial resolution: 5 meters. Source: Rojas Villalobos. ..................................................................................................................................... 88 11 List of figures Figure 1. Location of Cuauhtemoc Basin and Bustillos Lagoon (Source: Rojas Villalobos with data retrieved from Google Maps and National Institute of Statistics and Informatics – INEGI, 2016). ............................................................................................................................................ 16 Figure 2. The study area of Laguna de Bustillos, Chihuahua. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). .................................................................. 21 Figure 3. Schematic of the workflow to generate the 3D model. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). .......................................................... 22 Figure 4. Components to calculate the height of the lake bottom above sea level. Source: Rojas Villalobos. ..................................................................................................................................... 23 Figure 5. Demonstration of matching areas between water surface extracted from a multispectral satellite image and the 3D model at the same reference level. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). .................................................................. 29 Figure 6. Map showing bathymetry, GPS points, derived curves from multispectral RS, and regional contours (INEGI). Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). ................................................................................................................ 32 Figure 7. Triangulated Irregular Network is representing the topobathymetric 3D model of Laguna de Bustillos. Source: Rojas Villalobos............................................................................. 33 Figure 8. 3D perspective of Laguna de Bustillos (5 times height exaggeration for better visualization). Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). .......................................................................................................................................... 33 Figure 9. Graph showing the behavior of the intersection percentage between the surfaces of the 3D model and the areas of RS images along elevation. ................................................................ 36 Figure 10. Graphs of the surface and volume equations adjusted to the 3D model. .................... 38 Figure 11. RS time series contours. The dark contour delimits the outer areas with greater 3D model performance and the internal area with less accuracy. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). .................................................................. 39 Figure 12. Schematic flow chart of the process of comparing the REEM and EEFlux models to obtain E estimations of water bodies by comparing the S-Penman equation. Source: Rojas Villalobos 2019. ............................................................................................................................ 56 Figure 13. Location of the Bustillos Lagoon and the agro-meteorological station. Source: Rojas Villalobos with data retrieved from INEGI (2019)....................................................................... 57 Figure 14. Evaporation values of S-Penman, REEM, and EEFlux during the 2017 agricultural season for the Bustillos Lagoon. Source: Rojas Villalobos with data retrieved from UNIFRUT (2019), USGS LandsatLook Viewer (2019), and EEFlux (2019). ............................................... 61 Figure 15. Comparative graphic of residuals predicted E on RS map models versus observed E (S-Penman). Source: Rojas Villalobos. ........................................................................................ 64 Figure 16. Seasonal evaporation comparison of RS models versus S-Penman data from April 4, 2017 to September 14, 2017. Source: Rojas Villalobos. .............................................................. 67 12 Figure 17. ET (crop fields) and evaporation (lagoon) comparison maps of REEM and EEFlux models in the Cuauhtemoc Valley for June 17, 2017. Source: Rojas Villalobos with data retrieved from USGS (2019) and EEFlux (2019). ........................................................................ 68 Figure 18. Study area where the algorithm was applied. The Bustillos Lagoon in Chihuahua. Source: Rojas Villalobos with data retrieved from INEGI (2019). .............................................. 82 Figure 19. DTM of the Bustillos Lagoon. Source: Rojas Villalobos with data retrieved from Rojas-Villalobos et al. (2018). ...................................................................................................... 83 Figure 20. The schematic diagram shows single-inputs and multiple-output data for iterative algorithm. Source: Rojas Villalobos. ............................................................................................ 84 Figure 21. Equations to calculate Average Depth and Maximum Depth. .................................... 84 Figure 22. Flowchart of the iterative algorithm to compute hydrologic characteristics using single-input data. Source: Rojas Villalobos. ................................................................................. 86 Figure 23. Water surface coverage map at different heights above sea level of the Bustillos Lagoon. Sources: Rojas Villalobos with data from Rojas-Villalobos (2018). ............................. 89 Figure 24. Comparative graph of volume, surface area, average depth, and maximum depth according to the height above sea level. Source: Rojas Villalobos. ............................................. 90 13 Chapter 1 Introduction The base of the economy of the Cuauhtemoc region in the Mexican state of Chihuahua is based on the apples, forage, and dairy products. In the early 20th century, one of the concerns of the Mexican government was the sparse population in the northwest, which was isolated after the Mexican Revolution in 1910, so the government implemented an immigration program to attract foreigners interested in agriculture. In March 1922, the Canadian Mennonite exodus began to populate the Cuauhtemoc region under the protection of Mexican government that offered: no conscription, no oath to the country, and no restrictions to exercise their religious principles. They would be allowed to create their schools with their teachers, and they would have an independent economic regime. When this community arrived at Cuauhtemoc, they had to change the agriculture techniques to be able to plant; therefore, they studied the soils to choose which kind of crop they would apply. At that time, alfalfa, apple orchards, barley, beans, corn, cotton, oats, wheat, and other fruit trees were planted on the soil with more humidity; the land with saline soils was utilized for grazing. Currently, most of the land with the best soil to plant, is owned by private, and ejidos own some of the reminder. The big private owners are Mennonites, and they have many financial resources to buy irrigation technology. On the other hand, many farmers are Mennonites and mestizos that have not yet modernized their cultivation techniques (75% of irrigation areas). These farmers still use flood irrigation since they arrived in the region; also, this type of irrigation is used to cover areas as long as 4,500 feet, causing loss of water by the hydrologic 14 wedge effect. These traditional irrigation methods, combined with the type of crops and intensive agricultural productivity employed in this region, have had a high impact on the static water table of the aquifer. By 2000, the National Water Commission (CONAGUA) had forecasted that the water table would decrease by 15 - 20 meters (50 - 65 feet) by 2030 if the current conditions of recharge and extraction remain (Ibañez Hernandez, 2010). However, another study made by CONAGUA in 2004 indicated that the aquifer depletion was 2.4 meters per year. The first effect on the city was scheduled water shortages in the entire city, especially in lands at higher elevation. In order to cover the water demand, two public wells were extended 60 meters depth (197 feet) in 2014. Studies in recent years have demonstrated the rapid depletion of the Cuauhtemoc aquifer due to various factors that come together in this water resource: the massive amount of extracted water, the low aquifer recharge rates, several droughts, and the irrigation techniques. The extractions of water in the Cuauhtemoc basin are heterogeneous according to land use. Excessive use of water is associated with large agricultural areas distributed throughout the basin, and in some other areas, the aquifer level is likely stable. These variations are due to the different velocities of groundwater flows and physical soil conditions (Díaz Caravantes, Bravo Peña, Alatorre Cejudo, & Sánchez Flores, 2014). According to the groundwater balance of Cuauhtemoc basin published in the Official Journal of the Federation (DOF, 2015), the groundwater inflow is 51 Mm3, vertical recharge precipitation is 41.5 Mm3, water return from irrigation is 22.7 Mm3, and extraction is 311.2 Mm3. The official document indicates a yearly deficit of 196 Mm3 of water across the basin. 15 Moreover, technical data of the study area is almost non-existent, the government are opaque in how they generate the water information, and the information that exists is inconsistent. For instance, Ortiz and Amado (2001) cite a document from the National Water Commission in 1989 where the Bustillos Lagoon has an area of 200 km2, but Landsat 8 image (Nov 03, 2015) showed that the area of the waterbody is 117 km2. Amado et al. (2016) cited a portal web of CONAGUA where the basin of Bustillos Lagoon is 4,072 km2, but actually, the area is 3,259 km2. In some other official documents, government official only describe the physiography of the lagoon location (INEGI, 2003). Given this challenging panorama, the analysis of water balance in a region involves a complex interaction of natural and anthropogenic processes that affect the quality of groundwater and long-term availability. That is why in the next three chapters of this research, three factors that assist in planning and forecasting the future state of the regional aquifer are addressed: the storage capacity of the Bustillos Lagoon, the evaporation occurring in the Bustillos Lagoon, and techniques to estimate quickly and efficiently the geomorphological variables (volume, area, and depth) of this water body. Research location The enclosed basin of Bustillos Lagoon is in the municipality of Cuauhtemoc and situated in the central-west region of the state of Chihuahua in the transition zone between the plateau and the mountains, with an area extent of 3,259 km2, as depicted in Figure 1. It is at 28º 24' 18'' north, 106º 52' 00'' west, and an altitude of 2,063 meters above sea level. Cuauhtemoc municipality is bounded by the municipalities of Namiquipa to the north, Riva Palacio to the 16 east, Cusihuiriachi and Gran Morelos to the south, and Bachiniva and Guerrero to the west (Figure 1)(INEGI 2014). The climate is warm semi-dry since it is in a transition zone between the semi-humid climate of the mountains and the desert of Chihuahua (García 1964). The geology is composed of extrusive igneous rock: rhyolite-tuff acid (29.3%), basalt (16.6%), andesite (0.1%), and volcanoclastic. The plains are composed of conglomerate (40.8%) and sandstone-cluster (0.3%). The average annual temperature is between 12° C and 20° C. The average annual rainfall varies between 300 and 500 mm per year (INEGI 2010). Figure 1. Location of Cuauhtemoc Basin and Bustillos Lagoon (Source: Rojas Villalobos with data retrieved from Google Maps and National Institute of Statistics and Informatics – INEGI, 2016). 17 Chapter 2 Topobathymetric 3D model reconstruction of shallow water bodies through remote sensing, GPS, and bathymetry Article submitted: January 6, 2018. Accepted: November 11, 2018. Journal: Tecnociencia Chihuahua. http://tecnociencia.uach.mx/numeros/v12n1/data/Topobathymetric_3D_model_reconstruction_L aguna_de_Bustillos.pdf Abstract Since there are no mathematical models that can calculate the Laguna de Bustillos’ water storage levels, water balance requires this data to understand the connectivity between this water body and the Cuauhtemoc aquifer. This article presents a new three-dimensional reconstruction technique based on a time series of multispectral remote sensing images, bathymetry, a topographic survey with high precision GPS, and regional contours. With the images of Landsat ETM+/OLI and Sentinel 2A from 2012 to 2013, 2016, and 2017, the contours of the water surface were extracted using the MNDWI and were associated with an elevation received from GPS. An Autonomous Surface Vehicle was also used to obtain the bathymetry of the lake. A topographic survey was carried out using GPS in populated areas, and the contour lines extracted from the INEGI Continuous Elevations Model 3.0. A DEM was constructed using ArcGIS 10.5.1, and surfaces and volumes were calculated at different elevations and compared with 16 Landsat TM/ETM+/OLI multispectral images from 1999 to 2018. The results showed that the mean of the average intersection area between the test images and the area extracted from the 3D http://tecnociencia.uach.mx/numeros/v12n1/data/Topobathymetric_3D_model_reconstruction_Laguna_de_Bustillos.pdf http://tecnociencia.uach.mx/numeros/v12n1/data/Topobathymetric_3D_model_reconstruction_Laguna_de_Bustillos.pdf 18 model is above 90.9% according to the confidence interval, kappa overall accuracy 95.2–99.7%, and a coefficient 89.9–99.3%. This model proved to be very accurate on a regional scale when the water level exceeded 1971.32 meters above mean sea level and useful to evaluate and administer water resources. Introduction The Laguna de Bustillos is in a region that has a high demand for groundwater for the agricultural industry, making the Cuauhtemoc aquifer the largest over-exploited aquifer in northwest Mexico (Comisión Nacional del Agua, 2016). It is necessary to provide updated data to the water balance of the basin to improve water management in the region. Because there are no known mathematical models that calculate water storage, it is imperative to develop a new technique or method that allows us to estimate the water volume contained in water bodies. The calculation of water storage of shallow water bodies requires the construction of 3D models of the terrain including the surrounding areas. Integrating techniques based on sound, spectral analysis of satellite imagery, and GPS allow researchers to increase the accuracy of the existing 3D models and expand them from the reservoir representation to a topobathymetric integrated model. Topobathymetry is a geospatial concept that integrates bathymetric and topographic data from different spatial scales, time, and sensors. The terrain model is applied to monitor coastal erosion, sea level rise, flood impact reduction programs, and coral barrier studies (Gesch et al., 2016). Digital terrain models, topography, bathymetry, and the use of water body contours are essential sources for integrating this model. Some research tried to get 3D models, but only one 19 or two data sources were used in comparison with those applied in this research. The delimitation of water bodies is an indirect way of getting contour lines through differentiating the spectral response between the green band (G) and the bands near infra-red (NIR) or the infra-red short- wave band (SWIR). The Normalized Difference Water Index (NDWI) (McFeeters, 1996) and the Modification of Normalized Difference Water Index (MNDWI) (Xu, 2006) have been used to monitor (Lu et al., 2013) changes in the extent of the lakes (Ma et al., 2007), and the location of water bodies (Rana and Neeru, 2017). Sonar is a technique that uses sound waves to calculate water depth (Knott and Hersey, 1957) and has advantages such as high accuracy (± 0.1 m), low cost, and the device can be mounted on any boat. Several types of research have used sound for mapping water bodies (McPherson et al., 2011; Popielarczyk and Templin, 2014; Giordano et al., 2015). Leon and Cohen (2012) modeled the volume of Lake Eyre in Australia using bathymetry and remote sensing. The authors used surveys realized in 1974 and 1976 with the precision of ± 0.3 m in the vertical component and up to ± 500 m in the horizontal component, which proved to be a very limited and inaccurate method. Water storage has two components: groundwater and surface water (lakes, ponds, or reservoirs) (Brooks et al., 2012). Some variations in water storage in the reservoirs are due to an underground hydraulic connection between aquifers and water bodies (Isiorho et al., 1996; Winter, 1999). These variations in the volume of water can be so drastic that large reservoirs dry up in a short time like Laguna de Bustillos had in years 2002 to 2006 and 2013 (NASA, 2017). Although there is a geohydrological study that supports recharge deficit in the aquifer, there is no information about the storage capacity of Laguna de Bustillos. The lack of information 20 encourages the main objective of this research to generate a new technique to generate a 3D topobathymetric model that contributes to the generation of updated data, which allows the deduction of variables, such as underground infiltration from the catchment area of Laguna de Bustillos. Despite these models of volumetric estimation of water bodies, the combination of more than two different topobathymetric measurement techniques had not been explored. This document proposes a unpublish new method integrating three methodologies to generate a more robust and accurate three-dimensional model. Materials and methods This study was conducted between 2016 and the first semester of 2018 in the Spatial Applications and Research Center at the New Mexico State University. The study area of Laguna de Bustillos is in the quadrant between the coordinates 28°38'51''N – 28°28'27''N and 106°57'3''W – 106°38'50''W in the municipality of Cuauhtemoc, in the state of Chihuahua (Figure 2). 21 Figure 2. The study area of Laguna de Bustillos, Chihuahua. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). This region’s climate is warm and semi-arid since it is in a transition zone between the semi-humid climate of the mountains and the Chihuahua desert (Kottek et al., 2006). The average annual temperature ranges from 6.9 to 21°C, with an average annual rainfall of about 528 mm per year (Servicio Meteorológico Nacional, 2017). The authors designed a new four-stage method to develop a 3-D topobathymetric model for the purpose of determining water storage: i) extract contour lines through a time series of remote sensing; ii) determine bathymetry; iii) perform a topographic survey (GPS-RTK); and iv) 22 extract contours from the regional terrain digital model. Also, it was included a regression analysis in determining the two equations that provide the volume and surface area using water height. The flowchart below (Figure 3) shows the modeling process. Figure 3. Schematic of the workflow to generate the 3D model. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). Bathymetry The New Mexico Water Resources Research Institute (WRRI) funded a project to build an Autonomous Surface Vehicle (ASV) to generate bathymetric data for shallow water bodies. A PVC center frame was attached to a two-hulled catamaran boat, propelled by two motors, and equipped with a GPS on the top to receive signals via satellite to provide the direction and location. An Ardupilot® system automated the catamaran navigation through an Arduino® MEGA 2560 board to receive the GPS signal while the sonar data bus decoded and recorded the 23 information on an SD card. Subsequently, the recorded points were downloaded to a computer for processing. The transducer was a Garmin® Intelliducer Thru-Hull NMEA-0183, which does not require the previous calibration and can measure from 60 cm to 200 m with a 0.1 m accuracy (Rojas-Villalobos, 2016). To construct a 3D model of the region including the bottom of the lake, the bathymetry data (depth) was transformed into topographic data (height). Figure 4 shows the schematic of the surveying process to transform to the correct topographic points. Figure 4. Components to calculate the height of the lake bottom above sea level. Source: Rojas Villalobos. The following equation ( 1 ) was developed to calculate the altitude above sea level for each bathymetric point: 24 ( 1 ) Here, ABP is the height of the bathymetric point, ASNM is the altitude above sea level of the reference level, PS is the depth of the sonar, and PR is the recorded depth. The bathymetry consisted of 5 trajectories, and the data were adjusted through the above equation using the reference levels of the survey days. A GPS-RTK was used to establish the fixed reference point corresponding to the height of the lake contour and was linked to the bathymetry obtained that day. Contour extraction from remote sensors Since the spatial resolution of remote sensing is the most important factor for delineating the contours of water bodies, Landsat ETM +, Landsat OLI (Operational Land and Imager), and Sentinel 2A (Table 1) were chosen to build the MNDWI. Table 1. Collection of remote sensing data used in this article. Sensor Acquisition date Bands (µm) Spatial resolution (m) Landsat ETM+ (USGS, 2017a) 19 May 2012; 4 June 2012; 20 June 2012; 14 January 2013; 6 February 2013; 22 February 2013 2 (Green 0.52-0.60) 5 (SWIR-1 1.55-1.75) 8 (Panchromatic) 30 30 15 Landsat OLI (USGS, 2017a) 2 August 2013; 14 June 2016, 29 August 2017; 5 September 2017 3 (Green 0.533-0.590) 6 (SWIR1 1.566-1.651) 8 (Panchromatic 0.503- 0.676) 30 30 15 Sentinel 2A (ESA, 2017) 20 March 2017; 8 June 2017 6 August 2017 3 (Green 0.542-0.577) 11 (SWIR1 1.568-1.658) 10 20 25 These images are available for free on the LandsatLook Viewer websites of the United States Geological Survey (USGS, 2017a) and the Copernicus Open Access Hub of the European Space Agency (ESA, 2017). Seven images were selected with the lowest possible cloudiness over the study area during the time the lake had gradually dried (March 2012 – August 2013). Also, six recent images were downloaded to establish the maximum lagoon extent and baseline curves for the bathymetry data (June 2016 – September 2017). Using the Semiautomatic Classification extension (Congedo, 2013) in QGIS®, atmospheric correction was applied to the images using the method of Subtraction of Dark Objects 1 (Chavez, 1996). Then, a fusion of images was performed with the panchromatic band (ETM + and OLI) using the Brovey transformation (Johnson et al., 2012) to increase the spatial resolution to 15 m before the MNDWI construction. The Normalized Difference Water Index (NDWI) was created to identify Landsat water bodies. The high relative reflectance of green (G) in the electromagnetic spectrum contrasts with the high absorption of the NIR in clear water (McFeeters, 1996). Excessive suspended matter in the water increases reflectance measurements in the NIR band (Ruddick et al., 2006), thus dramatically reducing the difference between the G-NIR bands, which makes it difficult to distinguish between water and non-water surfaces. Therefore, the NDWI method is not fit for Laguna de Bustillos due to the turbidity of water. The MNDWI suppresses this problem by replacing the NIR band with an infrared shortwave band (SWIR) because the water absorbs energy and the reflectance is low. The equation that determines the MNDWI (Xu, 2006) is: 26 ( 2 ) Where G is the green band of the electromagnetic spectrum and SWIR is the short-wave band of the infrared spectrum. The possible MNDWI values are from -1 to 1. In ArcGIS®, the raster calculator was used to apply the MNDWI equation to Landsat and Sentinel images. According to the MNDWI method, positive values represent water and negative values the surface without water. Therefore, the resulting raster was reclassified by assigning 1 to those values greater than 0 and 0 to values less than or equal to 0. From the reclassified images, the contours were extracted and examined through visual interpretation. This procedure ensures that the extracted contours correspond to the edge of the lake using false infrared color composite images and avoids errors due to the influence of the vegetation. A failure of the SLC (Scan Line Corrector) introduced strips with missing data in the Landsat ETM+ images captured on May 31st, 2003 (USGS, 2017b). Due to this error in the sensor, only segments were vectorized corresponding to the edge of the water surface. An orthometric height was assigned to the contours using the closest ABP to the contour line (<0.5 meters). When there were no bathymetry points near the line, points were selected in a buffer of 1 to 2 m on each side of each contour. The contour took the mean height following the Classic Central Limit Theorem (Erdös and Rényi, 1959; Dowdy et al., 2011). According to this theorem, when the sample size increases, the average sample will approximate a normal distribution. This procedure reduces the uncertainty and variability of bathymetric data due to 27 boat sway and sonar accuracy (Krause and Menard, 1965; Eltert and Molyneux, 1972; Schmitt et al., 2008). Topography The GPS points were measured using two SOKKIA GRX2 GNSS devices with a horizontal accuracy of 5 mm and 10 mm on the vertical axis. A GPS was established as base at the coordinate 28°27'25.1532"N and 106°47'24.9432"O at the height of 2069.08 on the WGS ellipsoid of 1984. 1006 topographic points were collected and transformed to the Mexican Gravimetric Geoid 2010 (GGM10) to generate altitude above the mean sea level (INEGI, 2015). Digital elevation model A contour was extracted at every meter from the Mexican Elevation Continuation 3.0 (CEM 3.0) of the National Institute of Statistics and Geography (INEGI, 2016). On September 5th, 2017, the water level of the lake was 1975.56 m above sea level (masl). For this reason, contour lines below 1976 m were eliminated from the regional DEM. Topobathymetric 3D model and volume estimation An MDE with a spatial resolution of 2 m was created using the four sources of elevation data using the Topo-to-Raster tool contained in the 3D analysis module of ArcGIS. This tool allows the creation of hydrologically correct lifting meshes based on the ANUDEM program (Hutchinson et al., 2011). Since the triangulated irregular network (TIN) generates more accurate volumetric calculations (Mi et al., 2007; Hanjianga et al., 2008), the DEM was converted into a TIN. The volume and water surface were calculated from 1970.50 m to 1978.9 masl every 1 mm using the ArcGIS Polygon Volume tool. 28 Statistical Evaluation Since there is no previous model to evaluate the lake storage, 16 areas of water coverage of different scenes were extracted through remote sensing (RS) when the lake was drying (real area) ( Table 2). Table 2. List of multispectral images used to compare 3D model contours. Sensor Acquisition date Water surface (km2) Landsat TM (USGS, 2017a) June-25-1999 May-26-2000 June-11-2000 March-17-2001 April-02-2001 Abril-27-2001 November-24-2002 September-15-2003 December-21-2006 99,585,900 92,322,000 90,583,200 77,341,500 70,556,400 64,676,700 41,630,400 64,507,500 109,260,000 Landsat ETM+ (USGS, 2017a) January-27-2000 May-05-2001 August-28-2002 104,792,438 61,820,100 68,073,300 Landsat OLI (USGS, 2017a) May-01-2014 June-02-2014 August-28-2014 October-08-2014 87,509,700 79,517,700 109,547,000 118,406,700 29 The area of each scene was used to extract the corresponding contour line from the 3D model and generate the area. Using ArcGIS, the intersection of the two layers was the area of a coincidence that was statistically evaluated (Figure 5). Figure 5. Demonstration of matching areas between water surface extracted from a multispectral satellite image and the 3D model at the same reference level. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). Some Landsat ETM + and OLI images were replaced with recent Sentinel 2 images (early 2018) to distribute the extracted contours along the height through the 3D model ( 30 Table 3). This procedure is used to evaluate the model accuracy (reality vs. model). Table 3. List of multispectral images used to compare areas between reality and 3D model. Added images are identified with *. Sensor Acquisition date Water surface (km2) Landsat TM (USGS, 2017a) June-25-1999 May-26-2000 June-11-2000 March-17-2001 April-02-2001 99,585,900 92,322,000 90,583,200 77,341,500 70,556,400 Landsat ETM+ (USGS, 2017a) January-27-2000 May-05-2001 August-28-2002 104,792,438 61,820,100 68,073,300 Sentinel 2 (ESA, 2017) May-04-2016* July-23-2016* January-14-2018* April-04-2018* 109,321,000 105,033,000 133,912,000 131,504,000 Because of the surface area changes according to the elevation of the water surface, it is not possible to evaluate the efficiency of the model directly. For this reason, the relationship between the coincidence surface and the reference area of the satellite image were used. The maximum possible relation between both areas is 100% because the level curves obtained from the 3D model are directly related to the waterbody contours. The water/non-water coverage maps of the model and the satellite images of each year ( Table 3) were analyzed using the Kappa statistic (K-hat) through QGIS (QGIS, 2018) and Semi- Automated Classification Plugin (Congedo, 2013). The Kappa coefficient and overall accuracy 31 allows us to know the degree of agreement between the 3D model and the water body surface (Card, 1982; Jensen, 2007; Congalton and Green, 2008; Lillesand et al., 2014) . Also, the t- statistical distribution was applied to find the lower limit of the 95% Confidence Interval and estimated the range of acceptable match surface values (from Table 2) according to the sample mean (Dowdy et al., 2011) ( 3 ). ( 3 ) Where X is the mean of the sample, α is the level of significance, ν is the degrees of freedom (n -1), s is the standard deviation, and n is the sample size. Finally, two equations were generated representing the area of the water surface and the volume contained in the lake according to the elevation of the water surface. Results and discussion Figure 6 shows the sources of data used for the reconstruction of the topobathymetric model: 13 contours from remote sensors, 29,715 bathymetry points, 1,006 GPS points, and INEGI contours. 32 Figure 6. Map showing bathymetry, GPS points, derived curves from multispectral RS, and regional contours (INEGI). Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). As a result of the reconstruction data process, Figures 7 and 8 show the 3D topobathymetric model and a 3D perspective of the Laguna de Bustillos. The results show that the deepest point of the lake is at 1970.215 masl, the maximum depth is 3.785 m when the water level reaches the 1974 masl, the water storage is 324.4 Mm3, and the average depth is 1.37 m. 33 Figure 7. Triangulated Irregular Network is representing the topobathymetric 3D model of Laguna de Bustillos. Source: Rojas Villalobos. Figure 8. 3D perspective of Laguna de Bustillos (5 times height exaggeration for better visualization). Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). 34 Since the Kappa statistic shows the difference between classified values of the satellite image (reference data) and the surface of water body generated by the 3D model, the coincidence is expected to be high. Typically, Kappa values greater than 0.80 represent a strong match between the compared data. The result of the comparison shows an overall accuracy higher than 95.21% and the K-hat coefficients above 0.899. Table 4 shows the increase of the values of overall accuracy and the Kappa coefficient when the water level is higher. Table 4. Kappa coefficient values and overall accuracy between imagery (reality) and simulation (3D model). Date Sensor Surface (km2) Elevation (m) Depth Average (m) K-hat Overall Accuracy (%) 25/06/1999 Landsat TM 99.59 1971.713 0.710 0.9627 98.150 27/01/2000 Landsat ETM+ 104.86 1972.034 0.987 0.9669 98.383 26/05/2000 Landsat TM 92.32 1971.460 0.501 0.9347 96.738 11/06/2000 Landsat TM 90.58 1971.442 0.493 0.9316 96.582 17/03/2001 Landsat TM 77.35 1971.284 0.409 0.9079 95.499 02/04/2001 Landsat TM 70.56 1971.235 0.397 0.8993 95.212 05/05/2001 Landsat ETM+ 61.82 1971.168 0.383 0.8993 95.480 28/08/2002 Landsat ETM+ 68.07 1971.228 0.405 0.9289 96.669 04/05/2016 Sentinel 2 109.32 1972.809 1.709 0.9939 99.710 23/07/2016 Sentinel 2 105.03 1972.050 1.000 0.9883 99.430 14/01/2018 Sentinel 2 133.91 1975.857 4.195 0.9787 99.171 04/04/2018 Sentinel 2 131.55 1975.554 3.955 0.9687 98.750 35 It is observed that the values of elevation that are between 1971.168 and 1971.284 have a value of K-hat less than 0.9289 and are associated with water coverage less than 80 km2. When the water level rises above 1971.284 m, the Kappa indicator increases its value above 0.93, reaching levels of 0.99. Also, low K-hat values (0.8993 – 0.9289) are associated with low depth averages (<0.41 m) in contrast to those K-hat values above 0.96 that are in depth averages greater than 0.71 m. Conversely, with a confidence level of 95%, the mean of the percentage of matching areas between the satellite images and the 3D model is greater than 90.9% (Table 5). Table 5. Confidence Interval analysis for the percentage of the matching area between the three- dimensional model and the sample images. Mean 0.934663471 Degree of freedom 15 Standard Error 0.0145891 α 0.05 Median 0.954318 t0.05,15 1.753 Standard deviation 0.0583564 t0.05,15 Std. Error 0.025574712 Simple variance 0.0034055 IC0.95: δ - t0.05,15 Std. Error 0.9090888 Sample size 16 Below the contour 1971.325 m, four of the six comparisons are below the lower limit of the confidence interval (Figure 9). 36 Figure 9. Graph showing the behavior of the intersection percentage between the surfaces of the 3D model and the areas of RS images along elevation. The mean area intersected below the reference level is only 88.91%, while in the upper range, it is 97.01%. Two equations were generated that estimated the area of water coverage according to the depth of the lake. The first equation calculated the volume below the 1971.325 masl and the second equation calculated the remaining volume above it. Similarly, two other equations were generated estimating the amount of water in the lake. The determination coefficients (R2) for the estimated equations are greater than 0.9882; this indicates that the equations obtained are suitable for the topobathymetric model within the extent limits of the lake (Figure 10). 37 38 Figure 10. Graphs of the surface and volume equations adjusted to the 3D model. Conclusions Four different techniques, such as bathymetry, GPS-RTK points, and contour lines extracted from the remote sensors, were decisive in creating this new three-dimensional modeling methodology for water bodies. Its efficiency is demonstrated after the statistical analysis applied. According to the results obtained in the Kappa analysis and the confidence interval, the 3D model is a robust and precise model (Kappa>0.80). Three processes were important in the construction of the model: • The use of high precision GPS helped in fixing the reference height points of the contours of the most recent satellite images (2015 – 2018) with great precision and accuracy. • The bathymetric points linked to the current height of the water level of the lake were instrumental in establishing the height above sea level at the bottom of the lagoon. • The related height between the bathymetric points and the levels closest to the bottom of the lake was extracted from the satellite images (1999 – 2002). 39 Additionally, it was observed that the segments of the contours extracted from the Landsat ETM+ images with an error in the SLC (USGS, 2017b) influenced the relative low efficiency (0.8993 < Kappa <0.9079) of the model below 1971.325 masl. On the other hand, effectiveness in the top height ranges from 1971.5 to 1974 masl was as a result of the spatial resolution of the satellite images of Landsat OLI (15 m panchromatic) and Sentinel 2 (10 m) (Figure 11). Figure 11. RS time series contours. The dark contour delimits the outer areas with greater 3D model performance and the internal area with less accuracy. Source: Rojas Villalobos with data retrieved from LandsatLook Viewer (USGS, 2017a). 40 Although this 3D hydrological model is very robust to be used in the administration of water in the basin, special care must be taken in forecasting floods in rural-urban areas. The model simulates much of the flooded areas of Mennonite farmers, but the 3D model should not be used to prevent flood risks due to the topographic complexity with dams and ditches. In future work, researchers should continue the bathymetric survey with greater data density using a sonar with increased accuracy to further the model’s efficiency. The acquisition of more bathymetric data will allow replacing contours extracted from the oldest images such as Landsat ET and ETM +. Additionally, the photogrammetric triangulation could be of great benefit in urban and agricultural zones to delineate more accurate topography. This development is the first step to estimate the volume of water in the Laguna de Bustillos as this work produces estimates that approximate the actual values and such research is relevant to water management in the region. 41 References Brooks, K.N., P.F. Ffolliott, and J.A. Magner 2012. Hydrology and the Management of Watersheds. John Wiley & Sons. Iowa. 533 p. Card, D.H. 1982. Using known map category marginal frequencies to improve estimates of thematic map accuracy. Photogramm. Eng. Remote Sens 48: 431–439. Chavez, P.S. 1996. Image-based atmospheric corrections-revisited and improved. Photogramm. Eng. Remote Sens 62: 1025–1035. Comisión Nacional del Agua, 2016. Estadísticas del agua en México. Secretaría de Medio Ambiente y Recursos Naturales, México. Congalton, R.G., Green, K., 2008. Assessing the accuracy of remotely sensed data: principles and practices. CRC press. Congedo, L., 2013. Semi-automatic classification plugin for QGIS. Sapienza Univ. Rome. Dowdy, S., Wearden, S., Chilko, D., 2011. Statistics for Research. John Wiley & Sons. Eltert, J.F., Molyneux, J.E., 1972. The long-distance propagation of shallow water waves over an ocean of random depth. J. Fluid Mech 53: 1–15. Erdös, P., Rényi, A., 1959. On the central limit theorem for samples from a finite population. Publ Math Inst Hung. Acad Sci 4: 49–61. ESA, 2017. Copernicus Open Access Hub [WWW Document]. URL https://scihub.copernicus.eu/ (accessed 8.23.17). Gesch, D.B., Brock, J.C., Parrish, C.E., Rogers, J.N., Wright, C.W., 2016. Introduction: Special Issue on Advances in Topobathymetric Mapping, Models, and Applications. J. Coast. Res 76: 1–3. Giordano, F., Mattei, G., Parente, C., Peluso, F., Santamaria, R., 2015. MicroVeGA (micro vessel for geodetics application): A marine drone for the acquisition of bathymetric data for GIS applications. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci 40: 123. Hanjianga, X., Limina, T., Longa, S., 2008. A strategy to build a seamless multi-scale TIN-DEM database. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci 37. Hutchinson, M.F., Xu, T., Stein, J.A., 2011. Recent progress in the ANUDEM elevation gridding procedure. Geomorphometry 2011, 19–22. INEGI, 2016. Continuo de Elevaciones Mexicano 3.0 (CEM 3.0) [WWW Document]. Datos Relieve. URL http://www.inegi.org.mx/geo/contenidos/datosrelieve/continental/continuoelevaciones.as px (accessed 3.16.16). INEGI, 2015. El geoide gravimétrico mexicano 2010. Instituto Nacional de Estadística y Geografía, México. Isiorho, S.A., Matisoff, G., Wehn, K., 1996. Seepage relationships between Lake Chad and the Chad aquifers. Groundwater 34: 819–826. Jensen, J.R., 2007. Remote Sensing of the Environment: An Earth Resource Perspective 2/e, 2nd ed, Prentice Hall Series. Pearson Education. 42 Johnson, B.A., Tateishi, R., Hoan, N.T., 2012. Satellite image pansharpening using a hybrid approach for object-based image analysis. ISPRS Int. J. Geo-Inf 1: 228–241. Knott, S.T., Hersey, J.B., 1957. Interpretation of high-resolution echo-sounding techniques and their use in bathymetry, marine geophysics, and biology. Deep Sea Res 1953 (4): 36–44. Kottek, M., Grieser, J., Beck, C., Rudolf, B., Rubel, F., 2006. World map of the Köppen-Geiger climate classification updated. Meteorol. Z. 15, 259–263. Krause, D., Menard, H., 1965. Depth distribution and bathymetric classification of some sea- floor profiles. Mar. Geol. 3, 169–193. Leon, J.X., Cohen, T., 2012. An improved bathymetric model for the modern and palaeo Lake Eyre. Geomorphology 173, 69–79. Lillesand, T., Kiefer, R.W., Chipman, J., 2014. Remote sensing and image interpretation. John Wiley & Sons. Lu, S., Ouyang, N., Wu, B., Wei, Y., Tesemma, Z., 2013. Lake water volume calculation with time series remote-sensing images. Int. J. Remote Sens. 34, 7962–7973. Ma, M., Wang, X., Veroustraete, F., Dong, L., 2007. Change in area of Ebinur Lake during the 1998–2005 period. Int. J. Remote Sens. 28, 5523–5533. McFeeters, S.K., 1996. The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int. J. Remote Sens. 17, 1425–1432. McPherson, K.R., Freeman, L.A., Flint, L.E., 2011. Analysis of methods to determine storage capacity of, and sedimentation in, Loch Lomond Reservoir, Santa Cruz County, California, 2009 (USGS Numbered Series No. 2011–5141), Scientific Investigations Report. U.S. Geological Survey, Reston, VA. Mi, H., Zai, J., Jiang, X., 2007. Contrast and Analysis of Reservoir Storage Calculation Methods [J]. Surv. Mapp. Geol. Miner. Resour. 2, 000. NASA, 2017. History Landsat Science [WWW Document]. URL https://landsat.gsfc.nasa.gov/about/history/ (accessed 9.1.17). Popielarczyk, D., Templin, T., 2014. Application of integrated GNSS/hydroacoustic measurements and GIS geodatabase models for bottom analysis of Lake Hancza: the deepest inland reservoir in Poland. Pure Appl. Geophys. 171, 997–1011. QGIS Development Team (2018). QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org Rana, H., Neeru, N., 2017. Water Detection using Satellite Images Obtained through Remote Sensing. Adv. Comput. Sci. Technol. 10, 1923–1940. Rojas-Villalobos, H.L., 2016. 3D bathymetric model of a shallow lagoon measured by a solar powered low-cost autonomous surface vehicle prototype in Cuauhtémoc, Chihuahua, Mexico. (Academic Research), Student Water Research Awards 2015-2016. New Mexico Water Resources Research Institute, Las Cruces, NM. Ruddick, K.G., De Cauwer, V., Park, Y.-J., Moore, G., 2006. Seaborne measurements of near infrared water-leaving reflectance: The similarity spectrum for turbid waters. Limnol Ocean. 51, 1167–1179. 43 Schmitt, T., Mitchell, N.C., Ramsay, A.T.S., 2008. Characterizing uncertainties for quantifying bathymetry change between time-separated multibeam echo-sounder surveys. Cont. Shelf Res. 28, 1166–1176. Servicio Meteorológico Nacional, 2017. Información Climatológica (1981-2010) [WWW Document]. URL http://smn.cna.gob.mx/es/informacion-climatologica-ver- estado?estado=chih (accessed 9.4.17). USGS, 2017a. LandsatLook Viewer [WWW Document]. URL https://landsatlook.usgs.gov/ (accessed 8.22.17). USGS, 2017b. SLC-off Products: Background | Landsat Missions [WWW Document]. URL https://landsat.usgs.gov/slc-products-background (accessed 8.28.17). Winter, T.C., 1999. Relation of streams, lakes, and wetlands to groundwater flow systems. Hydrogeol. J. 7, 28–45. Xu, H., 2006. Modification of normalized difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens. 27, 3025–3033. 44 Chapter 3 Comparison of evaporation estimates from REEM and EEFlux models in a shallow water body. Case: Bustillos Lagoon, Chihuahua, Mexico. Article submitted: November 19, 2019. Status: In review. Journal: Tecnología y Ciencias del Agua (IMTA). Abstract Water body evaporation (E) within endorheic basins in semiarid areas is a critical factor for the determination of the water balance. Unfortunately, the Bustillos Lagoon has dried up completely six times during this century, and there are no records of the evaporation rate. Furthermore, accurate E measurements can provide valuable information for the sustainable management of water resources for protecting wild habitat in the face of climate change scenarios. Evaporation can be estimated, however, through methods as efficient as Penman using variables from agroclimatic stations, such as wind velocity, net radiation, relative humidity, and air temperature, which have a spatiotemporal variability. Within the evaporation models based on remote sensors (RS) is the surface energy balance model (SEB), which has been applied to different methodologies and extends the measurements of evapotranspiration (ET) at a regional level. SEB-based methodologies use physical principles with minimal weather data requirements to estimate ET. Hence, this article compares two methodologies that estimate evaporation using RS: The Regional Evapotranspiration Estimate Model (REEM) and the Earth Engine Evapotranspiration Flux (EEFlux). The comparison of ET measurements obtained from REEM and EEFlux for seven Landsat OLI scenes in the agriculture cycle of April to September applied 45 against to the simplified Penman equation showed that the REEM performed better (d=94%) than the EEFlux (d=68%) for the indicated period. Although the comparison of REEM and EEFlux shows accurate E measurements (REEM), gridded weather data (EEFlux) need to be improved by increasing the scale using local information. Introduction The Bustillos Lagoon is the largest water body (~100 km2) in the Cuauhtemoc Valley (in Chihuahua, Mexico), which is in an endorheic basin (3,302 km2).The climate is semiarid, and agriculture is intensive. High competition for water resources among stakeholders (Díaz Caravantes, Bravo Peña, Alatorre Cejudo, & Sánchez Flores, 2014) has exerted high pressure on the aquifer. According to Mexican authorities, this phenomenon has caused the aquifer to be overexploited (Diario Oficial de la Federación, 2015). For this reason, farmers have made dams and ditches to divert and retain a small part of the tributary flows before they reach the Bustillos Lagoon. These practices, however, limit the source of water that supplies it. The Bustillos Lagoon, like any water body, is essential for its thermoregulatory climate function in the region as it absorbs heat fluxes and releases moisture (Rooney & Bornemann, 2013; Subin, Murphy, Li, Bonfils, & Riley, 2012). In addition, it is ecologically important as a resting place for migratory waterbirds (Mireles & Mellink, 2017). Aquatic systems in semiarid areas are susceptible to drastic variations in water levels, which affects the aquatic life (Amado-Álvarez, Pérez Cutillas, Ramírez Valle, & Alarcón Cabañero, 2016) that feeds waterbirds. If water resources are not correctly managed, regional sustainability will be jeopardized, causing the desiccation of water bodies such as the Aral Sea between Kazakhstan and Uzbekistan (Gross, 2017), Lake Chad in 46 the borders of Niger, Nigeria, Cameroon and Chad (Okpara, Stringer, Dougill, & Bila, 2015) and Lake Urmia in Iran (AghaKouchak et al., 2015). These water bodies are drying up because of the diversion of tributary rivers to agricultural fields, droughts, and upstream competition for water. Evaporation data of the lagoon are required to establish administrative water resource policies to avoid catastrophic scenarios and to conserve the water balance in the Cuauhtemoc Basin Evapotranspiration (ET) is a process that combines the evaporation of water surfaces, the evaporation of soil moisture, and the transpiration of vegetation (Erickson et al., 2008). Evaporation is part of ET, which is governed by aerodynamic and energy equations (Penman, 1948). Under this approach, it is possible to estimate the evaporation of a water body through the calculation of ET using remote sensing techniques. The most effective (and costly) techniques for measuring evapotranspiration are lysimeters or eddy covariance flux stations (Hirschi, Michel, Lehner, & Seneviratne, 2017), which do not exist in the study area. Because of this situation, it is necessary to explore emerging alternative methodologies for measuring ET. Rohwer (1931) developed evaporation coefficients (Kpan) for the evaporation pan method (U.S. Class A pan) for each month of the year. The problem with this approach is that the method used lakes in the state of Colorado as research sites. These sites contained clear water, and the physical aspects of the metal pan container affected evaporation measurements (Fu, Charles, & Yu, 2009; Rayner, 2007). In addition, a pan coefficient is a function of depth and surface area of the lake that is being estimated. The Bustillos Lagoon has particular characteristics that make it different from other lagoons and lakes. For example, in addition to being a shallow lagoon, 47 turbidity is high, caused by the content of suspended material (Álvarez, Cutillas, Valle, & Cabañero, 2016; Amado-Alvarez et al., 2019). Radiation flux from the sun penetrates deeply into the water column in clear water conditions, absorbing energy (Smith & Tyler, 1967). Under conditions of turbidity and low depth (<3 m) (Rojas-Villalobos, Alatorre-Cejudo, Stringman, Samani, & Brown, 2018), solar radiation is scattered by suspended particles on the surface layer. Therefore, the water temperature is increased, resulting in more evaporation (Kirk, 1985). Under these conditions, it is not possible to apply pan evaporation coefficients, since the physical characteristics change in each lake or lagoon. The methods for calculating evaporation can be classified into those based on: daytime air temperature range such as that of Papadakis (Papadakis, 1965); air temperature and day length such as Hamon (Hamon, 1960), and Blaney-Criddle (Blaney & Criddle, 1957); solar radiation and air temperature such as Jensen-Haise (Jensen & Haise, 1963), Makkink (Makkink, 1957), and Stephens-Stewart (Stephens & Stewart, 1963); heat flux and water vapor flux (combination) such as Priestley-Taylor (Priestley & Taylor, 1972), De Bruin-Keijman (De Bruin & Keijman, 1979), Penman (Penman, 1948), Brutsaert-Stricker (Brutsaert & Stricker, 1979), and De Bruin (De Bruin, 1978). Although these methods can offer good evaporation approximations, estimates are local at the point of the reference weather station. Given this limitation, remote sensing (RS) techniques expand measurements to the regional scale in a cost-effective way. There are different satellite-based methods established on physical relationships and theoretical foundations. Zhang, Kimball, & Running (2016) classified ET retrieval methods in eight groups: i) Penman-Monteith (PM) (Cleugh, Leuning, Mu, & 48 Running, 2007; Li et al., 2017); ii) Priestley-Taylor (PT) (Martínez Pérez, García-Galiano, Martin-Gorriz, & Baille, 2017); iii) water-carbon linkage (WCL)(Fisher et al., 2018); iv) water balance (WB) (Reitz, Senay, & Sanford, 2017); v) maximum entropy production (MEP)(H. Wang, Tetzlaff, & Soulsby, 2017); vi) surface energy balance (SEB)(Senkondo, Munishi, Tumbo, Nobert, & Lyon, 2019); vii) Ts-VI space (TVI) (Zhu, Jia, & Lv, 2019); and viii) empirical and other methods (EO). Each physical-theoretical basis reported by these groups has advantages and restrictions. For instance, PM models have a robust physical base, but on the other hand, the forcing of meteorological variables induces and propagates uncertainty in the evaporation estimate. The simplified PM model is the theoretical basis of PT as a primary governing equation by adding semiempirical equations. The estimations of the water-carbon linkage method use the advantages of carbon processes, which increases uncertainty in carbon fluxes caused by forcing climatological data. The theory of nonequilibrium thermodynamics is the basis of the MEP model, which requires few enforced climatological variables but requires continuous surface temperature measurements. The SEB models require minimum local weather data and RS, but they are susceptible to temperature deviations and need clear-sky conditions. TVI models have low-temperature sensitivity but require clear-sky conditions, and they oversimplify TVI space relationships. A weak theoretical base of empirical models does not make them a robust option for application in water management policies. Within the SEB classification, there are two methodologies with a strong physical- theoretical bases: the regional evapotranspiration estimate model (REEM) (Hewitt, Fernald, & 49 Samani, 2018; Kıvrak, Bawazir, Samani, Steele, & Sönmez, 2019; A. Samani & Bawazir, 2015; Z. Samani, Bawazir, Bleiweiss, et al., 2007; Z. Samani, Skaggs, & Bleiweiss, 2005) and the Earth Engine Evapotranspiration Flux (EEFlux) (Allen et al., 2015; Ayyad, Al Zayed, Ha, & Ribbe, 2019), which is a version of mapping evapotranspiration at high resolution with internalized calibration (METRIC) (Allen, Tasumi, & Trezza, 2007; Allen, Tasumi, Trezza, et al., 2007). REEM and METRIC use the same physical basis of Surface Energy Balance Algorithms for Land (SEBAL) (Bastiaanssen, Menenti, Feddes, & Holtslag, 1998; Bastiaanssen, Pelgrum, et al., 1998) but with some differences in sensible heat flux (H) estimation and net radiation (Rn). EEFlux is an integration of the METRIC model in the Google Earth Engine platform, which uses Landsat satellite images, NLDAS and CFSv2 gridded weather data (the United States and rest of the world, respectively) for calibrating the METRIC model (Allen, Tasumi, & Trezza, 2007; Allen, Tasumi, Trezza, et al., 2007; Irmak et al., 2012). Also, remote sensors can estimate the evaporation of water bodies through the relationship with the reference evapotranspiration of agroclimatological stations and thus have the basis for establishing policies about consumptive water use. Because of the particular semiarid climatic conditions of the Cuauhtemoc Valley, as well as the turbidity and shallowness of the Bustillos Lagoon, the objective of this paper is to examine the effectiveness and performance of two evapotranspiration models based on remote sensors (REEM and EEFlux) against the simplified Penman (S-Penman) equation (Valiantzas, 2006)(derived from the Penman equation) to estimate the E of this water body. 50 Background Reference evapotranspiration model A well-known proven method for estimating evaporation from a free water surface is the Penman (Penman, 1948) equation, which is widely used around the world (Bozorgi, Bozorg- Haddad, Sima, & Loáiciga, 2018; Cabrera, Anache, Youlton, & Wendland, 2016; B. Wang, Ma, Ma, Su, & Dong, 2019). This research used a simplified version of the Penman (S-Penman) equation, which uses standard climatological records, such as solar radiation, air temperature, relative humidity, and wind speed at a 2-m height above the ground surface (Valiantzas, 2006), as noted below: (1) where E is the evaporation (mm d-1), α is albedo (0.08 for water), Rs is the solar radiation data estimated from measured daytime hours (MJ m-2 d-1), T is the mean air temperature (ºC), Ra is the extraterrestrial solar radiation (MJ m-2 d-1), RH is the mean relative humidity (%), au is equal to 1 when the wind function is used from the original Penman equation(1948), and u is the average wind velocity (m s-1). Brief description of remote sensing models The surface energy balance equation (Bastiaanssen, Menenti, et al., 1998; Bastiaanssen, Pelgrum, et al., 1998) is the foundation of ET models based on remote sensors such as REEM (Samani et al., 2007; Samani, Bawazir, Bleiweiss, Skaggs, & Schmugge, 2006) and EEFlux (Allen et al., 2015; Allen, Tasumi, & Trezza, 2007), which determines the latent heat flux that 51 represents the residual of the surface equation of energy used in the process of evapotranspiration. The equation can be expressed as (2) where LE is the latent heat flux of vaporization, Rn is the net radiation at the surface, G is the soil heat flux, and H is the sensible heat flux into the air. For REEM, the units of the surface energy balance equation are in MJ m-2 day-1), while EEFlux uses Wm-2. The different components of the equation can be solved separately through energy flux models. Below are the fundamental descriptions of the REEM, METRIC, and EEFlux models and their main features. REEM Samani and other researchers developed a methodology to calculate Rn (Z. Samani, Bawazir, Skaggs, et al., 2007; Z. Samani, Bawazir, Bleiweiss, et al., 2007): (3) where Rn is the daily net radiation (MJ m-2 day-1), Rni is the instantaneous clear sky net radiation (W m-2), Rs the daily shortwave solar radiation (MJ m-2 day-1), Rsi the instantaneous shortwave solar radiation (W m-2), Ta is the mean daily temperature (°K), and Ti is the instantaneous air temperature (°K). The instantaneous net radiation is the difference between incoming and outgoing fluxes and is estimated (Bastiaanssen, 1995) as 52 (4) where Rni is the instantaneous net radiation (W m-2), α the surface albedo (nondimensional), Rsi is the instantaneous incoming shortwave radiation (W m-2), RL↓ is the instantaneous incoming longwave radiation (W m-2), RL↑ is the instantaneous outgoing longwave radiation (W m-2), and ε0 is the surface emissivity (nondimensional). The detailed process to obtain Rni is outlined by Samani et al. (2007b). The instantaneous soil heat flux (Gi) was calculated at the time when the satellite overpassed the study site using a normalized difference vegetation index (NDVI) (Z. Samani, Sammis, Skaggs, Alkhatiri, & Deras, 2005) by the next equation: (5) The aerodynamic equation (Bastiaanssen, 1995) and the Monin–Obukhov similarity theory (Monin & Obukhov, 1954) were combined to estimate instantaneous sensible heat flux. The aerodynamic equation is expressed as (6) where ρa is the air density (kg m-3), Cp is the specific heat of air (1,004 J (kg -1 K-1)), T0 is the aerodynamic surface temperature (°K), Ta is the air temperature (°K), rah is the aerodynamic surface resistance, and dT is the air temperature gradient calculated through a Bastiaanssen (2005) equation. Moreover, dT needs a and b constants for a calibration, for which they were 53 empirically computed by selecting two pixels called “hot and cold pixels” taken from the image. These pixels represent extreme conditions: one of aridity (latent heat flux close to zero for dry soil) and the other of humidity (sensible heat flux close to zero for well-irrigated crop), respectively. The Hi equation was used to calculate dT. The cold pixel took the sensible heat value. Because there is no ET on dry bare soil, instantaneous latent heat was set to zero, and Rni and Gi could be calculated. The hot pixel was estimated as the Hi value by calculating the residual of the energy balance: (7) The ground surface wind speed data (2 m) was extrapolated to 200 m, and an iterative stability correction model based on the Monin–Obukhov similarity theory was used to estimate the aerodynamic resistance (rah) (Allen, Tasumi, & Trezza, 2007; Bastiaanssen, 1995) for each pixel. The Hi and dT were calculated for each pixel after calibration constants were estimated. The Gi and Rni were calculated at the time of the satellite overpass for the study area. The detailed process for obtaining the ET is outlined by Samani et al. (2006, 2007a, 2007b). METRIC The net radiation (Rn) is the balance of all outgoing radiant fluxes and all incoming radiant fluxes, including solar radiation and radiation in the thermal band. METRIC uses the same Rn equation as the REEM: (8) 54 where the net radiation is in W m-2, RS is the incoming solar radiation (W m-2), α is the albedo of surface (nondimensional), RL↓ is the incoming longwave radiation (W m-2), RL↑ is the outgoing longwave radiation (W m-2), and ε0 is the thermal emissivity of the surface (nondimensional). METRIC uses the same algorithms to compute Rn as the REEM. The process to determine Rn is detailed by Allen et al. (2007a, 2007b). In METRIC, G is estimated by the following equations defined by Tasumi et al. (2003), which depend on the net radiation and the leaf area index (LAI) vegetation: (9) (10) where Ts is the temperature on the near surface (°K). In addition, “cold” and “hot” pixels are used in METRIC, which employs the same algorithm to calculate H in Eq. 5 but with differences in pixel selection. Because surface wetness has higher values than other surrounding vegetation crops, the cold pixel assumes 1.05 times ETref, which is calculated from the standardized ASCE Penman-Monteith equation (ASCE– EWRI, 2005). As in REEM, the hot pixel is anchored to a dry agricultural surface free of vegetation, which assumes that latent heat flux is equal to 0. A detailed process can be found in Allen et al. (2007a, 2007b). 55 EEFlux The algorithms used in METRIC were adapted to the EEFlux using JavaScript and Python APIs to compute the ET automatically. While METRIC uses a weather station to calibrate the model at the runtime, the EEFlux uses gridded weather data sets from external sources to estimate at-surface reflectance, autocalibration, and the daily soil-water evaporation process. These sources are the NLDAS (with a 12-km grid size), the GridMET, and Daymet data sets for the United States. Furthermore, CFSv2 (with a 10-km grid size) provides gridded weather data for the rest of the world. Irmak et al. (2012) and Allen et al. (2015) outlined the implementation of the METRIC equations in EEFlux. Material and methods The processes that integrate the methodology for comparing the performance of REEM and EEFlux with the S-Penman are noted below in Figure 12: 56 Figure 12. Schematic flow chart of the process of comparing the REEM and EEFlux models to obtain E estimations of water bodies by comparing the S-Penman equation. Source: Rojas Villalobos. Study area This study was conducted during the agricultural cycle from April 2017 to September 2017 in the Cuauhtemoc Valley. The Bustillos Lagoon is a shallow endorheic freshwater body in the municipality of Cuauhtemoc, in the Mexican state of Chihuahua. The lagoon is in latitude 28°33’59.36” N and longitude 106°46’7.33” W. The lagoon has an approximately oval shape, of which the major axis is 17 km, and the minor is 8 kilometers with an average depth of 1.7 m. In addition, the area can fluctuate to around 100 km2 (Figure 13)(Rojas-Villalobos et al., 2018). Currently (August 2019), the surface of Bustillos Lake is 116.7 km2; moreover, it stores 312.7 57 hm3 and has an average depth of 2.68 m. The turbidity of the lagoon water is closely related to the shallow depth and high concentrations of sediment carried by the tributaries. Additionally, surface water erosion in the region is mainly due to extensive agriculture, sparse riparian vegetation, and the deforestation of the slopes of the mountain ranges that delimit the basin (Álvarez et al., 2016; Amado-Alvarez et al., 2019). Figure 13. Location of the Bustillos Lagoon and the agro-meteorological station. Source: Rojas Villalobos with data retrieved from INEGI (2019). 58 Agro-meteorological data An agroclimatic station, ADCON™, located 4.5 kilometers west of the Bustillos Lagoon at 28°34'11.5"N, 106°54'29.4"W and 2004 m.a.s.l provided hourly meteorological data that REEM required to calculate the ET for each date from downloaded Landsat 8 OLI satellite images. In addition, the agroclimatic station provided data for computing E by using the standardized S-Penman equation (Valiantzas, 2006) through TR1 Combi sensors for temperature and relative humidity, as well as pyranometers (SP Lite and CMP3), and wind speed. Landsat 8 OLI selection Seven Landsat 8 OLI images (Table 6), from two different Paths were chosen for continuity in the temporal and geographical space between the beginning (April) and the end (September) of the agricultural cycle in the Cuauhtemoc Basin. Additionally, the images met no cloud criteria (clear-sky) in the study area. For this reason, the intersection strip between Path 32 and Path 33 was used to estimate the ET. Table 6. Landsat 8 OLI imagery used to estimate ETa through REEM and EEFlux. Source: USGS (2019). Date DOY Overpass time (local time) Scene 07-04-2017 97 10:33:49 LC80320402017097LGN00 23-04-2017 113 10:33:40 LC80320402017113LGN00 30-04-2017 120 10:39:46 LC80330402017120LGN00 09-05-2017 129 10:33:40 LC80320402017129LGN00 16-05-2017 136 10:39:56 LC80330402017136LGN00 17-06-2017 168 10:40:12 LC80330402017168LGN00 14-09-2017 257 10:34:24 LC80320402017257LGN00 59 REEM and EEFlux raster The satellite images were radiometrically calibrated and atmospherically corrected using the ENVI® software through the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH™) tool. Once the satellite images were processed for obtaining the ETa through the REEM, the ETa layers of the EEFlux model were downloaded from the web portal (http://eeflux- level1.appspot.com/). Lagoon delineation The sampling was carried out through a lagoon polygon that was created using the Modified Normalized Difference Water index (MNDWI), which discretizes the water surface from the rest of the image (Xu, 2006). The outline of the polygon of the lagoon was adjusted by 50 meters to reduce water detection errors on the shore caused by expanding and contracting throughout the agricultural cycle. Statistical evaluation Statistical comparison was performed using the relationship between the observed values (Oi)(S-PENMAN) and the estimated or predicted values (Pi) (REEM and EEFlux). A set of statistical indicators were applied to evaluate the performance of each model. A linear regression analysis (y=ax+b) was applied to obtain the (a) slopes and (b) intercept variables; moreover, a residual analysis was performed to see if there were atypical values that affect the models. According to Chai and Draxler (2014), it is a good practice to include mean absolute error (MAE) and root mean square error (RMSE), because they are indicators that integrate the main differences between observed and estimated values. The variance (Sd2) was calculated to know 60 how much difference there was between observed and predicted values. The mean bias error (MBE) was included to find if there was a systematic error. The consistent error between the distance of linear regression and the 1:1 line is known as systematic RMSE (RMSEs). Unsystematic RMSE (RMSEu) is when the error is randomized, caused by an unknown source. When an unsystematic RMSE has low values, and the systematic RMSE value is close to RMSE, the model can be considered valid (Willmott et al., 1985). The efficiency model (EF) was applied by using the predicted and observed measured variations (Greenwood, Neeteson, & Draycott, 1985; Nash & Sutcliffe, 1970). Finally, an agreement index (d) (Willmott, 1981, 1982; Willmott & Wicks, 1980) was estimated for comparing between hydrological models. Lower is better Lower is better closer to 0, better closer to 0, better closer to 1, better 61 closer to 1, better where Oi +is the observed value (S-Penman) in the record i, Pi is the predicted value from the REEM and EEFlux models in area i, N is the number of observations (7), and n is the number of season days (256). Furthermore, P'i and O'i were obtained as Results The plotted results of E (S-Penman), mean E from the REEM and the EEFlux for the Bustillos Lagoon are shown in Figure 14. Figure 14. Evaporation values of S-Penman, REEM, and EEFlux during the 2017 agricultural season for the Bustillos Lagoon. Source: Rojas Villalobos with data retrieved from UNIFRUT (2019), USGS LandsatLook Viewer (2019), and EEFlux (2019). 62 According to Table 7, the EEFlux had significant variations at the beginning of the season on April 7 and April 23 (24.2% and -51.8%, respectively),as well as at the end of the cycle on June 17 and September 14 (-36.7% and -74.2%), while the REEM had sensitive variations on September 14 (13.6%). In the case of the REEM, the percentage variations represented a difference of less than 0.7 mm of evaporation except for May 9 and June 17, which were 0.98 and 1.11 mm, respectively. The EEFlux presented variations greater than 3.1 mm of evaporation for 3 of the seven days. For April 7, April 30, May 9 and May 16, the variations for the models tested were between 1.15 and 1.57 mm of the reference model. Table 7. Comparative table of errors between the reference evaporation and the models based on remote sensors (REEM and EEFlux). Source: Rojas Villalobos. Date DOY E Reference (mm) REEM EEFlux mm Error (%) mm Error (%) Apr-07-2017 96 6.3 6.6 4.0 7.9 24.2 Apr-23-2017 112 6.0 6.2 2.4 2.9 -51.8 Apr-30-2017 118 6.8 6.7 -1.5 5.6 -16.9 May-09-2017 128 8.4 7.4 -11.6 7.0 -16.5 May-16-2017 135 9.4 8.8 -6.6 7.8 -16.6 Jun-17-2017 167 9.0 7.9 -12.3 5.7 -36.7 Sep-14-2017 256 4.8 5.5 13.6 1.2 -74.2 Average -1.7 -26.9 63 Although the coefficient of determination (R2) was relatively high (0.953) to indicate that the REEM model produces evaporation values close to observed ones, , the slope (a = 0.6374) of the regression line does not ensure continuous linearity of predictions with the reference line. The intercept coefficient (b=2.3779) indicates overestimation of modelled data over observed values. The slope of the EEFlux (a=1.057) regression line closely matches the 1:1 reference of the observed data (S-Penman). Furthermore, the interception coefficient is negative (b = - 2.2123), and R2 is low (0.5105), which suggests an underestimation and high variance of the values predicted by the model. Both models concentrate on underestimation and overestimation values (EEFlux and REEM, respectively) in the range of 4.9 to 6.2 mm of daily evaporation. Figure 15 shows that the variance of the EEFlux model is not constant: while predicted evaporation values were low, the residual values were atypically high. In the residual analysis, evaporation is related to time. In other words, in April and September, the net radiation and temperatures were low, and as a result, there was less evaporation than that determined between May and August. When comparing the residuals between the two ET models, the REEM errors concentrate on the strip of ± 0.55 mm, which is quite acceptable, while more than 50% of the EEFlux residuals approximately exceed the range of ±1.37 mm and ±3.5 mm. 64 Figure 15. Comparative graphic of residuals predicted E on RS map models versus observed E (S-Penman). Source: Rojas Villalobos. The regression and residual analysis did not provide enough information to measure and compare the performance of the models studied. A more in-depth analysis was required for determining substantial differences between the comparison of the data of the predictive models with the reference ones. Table 8 shows the ranked analytical results for comparing the performance of ET models. For statistical indexes in complex evaluation systems, a weighting coefficient separately calculated is required. 65 Table 8. Summary of the ranked results of the comparative statistical indicators applied to the REEM and EEFlux versus S-Penman. Source: Rojas Villalobos. Index REEM (rank) EEFlux (rank) MAE (mm d-1) 0.55 (1) 2.23 (2) RMSE (mm d-1) -0.66 (1) 2.43 (2) Sd 2 (mm d-1) 0.44 (1) 3.14 (2) MBE (mm d-1) -0.25 (1) -1.79 (2) RMSEu (mm d -1) 0.60 (1) 4.26 (2) RMSEs (mm d -1) 0.41 (1) 2.14 (2) EF 0.82(1) -1.36 (2) R2 0.95 (1) 0.51 (2) d 0.94 (1) 0.68 (2) a (intercept) 1.75 (2) 1.62 (1) b (slope) 0.79 (2) 1.07 (1) The RMSE has been criticized for being inappropriate and misinterpreted in environmental and climate analyses (Willmott & Matsuura, 2005), but the results of the RMSE and MAE enrich the interpretation of the evaluated models (Chai & Draxler, 2014). In this study, the MAE and RMSE indicators agreed that the REEM presented a lower average error (MAE = 0.55 and RMSE = -0.66 both in mm d-1) among the data. Sd2 confirms the high variability that the EEFlux had (3.14 mm) in predicting the daily ETa in comparison to the REEM (0.44 mm). The bias indicator (MBE) agreed with the initial linear regression analysis as it showed a slight underestimation of the values calculated by REEM (-0.25 mm) in comparison with the higher underestimation of the values predicted by the EEFlux (-1.79 mm). The RMSEu results suggested that noise from an unknown source promoted a poor performance of the EEFlux model (4.26 mm). In contrast, the same index showed a lower influence of unknown variables in the REEM model (0.60 mm). According to the EF index, values close to 1 correspond to a model that predicts values close to the observed data. If the 66 index is less than 0, the mean observed data is a better predictor than the values estimated from the ET model (Nash & Sutcliffe, 1970; Pushpalatha, Perrin, Moine, & Andréassian, 2012). Therefore, according to the above, REEM (EF=0.82) had a higher performance than EEFlux (EF= -1.36). The statistical indicator of agreement "d" indicates the tendency of the previous indexes by suggesting that the REEM (0.94) is a better predictor of ETa than the EEFlux (0.68). The total E for the three models in the agricultural reference season was compared using daily estimations. In the case of the REEM and EEFlux, a linear interpolation technique was used to calculate the E between the dates of the seven available satellite images. The meteorological records of the aforementioned agroclimatic station were used for the computation of the daily E- reference through the S-Penman equation (Figure 16). The variability (SEE) was 3.2- and 3.4- mm day-1 for REEM and EEFlux, respectively. The total E for S-Penman was 968 mm, and 1137 mm for REEM, with 752 mm for EEFlux, which is equivalent to 115.29, 135.35, and 752 hm3 of water, respectively. 67 Figure 16. Seasonal evaporation comparison of RS models versus S-Penman data from April 4, 2017 to September 14, 2017. Source: Rojas Villalobos. Discussion Statistical results suggested a better predictive performance of the evaporation of water bodies of the REEM model versus the EEFlux model for the 2017 agricultural cycle. The residue analysis showed more considerable variability in the low ranges of E reference. This variability may be induced by solar radiation, air temperature, relative humidity, and wind because they are weather variables that have a strong influence on ET (Valipour, 2015). The METRIC model uses these variables to estimate H, employing the alfalfa reference ET by using the ASCE Penman-Monteith equation, and the model assumes that the cold pixel has a sensible heat flux (H) equal to zero. The REEM uses the same local variables by employing regression equations to calculate H and Rn. Figure 17 displays an ET and E comparison map of the REEM and EEFlux from agricultural fields and the Bustillos Lagoon (dated June 17, 2017). 68 Figure 17. ET (crop fields) and evaporation (lagoon) comparison maps of REEM and EEFlux models in the Cuauhtemoc Valley for June 17, 2017. Source: Rojas Villalobos with data retrieved from USGS (2019) and EEFlux (2019). The first difference between the application of METRIC within the EEFlux was that in EEFlux gridded weather data sets were used instead of climate data from the field. Point data such as from agroclimatic stations and interpolated data such as gridded data sets have 69 significant spatial differences. Although the interpolation models used to generate gridded weather data sets have improved, there is still a degree of uncertainty because of the distance between the meteorological reference stations. For instance, the Global Land Data Assimilation System (GLDAS-1), the North American Land Data Assimilation System (NLDAS-2), the Climate Forecast System Version 2 (CFSv2), GridMET, the Real-Time Mesoscale Analysis (RTMA) and the National Digital Forecast Database (NDFD) are gridded data sets with spatial resolution ranges between 4 to 12 km (Allen et al., 2015). Regarding gridded data, Blankenau (2017) found that there were biases and inconsistencies in the gridded climatic data potentially caused by the distances and the location of the interpolated points. The databases were built using weather stations located at airports, which do not represent the weather conditions of an agricultural area (colder and wetter). In addition, ET underestimations occurred because the gridded data did not integrate the effects of humidification and cooling near the surface when agricultural fields were irrigated. Since atmospheric conditions vary during the day, instantaneous weather data were obtained through linear regression from hourly values according to the time when the satellite passed over the study area. If the instantaneous data were generated from a large spatial resolution grid that integrates biases and errors, the uncertainty was propagated to the predicted data (ET) (Kauffeldt, Halldin, Rodhe, Xu, & Westerberg, 2013; Lobell, 2013; Phillips & Marks, 1996). The daily evaporation variability of the RS models and the value measured in the season was high since the coefficient of variation was 53.6% for REEM and 55.7% for EEFlux. The 70 daily E variability between the RS models and the value measured in the season was high since the coefficient of variation was 53.6 % for REEM and 55.7% for EEFlux. Similarly, the REEM overestimated E by 17.4 % when compared to reference values, while EEFlux underestimated E by 22.3% for the same period. In the segment from May 16 to September 14 (135-256 DOY), there were significant differences in the coefficients of variation when REEM obtained 70% and EEFlux 46%. The differences between the predicted values and the observed values were particularly high because of the large gaps between the dates of the acquired satellite images. Conclusions In this study, seven Landsat 8 images were used during the agricultural cycle from April to September 2017, when the REEM and EEFlux evapotranspiration models were compared with the reference ET to estimate the daily evaporation of the Bustillos Lagoon. ET estimation methods by remote sensors are sensitive to variations in weather conditions. In the interpolated grid of climatic parameters, there are regions where there are significant differences between observed and interpolated data. These regions are far from the interpolation source points, and the physical-environmental conditions are different. Gridded data should aggregate additional data source points where there are significant variations of the climatic parameters. An anchor weather station can improve the predictions of the evaporation of a water body as observed in the REEM model. The location of the weather station is a determining factor in the computation of the ET. In this study, an agroclimatic station located 4.5 km from the Bustillos Lagoon recorded weather conditions where the prevailing winds (SW-NW) pass before reaching the lake, which establishes the physical conditions for water evaporation. The temporal resolution of the satellite 71 scenes is a determining factor for the estimation of the total E since the gap between the dates of the images reduce data time uncertainty in order to obtain accurate values and a better performance of the RS models through interpolation methods. 72 References AghaKouchak, A., Norouzi, H., Madani, K., Mirchi, A., Azarderakhsh, M., Nazemi, A., … Hasanzadeh, E. (2015). Aral Sea syndrome desiccates Lake Urmia: Call for action. Journal of Great Lakes Research, 41(1), 307–311. doi: 10.1016/j.jglr.2014.12.007 Allen, R. G., Morton, C., Kamble, B., Kilic, A., Huntington, J., Thau, D., … Robison, C. (2015). EEFlux: A Landsat-based Evapotranspiration mapping tool on the Google Earth Engine. 2015 ASABE / IA Irrigation Symposium: Emerging Technologies for Sustainable Irrigation - A Tribute to the Career of Terry Howell, Sr. Conference Proceedings, 1–11. doi: 10.13031/irrig.20152143511 Allen, R. G., Tasumi, M., & Trezza, R. (2007). Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC)-Model. 133(4), 135–139. doi: 10.1061/(ASCE)0733-9437(2007)133:4(380) Allen, R. G., Tasumi, M., Trezza, R., Wright, J. L., Bastiaanssen, W., Kramber, W. J., … Robison, C. (2007). Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC)—Applications. Journal of Irrigation and Drainage Engineering, 133(4), 395–406. doi: 10.1061/(ASCE)0733-9437(2007)133:4(395) Álvarez, J. P. A., Cutillas, P. P., Valle, O. R., & Cabañero, J. J. A. (2016). Water resources degradation in a semiarid environment. Bustillos and Los Mexicanos lagoons (Chihuahua, México). Papeles de Geografía, (62), 107–118. Amado-Álvarez, J. P., Pérez Cutillas, P., Ramírez Valle, O., & Alarcón Cabañero, J. J. (2016). Degradación de los recursos hídricos en un ambiente semiárido. Las lagunas de Bustillos y de los Mexicanos (Chihuahua, México). Retrieved from https://digitum.um.es/digitum/handle/10201/51756 Amado-Alvarez, J., Pérez-Cutillas, P., Alatorre-Cejudo, L. C., Ramírez-Valle, O., Segovia- Ortega, E. F., & Alarcón-Cabañero, J. J. (2019). Multispectral analysis to estimate the turbidity as an indicator of the quality of water in reservoirs in Chihuahua State, Mexico. Revista Geográfica de América Central, 1(62), 49–77. doi: 10.15359/rgac.62-1.2 ASCE–EWRI. (2005). The ASCE standardized reference evapotranspiration equation. Reston, VA: ASCE–EWRI Standardization of Reference Evapotranspiration Task Committee Rep. Ayyad, S., Al Zayed, I. S., Ha, V. T. T., & Ribbe, L. (2019). The Performance of Satellite-Based Actual Evapotranspiration Products and the Assessment of Irrigation Efficiency in Egypt. Water, 11(9), 1913. doi: 10.3390/w11091913 Bastiaanssen, W. G. M. (1995). Regionalization of surface flux densities and moisture indicators in composite terrain: A remote sensing approach under clear skies in Mediterranean climates (SC-DLO). Retrieved from https://library.wur.nl/WebQuery/wurpubs/28279 Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., & Holtslag, A. A. M. (1998). A remote sensing surface energy balance algorithm for land (SEBAL). Part 1. Formulation. Journal of Hydrology, 212–213, 198–212. doi: 10.1016/S0022-1694(98)00253-4 73 Bastiaanssen, W. G. M., Pelgrum, H., Wang, J., Ma, Y., Moreno, J. F., Roerink, G. J., & Van der Wal, T. (1998). A remote sensing surface energy balance algorithm for land (SEBAL). Part 2: Validation. Journal of Hydrology, 212–213, 213–229. doi: 10.1016/S0022- 1694(98)00254-6 Blaney, H. F., & Criddle, W. D. (1957). Determining consumptive use and irrigation water requirements (Technical Bulletin No. 1275; p. 24). Retrieved from United States Department of Agriculture website: https://naldc.nal.usda.gov/download/CAT87201264/PDF Bozorgi, A., Bozorg-Haddad, O., Sima, S., & Loáiciga, H. A. (2018). Comparison of methods to calculate evaporation from reservoirs. International Journal of River Basin Management, 0(0), 1–12. doi: 10.1080/15715124.2018.1546729 Brutsaert, W., & Stricker, H. (1979). An advection-aridity approach to estimate actual regional evapotranspiration. Water Resources Research, 15(2), 443–450. doi: 10.1029/WR015i002p00443 Cabrera, M. C. M., Anache, J. A. A., Youlton, C., & Wendland, E. (2016). Performance of evaporation estimation methods compared with standard 20 m2 tank. 20(10), 874–879. doi: doi.org/10.1590/1807-1929/agriambi.v20n10p874-879 Chai, T., & Draxler, R. R. (2014). Root mean square error (RMSE) or mean absolute error (MAE)?–Arguments against avoiding RMSE in the literature. Geoscientific Model Development, 7(3), 1247–1250. doi: 10.5194/gmd-7-1247-2014 Cleugh, H. A., Leuning, R., Mu, Q., & Running, S. W. (2007). Regional evaporation estimates from flux tower and MODIS satellite data. Remote Sensing of Environment, 106(3), 285– 304. doi: 10.1016/j.rse.2006.07.007 De Bruin, H. a. R. (1978). A Simple Model for Shallow Lake Evaporation. Journal of Applied Meteorology, 17(8), 1132–1134. doi: 10.1175/1520- 0450(1978)017<1132:ASMFSL>2.0.CO;2 De Bruin, H. a. R., & Keijman, J. Q. (1979). The Priestley-Taylor Evaporation Model Applied to a Large, Shallow Lake in the Netherlands. Journal of Applied Meteorology, 18(7), 898– 903. doi: 10.1175/1520-0450(1979)018<0898:TPTEMA>2.0.CO;2 Diario Oficial de la Federación. (2015, July 6). ACUERDO por el que se da a conocer el resultado de los estudios técnicos de aguas nacionales subterráneas del Acuífero Cuauhtémoc, clave 0805, en el Estado de Chihuahua, Región Hidrológico Administrativa Río Bravo. Retrieved from http://www.dof.gob.mx/nota_detalle.php?codigo=5399497&fecha=06/07/2015 Díaz Caravantes, R. E., Bravo Peña, L. C., Alatorre Cejudo, L. C., & Sánchez Flores, E. (2014). Análisis geoespacial de la interacción entre el uso de suelo y de agua en el área peri- urbana de Cuauhtémoc, Chihuahua. Un estudio socioambiental en el norte de México. Investigaciones Geográficas, Boletín Del Instituto de Geografía, 2014(83), 116–130. doi: 10.14350/rig.32694 Erickson, A. J., Gulliver, J. S., Hozalski, R. M., Mohseni, O., Nieber, J. L., Wilson, B. N., & Weiss, P. T. (2008). Evaporation and Evapotranspiration | Stormwater Treatment: 74 Assessment and Maintenance. In Assessment of Stormwater Best Management Practices (p. 112). Retrieved from https://conservancy.umn.edu/bitstream/handle/11299/181358/cfans_asset_115795.pdf?se quence=1 Fisher, J. B., Melton, F., Middleton, E., Hain, C., Anderson, M., Allen, R., … Wood, E. F. (2018). The future of evapotranspiration: Global requirements for ecosystem functioning, carbon and climate feedbacks, agricultural management, and water resources. Reviews of Geophysics, 2618–2626. doi: 10.1002/2016WR020175@10.1002/(ISSN)1944- 9208.COMHES1 Fu, G., Charles, S. P., & Yu, J. (2009). A critical overview of pan evaporation trends over the last 50 years. Climatic Change, 97(1), 193. doi: 10.1007/s10584-009-9579-1 Greenwood, D. J., Neeteson, J. J., & Draycott, A. (1985). Response of potatoes to N fertilizer: Dynamic model. Plant and Soil, 85(2), 185–203. doi: 10.1007/BF02139623 Gross, M. (2017). The world’s vanishing lakes. Current Biology, 27(2), R43–R46. doi: 10.1016/j.cub.2017.01.008 Hamon, W. R. (1960). Estimating potential evapotranspiration (Thesis, Massachusetts Institute of Technology). Retrieved from https://dspace.mit.edu/handle/1721.1/79479 Hewitt, I. C., Fernald, A. G., & Samani, Z. A. (2018). Calculating Field-Scale Evapotranspiration with Closed-Chamber and Remote Sensing Methods. JAWRA Journal of the American Water Resources Association, 54(4), 962–973. doi: 10.1111/1752- 1688.12654 Hirschi, M., Michel, D., Lehner, I., & Seneviratne, S. I. (2017). A site-level comparison of lysimeter and eddy covariance flux measurements of evapotranspiration. Hydrology and Earth System Sciences, 21(3), 1809–1825. doi: https://doi.org/10.5194/hess-21-1809- 2017 Irmak, A., Allen, R. G., Kjaersgaard, J., Huntington, J., Kamble, B., Trezza, R., & Ratcliffe, I. (2012). Operational Remote Sensing of ET and Challenges. Evapotranspiration - Remote Sensing and Modeling. doi: 10.5772/25174 Jensen, M. E., & Haise, H. R. (1963). Estimating Evapotranspiration from Solar Radiation. Proceedings of the American Society of Civil Engineers, Journal of the Irrigation and Drainage Division, 89, 15–41. Kauffeldt, A., Halldin, S., Rodhe, A., Xu, C.-Y., & Westerberg, I. K. (2013). Disinformative data in large-scale hydrological modelling. Hydrology and Earth System Sciences, 17, 2845– 2857. doi: 10.5194/hess-17-2845-2013 Kirk, J. T. (1985). Effects of suspensoids (turbidity) on penetration of solar radiation in aquatic ecosystems. Hydrobiologia, 125(1), 195–208. doi: https://doi.org/10.1007/BF00045935 Kıvrak, C., Bawazir, S., Samani, Z., Steele, C., & Sönmez, B. (2019). Using Plant Phenology and Landsat-8 Satellite Data to Quantify Water Use by Onion Crop in the Mesilla Valley, New Mexico. International Journal of Crop Science and Technology, 5(1), 1–13. doi: 10.26558/ijcst.427945 75 Li, F., Cao, R., Zhao, Y., Mu, D., Fu, C., & Ping, F. (2017). Remote sensing Penman–Monteith model to estimate catchment evapotranspiration considering the vegetation diversity. Theoretical and Applied Climatology, 127(1), 111–121. doi: 10.1007/s00704-015-1628-2 Lobell, D. B. (2013). Errors in climate datasets and their effects on statistical crop models. Agricultural and Forest Meteorology, 170, 58–66. doi: 10.1016/j.agrformet.2012.05.013 Makkink, G. F. (1957). Testing the Penman Formula by Means of Lysimeters. Journal of the Institution of Water Engineers, (11), 277–288. Martínez Pérez, J. Á., García-Galiano, S. G., Martin-Gorriz, B., & Baille, A. (2017). Satellite- Based Method for Estimating the Spatial Distribution of Crop Evapotranspiration: Sensitivity to the Priestley-Taylor Coefficient. Remote Sensing, 9(6), 611. doi: 10.3390/rs9060611 Mireles, C., & Mellink, E. (2017). Use of Laguna de Bustillos, Chihuahua, by Waterbirds during the 2011–2012 Wintering Season. The American Midland Naturalist, 178(1), 82–96. doi: 10.1674/0003-0031-178.1.82 Monin, A., & Obukhov, A. (1954). Basic laws of turbulent mixing in the atmosphere near the ground. Tr. Akad. Nauk SSSR Geofiz. Inst, 24(151), 163–187. Nash, J. E., & Sutcliffe, J. V. (1970). River flow forecasting through conceptual models part I — A discussion of principles. Journal of Hydrology, 10(3), 282–290. doi: 10.1016/0022- 1694(70)90255-6 Okpara, U. T., Stringer, L. C., Dougill, A. J., & Bila, M. D. (2015). Conflicts about water in Lake Chad: Are environmental, vulnerability and security issues linked? Progress in Development Studies, 15(4), 308–325. doi: 10.1177/1464993415592738 Papadakis, J. (1965). Potential Evapotranspiration. Soil Science, 100, 76. doi: 10.1097/00010694-196507000-00039 Penman, H. L. (1948). Natural evaporation from open water, bare soil and grass. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 193(1032), 120–145. doi: 10.1098/rspa.1948.0037 Phillips, D. L., & Marks, D. G. (1996). Spatial uncertainty analysis: Propagation of interpolation errors in spatially distributed models. Ecological Modelling, 91(1), 213–229. doi: 10.1016/0304-3800(95)00191-3 Priestley, C. H. B., & Taylor, R. J. (1972). On the Assessment of Surface Heat Flux and Evaporation Using Large-Scale Parameters. Monthly Weather Review, 100(2), 81–92. doi: 10.1175/1520-0493(1972)100<0081:OTAOSH>2.3.CO;2 Pushpalatha, R., Perrin, C., Moine, N. L., & Andréassian, V. (2012). A review of efficiency criteria suitable for evaluating low-flow simulations. Journal of Hydrology, 420–421, 171–182. doi: 10.1016/j.jhydrol.2011.11.055 Rayner, D. P. (2007). Wind Run Changes: The Dominant Factor Affecting Pan Evaporation Trends in Australia. Journal of Climate, 20(14), 3379–3394. doi: 10.1175/JCLI4181.1 Reitz, M., Senay, G., & Sanford, W. E. (2017). Combining remote sensing and water-balance evapotranspiration estimates for the conterminous United States. Remote Sensing, 9(12), 117. doi: 10.3390/rs9121181 76 Rohwer, C. (1931). Evaporation from free water surfaces (No. 271; p. 107). Retrieved from United States Department of Agriculture website: https://naldc.nal.usda.gov/download/CAT86200265/PDF Rojas-Villalobos, H. L., Alatorre-Cejudo, L. C., Stringman, B., Samani, Z., & Brown, C. (2018). Topobathymetric 3D model reconstruction of shallow water bodies through remote sensing, GPS, and bathymetry. TECNOCIENCIA Chihuahua, 12(1), 42–54. doi: 10.5281/zenodo.3341734 Rooney, G. G., & Bornemann, F. J. (2013). The performance of FLake in the Met Office Unified Model. Tellus A: Dynamic Meteorology and Oceanography, 65(1), 21363. doi: 10.3402/tellusa.v65i0.21363 Samani, A., & Bawazir, S. (2015). Improving evapotranspiration estimation using remote sensing technology (TECHNICAL COMPLETION REPORT No. 125548; p. 28). Retrieved from https://nmwrri.nmsu.edu/wp- content/SWWA/Reports/Bawazir/June2015FINAL/report.pdf Samani, Z., Bawazir, A. S., Skaggs, R. K., Bleiweiss, M. P., Piñon, A., & Tran, V. (2007). Water Use by Agricultural Crops and Riparian Vegetation: An Application of Remote Sensing Technology. Journal of Contemporary Water Research & Education, 137(1), 8–13. doi: 10.1111/j.1936-704X.2007.mp137001002.x Samani, Z., Bawazir, S., Bleiweiss, M., Skaggs, R., Piñon, A., & Tran, V. D. (2007). Estimating Daily Net Radiation over Vegetation Canopy through Remote Sensing and Climatic Data. Journal of Irrigation and Drainage Engineering, 133(4), 291–297. doi: 10.1061/(ASCE)0733-9437(2007)133:4(291) Samani, Z., Bawazir, S., Bleiweiss, M., Skaggs, R., & Schmugge, T. (2006). Estimating riparian ET through remote sensing. Presented at the AGU Fall Meeting, San Francisco. Samani, Z., Sammis, T., Skaggs, R., Alkhatiri, N., & Deras, J. (2005). Measuring On-Farm Irrigation Efficiency with Chloride Tracing under Deficit Irrigation. Journal of Irrigation and Drainage Engineering, 131(6), 555–559. doi: 10.1061/(ASCE)0733- 9437(2005)131:6(555) Samani, Z., Skaggs, R., & Bleiweiss, M. (2005). Regional ET Estimation Model, REEM. Proc. 31st Intl. Symposium on Remote Sensing from Satellite. Presented at the Proc. 31st Intl. Symposium on Remote Sensing from Satellite, Tucson, AR. Senkondo, W., Munishi, S. E., Tumbo, M., Nobert, J., & Lyon, S. W. (2019). Comparing Remotely-Sensed Surface Energy Balance Evapotranspiration Estimates in Heterogeneous and Data-Limited Regions: A Case Study of Tanzania’s Kilombero Valley. Remote Sensing, 11(11), 1289. doi: 10.3390/rs11111289 Smith, R. C., & Tyler, J. E. (1967). Optical Properties of Clear Natural Water*. JOSA, 57(5), 589–595. doi: 10.1364/JOSA.57.000589 Stephens, J. C., & Stewart, E. H. (1963). A comparison of procedures for computing evaporation and evapotranspiration (pp. 123–133). Fort Lauderdle, Florida. Subin, Z. M., Murphy, L. N., Li, F., Bonfils, C., & Riley, W. J. (2012). Boreal lakes moderate seasonal and diurnal temperature variation and perturb atmospheric circulation: Analyses 77 in the Community Earth System Model 1 (CESM1). Tellus A: Dynamic Meteorology and Oceanography, 64(1), 15639. doi: 10.3402/tellusa.v64i0.15639 UNIFRUT. (2019). Resumen de las estaciones meteorológicas por municipios. Retrieved February 20, 2019, from http://unifrut.com.mx/rem/ Valiantzas, J. D. (2006). Simplified versions for the Penman evaporation equation using routine weather data. Journal of Hydrology, 331(3), 690–702. doi: 10.1016/j.jhydrol.2006.06.012 Valipour, M. (2015). Importance of solar radiation, temperature, relative humidity, and wind speed for calculation of reference evapotranspiration. Archives of Agronomy and Soil Science, 61(2), 239–255. doi: 10.1080/03650340.2014.925107 Wang, B., Ma, Y., Ma, W., Su, B., & Dong, X. (2019). Evaluation of ten methods for estimating evaporation in a small high-elevation lake on the Tibetan Plateau. Theoretical and Applied Climatology, 136(3), 1033–1045. doi: 10.1007/s00704-018-2539-9 Wang, H., Tetzlaff, D., & Soulsby, C. (2017). Testing the maximum entropy production approach for estimating evapotranspiration from closed canopy shrubland in a low- energy humid environment. Hydrological Processes, 31(25), 4613–4621. doi: 10.1002/hyp.11363 Willmott, C. J. (1981). On the validation of models. Physical Geography, 2(2), 184–194. doi: doi.org/10.1080/02723646.1981.10642213 Willmott, C. J. (1982). Some comments on the evaluation of model performance. Bulletin of the American Meteorological Society, 63(11), 1309–1313. doi: https://doi.org/10.1175/1520- 0477(1982)063<1309:SCOTEO>2.0.CO;2 Willmott, C. J., Ackleson, S. G., Davis, R. E., Feddema, J. J., Klink, K. M., Legates, D. R., … Rowe, C. M. (1985). Statistics for the evaluation and comparison of models. Journal of Geophysical Research: Oceans, 90(C5), 8995–9005. doi: https://doi.org/10.1029/JC090iC05p08995 Willmott, C. J., & Matsuura, K. (2005). Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance. Climate Research, 30(1), 79–82. Retrieved from JSTOR. Willmott, C. J., & Wicks, D. E. (1980). An empirical method for the spatial interpolation of monthly precipitation within California. Physical Geography, 1(1), 59–73. doi: https://doi.org/10.1080/02723646.1980.10642189 Xu, H. (2006). Modification of normalized difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14), 3025–3033. doi: doi.org/10.1080/01431160600589179 Zhang, K., Kimball, J. S., & Running, S. W. (2016). A review of remote sensing based actual evapotranspiration estimation. Wiley Interdisciplinary Reviews: Water, 3(6), 834–853. doi: 10.1002/wat2.1168 Zhu, W., Jia, S., & Lv, A. (2019). A Statistical Analysis of the Remotely Sensed Land Surface Temperature–Vegetation Index Method for the Retrieval of Evaporative Fraction Over Grasslands in the Southern Great Plains. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 1–8. doi: 10.1109/JSTARS.2019.2917183 78 Chapter 4 Single-input, multiple-output iterative algorithm for the calculation of volume, area, elevation, and shape using 3D topobathymetric models. Article submitted: November 19, 2019. Status: In review. Journal: Investigaciones Geográficas (UNAM). Abstract Most methods for estimating the morphometric values of water bodies use equations derived from hypsographic curves or digital terrain models (DTMs) that relate depth, volume (V), and area (A) and that model the uncertainty inherent in the complex underwater morphology. This research focuses directly on the use of topobathymetric models that include the bathymetry and topography of the surrounding area next to the water body. The projection of the water surface height (H) on each DTM pixel generates a water column with intrinsic attributes such as volume and area. The process is replicated among all cells and estimates the total area and volume of the water body. If the V or A is the input data, an algorithm that iterates height values is used to generate the new data, which is compared with the entered value that functions as a reference. If the difference between the reference value and the calculated value is less than an error threshold, the iteration stops, and the maximum and average depths are calculated. In addition, the raster and the shape that represent the body of water are created. The cross comparison of H-V-A showed that there is an error between 0.0034% and 0.000039% when any of the parameters are used as input data. Performance tests determined that pixel dimensions are directly proportional to the processing time for each iteration. The results of the 79 implementation of this algorithm were satisfactory since, for the DTM of Bustillos Lagoon, Chihuahua, Mexico, the simulation took less than 17 seconds in at most 22 iterations. Introduction Calculating physical characteristics of water bodies such as volume, surface area, and shape is a challenging process because of the complex underwater topography. The water bodies floor is usually irregular, with elevations and depressions that do not follow a specific pattern and therefore are difficult to model with mathematical equations. The height of the water level (H), volume (V), surface area (A), and the shape of the surface area are parameters that are not linearly interrelated. Determining these parameters, using a known value from the previously mentioned parameters (H, V, or A), will allow erosion modelers, hydrologists, geohydrologists, and ecologists, among others, to use morphological parameters in their simulations. Currently, geographic information systems (GIS) are able to deploy programming languages to develop tools that respond to complex problems. This approach generates information about the storage capacity in a dynamic way to support management policies dealing with flood risk zones and minimum water levels for maintaining ecologically healthy areas and other water resource issues. Several methodologies exist to calculate morphometric parameters, such as height, volume, area, maximum depth, and the average depth of water bodies. The first studies that relate volume-area-depth used the radius between average depth and maximum depth in addition to sinusoidal parameters (such as lake bottom profile) to do so (Lehman, 1975; Neumann, 1959). Sima and Tajrishy (2013) presented a model that relates volume-area-elevation using data from 80 remote sensors and analytical procedures such as a power model and a truncated pyramid model. This approach results in highly parameterized equations that relate morphometric values. Moreno-Amich and Garcia-Berthou (1989) used echo sounding to relate morphometric characteristics of depth-area measurements for developing hypsographic curves and generating a bathymetric map. Johansson et al. (2007) proposed two new mathematical models that interrelate morphometric values: the volume development, which is an equation based on the A-V relationship curve (Vd) and the hypsographic development parameter, which is the integration of A-depth and V-depth relation curves (Hd). These models require three inputs: V, maximum depth, and A. Recent methodologies that use autonomous aerial vehicles measure the height of the terrain through LIDAR (Laser Imaging Detection and Ranging) and surface water vehicles that measure the depth (bathymetry) through high-resolution echo sounding. These data sources are processed in GIS and generate accurate DTMs (Erena, Atenza, García-Galiano, Domínguez, & Bernabé, 2019). Regardless of the methodology, however, the equations that relate the morphometry variables inherit the uncertainty of the complexity and spatial variability of underwater topography (Rode et al., 2010). Chen et al. (2018) presented a method that uses 3D geometry of a dam with which the volume and floodplains are calculated. The algorithm uses precipitation and water stage values as input data. It then simulates the floods in two sections: the floodplain and the 3D model from which the morphological parameters are obtained. Despite being an efficient model, it is limited as to what input data it will accept. Thus, none of the 81 methods shown above are capable of offering solutions where morphological values interact with each other to respond to the needs of hydrologists. To address these uncertainties, this study develops a technique that fully estimates the H- V-A of water bodies using computational techniques through 3D models that include bathymetry and the surrounding terrain topography. This technique used the water column below the level of the water surface projected onto each pixel of the DTM; this calculation was applied to the entire study area to delineate the extension of the water body. The V or A was the reference variable deployed in an iterative algorithm until the error threshold was met. Study area The Bustillos Lagoon is in the endorheic Cuauhtémoc Basin, and the lagoon has an area of 3,298.15 km2. A mountain range, called Sierra Azul, surrounds it in the north-northeast; on the western flank are Mennonite colonies where the terrain slope is below 1%. In addition, the town of Anahuac is in the south. The Bustillos Lagoon is between the quadrant coordinates 28 ° 38 '51' 'N - 28 ° 28' 27 '' N and 106 ° 57 '3' 'W - 106 ° 38' 50 '' W in the Cuauhtémoc municipality in the Mexican state of Chihuahua (Figure 18). 82 Figure 18. Study area where the algorithm was applied. The Bustillos Lagoon in Chihuahua. Source: Rojas Villalobos with data retrieved from INEGI (2019). Because the Cuauhtémoc Basin is between the semi-humid climate of the Sierra Madre Occidental and the Chihuahuan Desert, the region's weather is warm and semi-dry (Kottek, Grieser, Beck, Rudolf, Rubel, 2006). The approximate elevation of the Cuauhtémoc region is 2,100 m above sea level (m.a.s.l), and the average annual temperature ranges from 6.9 to 21 ° C, with an average annual rainfall of about 528 mm y-1 (Servicio Meteorológico Nacional, 2019). Material and methods This methodology employed a personal computer with an Intel i3-8100 3.6 GHz processor, 48 GB RAM, and two SSD of 1 TB each. For this computational development in GIS, it is essential to have a DTM that includes the bathymetry of the study area for generating the 83 hydrological characteristics of the water body - the 3D topobathymetric model of the Bustillos Lagoon (spatial resolution of 5 meters)(Rojas-Villalobos, Alatorre-Cejudo, Stringman, Samani, & Brown, 2018) (Figure 19). Figure 19. DTM of the Bustillos Lagoon. Source: Rojas Villalobos with data retrieved from Rojas-Villalobos et al. (2018). The software used for GIS processing was ArcMap® version 10.6 of Environmental Systems Research Institute, ESRI (ArcMap, 2019). The 3D process tool called Surface Volume, which requires the DTM and the height of the water level as input data, performed V and A calculations (numerical results), and Python® version 2.7.13 was the language to encode the algorithm (Python Language Reference, 2019). The algorithm can capture one of the various input data, and as a result, it can generate the rest of the output data; for this reason, the algorithm is cataloged as a single-input, multiple- output data algorithm. The algorithm uses one of the following input values: the height of the 84 water level, the storage V, or the A of the water body (Figure 20). All calculations are in the International System of Units. Figure 20. The schematic diagram shows single-inputs and multiple-output data for iterative algorithm. Source: Rojas Villalobos. The algorithm was designed using the following criteria. The algorithm is divided into two sections and depends on the input data: i) water height in meters above sea level (m.a.s.l) and ii) V (m3) or A (m2). Some of the process used in the second section refers to procedures in the first section. The Surface Volume tool used the height value and the DTM to calculate the V and A of the water body. These two results were used to obtain the maximum and average depth (Figure 21) Figure 21. Equations to calculate Average Depth and Maximum Depth. 85 where AD is average depth (m), V is the volume (m3), A is the surface area of the water body (m2), MD is maximum depth (m), H is the height of water level (m), and Hmin is the height of the water body floor (m), which was extracted from the properties of the raster. The Map Algebra tool, included in ArcMap®, was used to perform the extraction of the raster that represents the water surface filtering of all pixels that were less than or equal to H. A conversion tool then saves the raster of surface water as a polygon (vector data) in a shapefile format or geodatabase. When the second process starts, the user captures (or by default) an error threshold (ET)(%) that is required to stop the iteration process and the value of V or A used to compute output information. V or A assumes the value of reference (Ref) that is used to compare with the new iterated values (V or A). The threshold limit (TL) is the value that stops the iteration and is the product of V or A, multiplied by ET. Three initial variables were as follows: Step equal to 1 used to increase or decrease H, H equal to 1 meter above the bottom of the lagoon, and direction (Dir) equal to “upward.” The iteration starts with the H and the Surface Volume tool that calculates new data (ND = V or A). If the absolute value of the difference between Ref and ND is less than TL, the algorithm proceeds to execute the procedures for calculating the output information such as raster image and polygon shape of the lagoon. Otherwise, H continues increasing and generating ND until the absolute difference between Ref and ND is less than TL (e.g., 50 m3 or 70 m2). If ND surpasses Ref, H decreases at the halfway point of the previous step (e.g., 0.5 m.) until the absolute difference between ND and Ref is less than TL. If the TL is not accomplished and the ND surpasses Ref, H starts to increase with a new step (e.g., 0.25 m.). 86 This iteration stops when the absolute difference is less than TL, and the algorithm calculates output data. The algorithm diagram is shown in Figure 22. Figure 22. Flowchart of the iterative algorithm to compute hydrologic characteristics using single-input data. Source: Rojas Villalobos. 87 Results and discussion Table 9 shows the results of three simulations with different data input. For the second and third models, the data resulting from the first simulation were used as input data (V and A) for the cross comparison since the V-A estimates are calculated directly from the height, and it is not necessary to iterate data. Table 9. Result of the calculations of the implementation of the algorithm in Python language. Study site: the Bustillos Lagoon, Chihuahua, Mexico. Error threshold = 0.01%. * Input data. Source: Rojas Villalobos. Error Threshold Area km3 Height masl Volume hm3 Area km2 Average depth m Maximum depth m Iterations Processing time s 0.01 1973.7* 289.1004 114.256 2.5302 3.4860 0 1.407 0.01 1973.69 289.100* 114.255 2.5302 3.4859 17 13.668 0.01 1973.69 289.0111 114.25* 2.5295 3.4852 14 11.466 The results of the cross comparison of H-V-A showed that the differences are 0.003, 0.0034, and 0.000039%, respectively. These values are below the established error threshold of 0.01% and represent a height differential of less than one micrometer, which is negligible in the lagoon modeling scale. The DTM covers an area of 246.00 km2, which contains the entire lagoon and the surrounding area in a buffer greater than 1000 meters. The lagoon has a storage capacity of 360.52 hm3 and a surface area of 122.8 km2 before extending to the floodplains at 1974.3 m.a.s.l. The processing time depends directly on the number of pixels of the DTM used in the modeling and not on the lagoon area itself. Two DTMs were modeled with different pixel dimensions: DTMa) 5.879 (width) x 3.925 (height) corresponding to 576.8 km2 and DTMb) 3.198 (width) x 3.077 (height) equivalent to 246.00 km2. Each pixel maintains a spatial resolution of 5 meters. The number of iterations varies between 16 and 22 because of variations 88 in the calculated volume, which does not exceed the reference volume in each of the iterations. These variations can decrease or increase the number of iterations and, consequently, the calculation time (Table 10). The dimension of the DTM is directly and linearly related to the processing time in each iteration. The DTMa model is 2.34 times larger than DTMb, and this ratio is repeated in the average runtime of 1.84 seconds per iteration for the DTMa and 0.78 seconds per iteration for DTMb. This advantage can be exploited by hydrological modelers that require real-time results because they do not have to consider the simulation area but rather the number of pixels contained in the DTM. Table 10. Iterative model processing times with various storage volume input values using two DTMs with different pixel dimensions. Pixel spatial resolution: 5 meters. Source: Rojas Villalobos. DTM Tested Area (km3) Volume (m3) Iterations Processing time (s) Seconds per Iteration 576.80 (DTMa) 100 16 29.78 1.86 576.80 150 22 39.85 1.81 576.80 200 18 33.46 1.85 576.80 250 21 38.61 1.83 576.80 300 16 29.82 1.86 246.00 (DTMb) 100 16 12.77 0.79 246.00 150 22 16.84 0.76 246.00 200 18 14.24 0.79 246.00 250 21 16.35 0.77 246.00 300 16 12.83 0.80 The advantage of the iterative model is that it uses three-dimensional models based on measurements such as bathymetry and topography that represent reality at a given spatial resolution. The accuracy of the algorithm results, as well as the raster and flood shape, depend on three factors: accurate data for constructing the DTM, the spatial resolution of the pixel, and the 89 selected error threshold. The complexity of underwater morphometry is shown in Figure 23. Four layers of the water surface are superimposed at every 25 centimeters in height in a stack to distinguish the nonlinearity of morphometric characteristics geographically. Figure 23. Water surface coverage map at different heights above sea level of the Bustillos Lagoon. Sources: Rojas Villalobos with data from Rojas-Villalobos (2018). The value of morphometric variables as the height of the water surface rises above the DTM does not show a constant pattern that can define a precise correlation between them (Figure 24). 90 Figure 24. Comparative graph of volume, surface area, average depth, and maximum depth according to the height above sea level. Source: Rojas Villalobos. The inflection points of the area and the volume in the previous graph, however, are in 1971.5 and 1972.0 m.a.s.l correspondingly. In this way, it is possible to establish equations by segments for each of the parameter, but not a system of equations that integrates the five variables as determined by the algorithm. Recommendations The iterative algorithm proved efficient in finding every one of the morphometric values of the Bustillos Lagoon within the proposed error threshold. The following recommendations, however, are listed. • The purpose of this document is not to evaluate the quality of DTM. To obtain accurate morphometric data and detailed and realistic coverage maps, however, the DTM must meet geographic accuracy (lowest error) in all three axes (X, Y, and Z). 91 • Use reasonable pixel dimensions of the study area. When there are more pixels, the processing time is greater. • This iterative model is not restricted to using a DTM; it is possible to replace with a triangulated irregular network (TIN), which is composed of triangles where the vertices are elevation points. • The algorithm can be implemented in any programming language that handles spatial components, such as Python-GDAL®, R® statistics, or Magik Smallworld®. • Despite the PC's computing capacity, the algorithm can be applied to any computer with minimum requirements: 4 GB RAM, HD with enough space to store the simulations (100-150 GB) and a fast processor (2.0 GHz +). 92 References ArcMap (Version 10.6.1) [English, Windows]. (2019). Retrieved from http://www.esri.com Chen, W., Nover, D., He, B., Yuan, H., Ding, K., Yang, J., & Chen, S. (2018). Analyzing inundation extent in small reservoirs: A combined use of topography, bathymetry and a 3D dam model. Measurement, 118, 202–213. https://doi.org/10.1016/j.measurement.2018.01.042 Erena, M., Atenza, J. F., García-Galiano, S., Domínguez, J. A., & Bernabé, J. M. (2019). Use of Drones for the Topo-Bathymetric Monitoring of the Reservoirs of the Segura River Basin. Water, 11(3), 445. https://doi.org/10.3390/w11030445 GE Corporation. (2019, September 9). MDT Magik Development Tools. Retrieved September 28, 2019, from https://www.mdt.net/ INEGI. (2019). Mapa Digital de México. Retrieved January 16, 2019, from https://www.inegi.org.mx/temas/mapadigital/ Johansson, H., Brolin, A. A., & Håkanson, L. (2007). New Approaches to the Modelling of Lake Basin Morphometry. Environmental Modeling & Assessment, 12(3), 213–228. https://doi.org/10.1007/s10666-006-9069-z Kottek, M., Grieser, J., Beck, C., Rudolf, B., & Rubel, F. (2006). World map of the Köppen- Geiger climate classification updated. Meteorologische Zeitschrift, 15(3), 259–263. https://doi.org/10.1127/0941-2948/2006/0130 Lehman, J. T. (1975). Reconstructing the rate of accumulation of lake sediment: The effect of sediment focusing. Quaternary Research, 5(4), 541–550. https://doi.org/10.1016/0033- 5894(75)90015-0 Moreno-Amich, R., & Garcia-Berthou, E. (1989). A new bathymetric map based on echo- sounding and morphometrical characterization of the Lake of Banyoles (NE-Spain). Hydrobiologia, 185(1), 83–90. https://doi.org/10.1007/BF00006070 Neumann, J. (1959). Maximum Depth and Average Depth of Lakes. Journal of the Fisheries Research Board of Canada, 16(6), 923–927. https://doi.org/10.1139/f59-065 Open Source Geospatial Foundation. (2019). GDAL — GDAL documentation. Retrieved September 20, 2019, from https://gdal.org/ Python Language Reference (Version 2.7.13) [English, Windows]. (2019). Retrieved from http://www.python.org R Foundation. (2019). R: The R Project for Statistical Computing. Retrieved September 28, 2019, from https://www.r-project.org/ Rode, M., Arhonditsis, G., Balin, D., Kebede, T., Krysanova, V., Griensven, A. van, & Zee, S. E. A. T. M. van der. (2010). New challenges in integrated water quality modelling. Hydrological Processes, 24(24), 3447–3461. https://doi.org/10.1002/hyp.7766 Rojas-Villalobos, H. L., Alatorre-Cejudo, L. C., Stringman, B., Samani, Z., & Brown, C. (2018). Topobathymetric 3D model reconstruction of shallow water bodies through remote sensing, GPS, and bathymetry. TECNOCIENCIA Chihuahua, 12(1), 42–54. https://doi.org/10.5281/zenodo.3341734 93 Servicio Meteorológico Nacional. (2019). Información Climatológica por estado. Retrieved January 25, 2019, from https://smn.conagua.gob.mx/es/informacion-climatologica-por- estado?estado=chih Sima, S., & Tajrishy, M. (2013). Using satellite data to extract volume–area–elevation relationships for Urmia Lake, Iran. Journal of Great Lakes Research, 39(1), 90–99. https://doi.org/10.1016/j.jglr.2012.12.013 94 Conclusions This document combines three lines of research that are directly interrelated. The 3D model of the Bustillos lagoon uses known variables such as volume, area, and depth, to estimate the volume of evaporated water according to the evaporation rates obtained by remote sensors. The iteration algorithm uses as a basis the 3D model to compute volume and area that, together with evaporation, indirectly estimate other water balance variables such as water infiltration into the aquifer from the lagoon. At the end, this document integrates useful tools and applicable knowledge in the real world. The databases generated for the region fill gaps of information that is necessary for the analysis of the water balance and the administration of water in the basin. The results obtained will be public for those researchers, government agencies, institutions of higher education, and people interested in these issues. When this dissertation was proposed to the doctoral committee, I was warned of how complex, demanding, and challenging it could be; they were not wrong. The developments and processes that took each of the chapters required knowledge and skills from areas as diverse as electronics, programming, hydrology, physics, mathematics, geography, autonomous aerial vehicles, geographic information systems, and remote sensing. It is crucial to establish that the skills mentioned above, and knowledge are the results of a cumulative learning process along many years of study, an intense desire to acquire knowledge, and a strong curiosity of how hydrologic process occur. 95 References Diario Oficial de la Federación. (2015, July 6). Acuerdo por el que se da a conocer el resultado de los estudios técnicos de aguas nacionales subterráneas del Acuífero Cuauhtémoc, clave 0805, en el Estado de Chihuahua, Región Hidrológico Administrativa Río Bravo. Retrieved from http://www.dof.gob.mx/nota_detalle.php?codigo=5399497&fecha=06/07/2015 Díaz Caravantes, R. E., Bravo Peña, L. C., Alatorre Cejudo, L. C., & Sánchez Flores, E. (2014). Geospatial analysis of land use and water interaction in the peri-urban area of Cuauhtémoc, Chihuahua: A socio-environmental study in northern Mexico. Investigaciones Geográficas, (83), 116–130. García, E. (1964). Modificaciones al sistema de clasificación climática de Koppen (2004th ed., Vol. 8). Mexico: UNAM. Ibañez Hernandez, O. (2010). Manejo del acuífero de Cuauhtémoc, Chih. por el Comité Técnico de Aguas Subterráneas. Brasilia, Brazil: Programa de las Naciones Unidas para el Medio Ambiente. Retrieved from http://bva.colech.edu.mx/xmlui/bitstream/handle/1/1369/ag0198.pdf?sequence=1 INEGI. (2010). Censo de Población y Vivienda 2010. Retrieved November 10, 2015, from INEGI website: http://www.inegi.org.mx/est/contenidos/proyectos/ccpv/cpv2010/Default.aspx INEGI. (2014). Anuario estadístico y geográfico de Chihuahua 2014. Retrieved from http://www.inegi.org.mx/prod_serv/contenidos/espanol/bvinegi/productos/anuario_14/70 2825065409.pdf INEGI. (2016). Continuo de Elevaciones Mexicano 3.0 (CEM 3.0). Retrieved March 16, 2016, from Datos de Relieve website: http://www.inegi.org.mx/geo/contenidos/datosrelieve/continental/continuoelevaciones.as px Ortiz Franco, P. O., & Amado Alvarez, J. P. (2001). Uso del Agua de la Laguna de Bustillos para la Producción de Maíz. Terra Lationamericana, 19(2), 183–189. Quintana S., V. M. (2014). Nuevo orden alimentario y disputa por el agua en el norte de México. Apuntes: Revista de Ciencias Sociales, 40(73), 175–202.