A remote sensing method for estimating regional reservoir area and evaporative loss

Journal of Hydrology 555 (2017) 213–227

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

A remote sensing method for estimating regional reservoir area and evaporative loss Hua Zhang a,⇑, Steven M. Gorelick b, Paul V. Zimba c, Xiaodong Zhang d a

School of Engineering and Computing Sciences, Texas A&M University-Corpus Christi, Corpus Christi, TX 78412, United States Department of Earth System Science, Stanford University, Stanford, CA 94305, United States c Center for Coastal Studies, Texas A&M University-Corpus Christi, Corpus Christi, TX 78412, United States d Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NW 87545, United States b

a r t i c l e

i n f o

Article history: Received 2 June 2017 Received in revised form 10 September 2017 Accepted 4 October 2017 Available online 7 October 2017 This manuscript was handled by Emmanouil Anagnostou, Editor-in-Chief, with the assistance of Xinyi Shen, Associate Editor Keywords: Evaporation Reservoir Google Earth Engine MNDWI Freshwater

a b s t r a c t Evaporation from the water surface of a reservoir can significantly affect its function of ensuring the availability and temporal stability of water supply. Current estimations of reservoir evaporative loss are dependent on water area derived from a reservoir storage-area curve. Such curves are unavailable if the reservoir is located in a data-sparse region or questionable if long-term sedimentation has changed the original elevation-area relationship. We propose a remote sensing framework to estimate reservoir evaporative loss at the regional scale. This framework uses a multispectral water index to extract reservoir area from Landsat imagery and estimate monthly evaporation volume based on pan-derived evaporative rates. The optimal index threshold is determined based on local observations and extended to unobserved locations and periods. Built on the cloud computing capacity of the Google Earth Engine, this framework can efficiently analyze satellite images at large spatiotemporal scales, where such analysis is infeasible with a single computer. Our study involves 200 major reservoirs in Texas, captured in 17,811 Landsat images over a 32-year period. The results show that these reservoirs contribute to an annual evaporative loss of 8.0 billion cubic meters, equivalent to 20% of their total active storage or 53% of total annual water use in Texas. At five coastal basins, reservoir evaporative losses exceed the minimum freshwater inflows required to sustain ecosystem health and fishery productivity of the receiving estuaries. Reservoir evaporative loss can be significant enough to counterbalance the positive effects of impounding water and to offset the contribution of water conservation and reuse practices. Our results also reveal the spatially variable performance of the multispectral water index and indicate the limitation of using scene-level cloud cover to screen satellite images. This study demonstrates the advantage of combining satellite remote sensing and cloud computing to support regional water resources assessment. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Reservoirs have been constructed for more than 4000 years primarily for drinking water and irrigation needs, as well as for modern functions of hydropower, flood control and navigation (Assouline et al., 2011). Evaporation from the water surface of a reservoir can significantly affect its function in ensuring the availability and temporal stability of water supply, exacerbating the problem of water scarcity. For example, the evaporation from agricultural irrigation reservoirs in a semiarid basin in Spain is equivalent to 27% of the total domestic water use in this basin (Martínez Alvarez et al., 2008). It is estimated that half of water stored in Australia’s reservoirs is lost due to evaporation (Rost et al., 2008). ⇑ Corresponding author. E-mail address: [email protected] (H. Zhang). https://doi.org/10.1016/j.jhydrol.2017.10.007 0022-1694/Ó 2017 Elsevier B.V. All rights reserved.

Quantifying reservoir evaporation requires reliable measurements of evaporation rate and reservoir area. Regarding evaporation rate, a variety of methods are available, ranging from local-scale approaches such as eddy covariance tower and scintillometer to regional-scale energy balance methods that are based on water surface temperature (Duan and Bastiaanssen, 2013; McGloin et al., 2014; Phillips et al., 2016). Nevertheless, the use of evaporation pans remains the most popular method for water resource managers, despite its sensitivity to local weather and landscape conditions (Lowe et al., 2009; Phillips et al., 2016). The satellite-derived MOD16 dataset has been widely used in studies on evapotranspiration, but it does not cover rivers, lakes, reservoirs, and wetlands (Mu et al., 2013; Mu et al., 2016). It is designed to describe evapotranspiration from vegetated land surfaces and its input includes Leaf Area Index, which is not applicable for water bodies.

214

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

Regarding reservoir area, the reservoir area-elevation curve is a widely-accepted approach to estimate the fluctuation of reservoir area (Gao, 2015). However, reservoir curves could be unavailable if the reservoir is in a data-sparse region, or questionable if longterm sedimentation has changed the original elevation-area relationship. The error in reservoir area estimates is largely neglected in analyzing the uncertainties of reservoir evaporation (Lowe et al., 2009). Remote sensing provides a direct means to observe variations of reservoir area. Established methods include various spectral indices based on visible and infrared bands of satellite sensors, as well as their combinations with image enhancing and classification approaches (Gao et al., 2012). The normalized difference water index (NDWI) and its variants (Feyisa et al., 2014; Ji et al., 2009; Xu, 2006) have been suggested as effective approaches as they are insensitive to subpixel vegetation component. The performance of these spectral water index methods depends on the determination of a threshold that differentiates water pixels from non-water pixels. This threshold is usually manually calibrated and its value reflects the combined effects of atmospheric absorption, lake water quality, shoreline land cover, and the spatial resolution of the satellite sensor (Duan and Bastiaanssen, 2013; Ji et al., 2009; Rokni et al., 2014; Yang et al., 2015). However, in prior research the spatial pattern of the water index threshold has been analyzed only based on a small group of water bodies or using coarse-resolution images (Du et al., 2016; Li et al., 2015; Li et al., 2013). Separating water from other surfaces in one satellite image is not a new problem, but consistent discrimination on the regional or global scale over multiple decades is a non-trivial challenge (Pekel et al., 2016). This would involve the processing of tens of thousands of satellite images. Such effort is too computationally expensive to be implemented using conventional approaches that rely on a couple of computers. Cloud computing platforms, such as Google Earth Engine (GEE), Amazon’s NASA NEX and Los Alamos National Laboratory’s Descartes Labs, have emerged as a highlyefficient technology for big data exploration and information extraction (Warren et al., 2015). They are particularly valuable in supporting large-scale spatiotemporal analysis with remote sensing, ranging from global forest cover change to transboundary freshwater crisis (Burke and Lobell, 2017; Hansen et al., 2013; Müller et al., 2016; Pekel et al., 2016). However, no study has been undertaken to employ such platforms to analyze the variations of reservoir area and evaporation. Therefore, the goal of this study is to develop a remote sensing framework for analyzing long-term reservoir evaporative loss at the regional scale. The framework combines GEE and local observations, addresses the spatial heterogeneity of water index thresholds, and evaluates the significance of reservoir evaporative loss to water resources management. Texas is well known for its diverse hydroclimatic conditions and its long history of reservoir development that has resulted in a network of over 3,400 reservoirs. This network provides a comprehensive context to evaluate the performance of satellite-based methods for reservoir evaporation estimation. The proposed framework is demonstrated through a case study of Texas’s 200 major reservoirs captured in 17,811 Landsat images over a 32-year period.

2. Material and methods 2.1. Study area Climate, geography, hydrology, and water resources management vary greatly across Texas, a 685,000-km2 area representative of both the drier western and wetter eastern regions of the United States (Wurbs and Ayala, 2014). In 2015, the total water

use in Texas was 15.2 billion m3, mainly distributed among irrigation (50%), municipal (34%), manufacturing (8%), power (3%), livestock (3%), and mining (1%) water use categories; over 40% of water withdrawal was from surface water sources (TWDB, 2017b). Between 2020 and 2070, Texas’s population is expected to increase 73% from 29.5 million to 51.0 million, while water supplies are expected to have a 10% decline (TWDB, 2017a). Growing population and declining groundwater reserves are resulting in intensifying demands on surface water resources (Wurbs and Ayala, 2014). The core component of Texas’ surface water resources management is a water right permit system that consists of over 3,400 reservoirs with conservation capacity greater than 2.5 million m3 (MCM). Conservation capacity is the volume of the reservoir’s conservation pool, which lies above the top of the dead pool and below the bottom of the flood pool. It is also referred as the active storage volume or capacity. Nearly all reservoirs in Texas are artificial lakes and most of them were built between 1960 and 1990 (Wurbs and Ayala, 2014). Our study only focuses on 200 major reservoirs with conservation capacity larger than 5.2 MCM (Fig. 1). All selected reservoirs are visible on Landsat imagery with cloud cover below 10%. Thus there is no uncertainty about the number of reservoirs and their locations. Their total conservation capacity is 40.9 billion m3 and accounts for 79% of Texas’s total capacity. The Government of Texas has continuous water area estimates only for one-third of these reservoirs. Monthly evaporation and precipitation rates of reservoirs are obtained from the Texas Water Development Board (TWDB) in one-degree latitude by one-degree longitude quadrangles (Fig. 2). The dataset is based on daily precipitation rates measured at 2,400 precipitation stations and daily evaporation rates at 76 pan evaporation stations across Texas (TWDB, 2016). The daily gross pan evaporation rates are first summed into monthly pan evaporation rates at each evaporation station. Then Thiessen polygons are developed based on the locations of all evaporation stations. The quadrangular pan evaporation rate is calculated as the sum of pan evaporation rates of all stations weighted by the areas of their Thiessen polygons in the quadrangle. It is further converted into a quadrangular reservoir evaporation rate by multiplying by pan-tolake coefficients, which are developed by the National Weather Service at the monthly level and are specific for each quadrangle. The quadrangular precipitation dataset is produced in a similar way except that there is no need for pan-to-lake precipitation coefficients. 2.2. Satellite data Satellite data used in this study are surface reflectance derived from Landsat 5 and Landsat 8 imagery. Landsat 5, the longest-operating Earth observation satellite, was launched in March 1984 and retired in January 2013. Landsat 5 TM images have six 30-m visible bands and a 120-m thermal band. Landsat 8 was launched in February 2013 as the newest member of the Landsat family. Landsat 8 OLI images have eight 30-m bands, two 100-m thermal bands and a 15-m panchromatic band. This study used a total of 17,811 Landsat images across 47 WRS-2 scenes (Fig. 2), including 15,748 Landsat 5 TM images from 1984 to 2011, and 2,063 Landsat 8 OLI images from 2013 to 2015. Due to the SLR-off malfunction and the lack of reliable algorithms to fix this problem, Landsat 7 ETM+ images were not included in this study, which resulted in a data gap in 2012. The surface reflectance data are produced by the U.S. Geological Survey (USGS) using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al., 2006) and accessed through GEE.

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

215

Fig. 1. Study area. Blue circles show 200 reservoirs in 20 basins in Texas. The size of the circle indicates the conservation capacity (active storage volume) of the reservoir. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2.3. Methods The proposed framework is shown in Fig. 3. It integrates GEE and local observations to analyze monthly evaporation volume of multiple reservoirs and to evaluate the significance of evaporative loss to regional water resources management. 2.3.1. Image preparation Initially, the designed boundary of each reservoir at its conservation capacity is extracted from the GIS data provided by TWDB. The boundary is then expanded by buffering the boundary polygon by 1,000 m. This adjusted boundary is sufficient to fully contain the variation of reservoir area over the study period. Confining image analysis to this adjusted boundary rather than the entire satellite image scene improves computation efficiency and reduces the presence of other landscape features that may affect water identification. If the adjusted reservoir boundary is fully covered in a single Landsat WRS-2 scene, images of this scene are used individually for water identification. Area estimates from multiple images within a month will be averaged to generate a monthly estimate. If the reservoir boundary spans two Landsat WRS-2 scenes, the involved images are first mosaicked into monthly images before water extraction. The mosaicking process deals with overlapping

area using a cloud mask to maximize the fraction of clear area in the mosaic image. The cloud mask is generated using the Fmask algorithm (Zhu et al., 2015) and contains pixel values of 0, 1, 2, 3, and 4 for clear, water, shadow, snow, and cloud conditions, respectively. For each selected single-scene image or mosaicked multiscene image, a cloud score is calculated as the ratio of the number of clear and water pixels to the total number of pixels. Pixels covered by cloud, shadow and snow conditions are excluded from further analyses. Images with cloud score higher than a specified threshold (P = 0.9) are selected for water identification. As this cloud score is generated at the scale of a single feature (reservoir), it is different from the cloud score provided by USGS, which is calculated over the entire scene. The comparison of feature-level and scene-level cloud scores in term of their effects on data availability is discussed in Section 3.4. 2.3.2. Water extraction and validation Water is identified in pre-processed satellite images using the modified normalized difference water index (MNDWI) (Xu, 2006):

MNDWI ¼

qGreen  qSWIR qGreen þ qSWIR

ð1Þ

216

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

Fig. 2. Ground observation data and satellite imagery of the study area: (a) quadrangle precipitation rate, (b) quadrangle evaporate rate, (c) total number of Landsat images (1984–2015), and (d) average cloud cover (1984–2015). Black dots indicate reservoirs. Colored polygons in (c) and (d) indicate Landsat WRS-2 scenes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where qGreen and qSWIR denote the reflectance of green and shortwave infrared (SWIR) bands, respectively. In Landsat 5 images, qGreen and qSWIR are bands 2 and 5, respectively; in Landsat 8 images, they are bands 3 and 6, respectively. This index is proven to be less sensitive to subpixel impurity than other spectral water indices (Duan and Bastiaanssen, 2013; Ji et al., 2009). The value of MNDWI varies from 1 to 1, with larger values indicating a higher fraction of water in the pixel. Water pixels have positive MNDWI values because of their higher reflectance in the green band than in the SWIR band, while non-water pixels (soil and vegetation) have negative MNDWI values due to their lower reflectance in the green band than in the SWIR band (Xu, 2006). Based on MNDWI values, water pixels are differentiated from non-water pixels by applying a threshold. In contrast to the constant threshold of zero for NDWI, the use of variable thresholds has been recommended for MNDWI. The adjustment of threshold to reflect subpixel land-cover components can be achieved using ground truth or a spectral library (Ji et al., 2009). This study assumes that the optimal threshold of MNDWI is spatially heterogeneous but temporally stable. Local observations of reservoir area are available at some reservoirs, generated using daily measurements of reservoir elevation and the reservoir elevation-area curve. For these reservoirs where observations are available, the satellite-derived monthly area estimates are

compared to local observations using three performance indices: the coefficient of determination (R2), the Nash–Sutcliffe model efficiency coefficient (NSE) and the Percent Bias (PBIAS). The R2 measures the percent of variation of the observed area explained by the satellite-estimated area, the NSE describes to what extent the remote sensing approach is better than using the mean of the observed area, and the PBIAS measures the deviation of the satellite-estimated area from the observed area. The comparison is made over a range of MNDWI thresholds from 0.5 to 0.1 at the interval of 0.05. Then values of each performance index are sorted by the degree of fitness from the best to the worst. An overall performance score is calculated at each threshold value as the weighted sum of rankings of three performance indices. The threshold with the lowest overall performance score is determined as the optimal threshold. This calibration procedure is implemented for each reservoir that has area observation records. The mean of optimal thresholds for calibrated reservoirs is then applied to all other reservoirs that do not have water area records. Efforts of image preparation and analysis are fully conducted through GEE. It is a computational infrastructure optimized for parallel processing of geospatial data, coupled with a dedicated programming interface and online access to a variety of global satellite imagery products including the Landsat archive (Pekel

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

217

Fig. 3. Research framework. The abbreviations LT5-SR, LC8-SR and MNDWI stand for Landsat 5 TM surface reflectance, Landsat 8 OLI surface reflectance and modified normalized difference water index, respectively.

et al., 2016). As a revolutionary geospatial analysis platform, GEE incorporates cloud computing to replace conventional offline efforts and enables efficient analysis of large-scale spatiotemporal analysis. The algorithms of this study are implemented through the JavaScript API of GEE. All satellite images are accessed, analyzed and managed online without any substantial requirement of local computing and data storage resources. 2.3.3. Evaporation estimation Monthly evaporation volume is calculated as the product of satellite-estimated reservoir area and observed reservoir evaporation rate. Monthly net evaporation volume is the product of satellite-estimated reservoir area and net evaporation rate (i.e.,

monthly reservoir evaporation rate minus reservoir precipitation rate). The net evaporation rate takes into the consideration precipitation to quantify water loss and gain through reservoir surface. For reservoirs with multiple area estimates within a month, the average area is used. Evaporation and precipitation rates are based on the 1°  1° quadrangle TWDB data. Each reservoir is assigned to the closest quadrangle based on the coordinates of the reservoir’s geometric center. Statistics of monthly evaporation and net evaporation volumes over 32 years are analyzed at both reservoir and basin levels. The estimates of evaporation volume are compared to reservoir conservation capacity and estuarine freshwater inflow to evaluate the significance of reservoir evaporative loss to water resources management.

218

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

3. Results and discussion 3.1. Estimates of reservoir area Fig. 4 shows satellite-derived monthly areas of six reservoirs from 1984 to 2015, compared to area observations derived from their elevation-area curves. Overall, satellite results agreed well

with the observations. The area variations shown in satellite results reflect both seasonal fluctuations of reservoir water storage in normal years and significant draw-downs caused by major droughts such as the record drought in 2010–2011. More importantly, the remote sensing methods extends the reservoir area estimates to the early 1980s when observations were not available for many reservoirs. Most reservoirs in the study area were built in the

Fig. 4. Satellite-based area estimates for six reservoirs. Column 1 compare satellite estimates (red circles) and elevation-based estimates (black lines), followed by an example of water extraction for each reservoir, including the satellite image (column 2), MNDWI (column 2) and water pixels (column 4). The reservoir locations are shown in Fig. 1 (white circles). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

219

Fig. 5. Calibration of MNDWI thresholds for 73 reservoirs: (a) optimal MNDWI thresholds, (b) coefficient of determination (R2), (c) Nash-Sutcliffe efficiency (NSE), and (d) percent bias (PBIAS).

1960s or 1970s and only one-third of them have consistent area observations. Even if the area observations are available, their records are typically less than 15 years long. Satellite estimates fill the information gap and provide a direct, continuous observation of 32 years for multiple reservoirs. Such information is important for analyzing the regional long-term patterns of reservoir area and storage changes. Fig. 5 shows the optimal MNDWI thresholds of 73 reservoirs and associated values of performance indices. The optimal thresholds were 0.17 ± 0.12, resulting in R2 of 0.78 ± 0.21, PBIAS of 3.97 ± 3.16%, and NSE of 0.66 ± 0.27, demonstrating an overall good performance of satellite estimates. The weights of performance indices were set as 0.5, 0.3 and 0.2, respectively, based on the best fit to the data of all 73 reservoirs available for calibration. As these indices agreed quite well in indicating the optimal thresholds, the selection of their weights had no substantial impact on the results of this study. A comparison of optimal thresholds identified under different weights is presented in the Supplementary Information (Fig. S1). The MNDWI threshold is a function of subpixel fractions of soil and vegetation and generally increases with vegetation fraction value (Ji et al., 2009). Positive MNDWI thresholds were identified for 14% of the reservoirs, mostly in the upstream areas of the Colorado, Brazos and Trinity Basins, where high vegetation cover is

present in the shoreline areas of these reservoirs. The average optimal threshold (0.17) was applied to 127 reservoirs that had no area observations. Based on visual inspection using 1m NAIP aerial imagery, this threshold yielded better results than using the optimal threshold from the nearest reservoir or the default threshold of zero. Fig. 6 shows the performance of MNDWI over a range of thresholds as measured by three performance indices. In this study, the responses of NSE and R2 are largely consistent. The NSE and R2 are widely used in hydrological studies, but they are sensitive to large values, at the expense of better performance for small values. The response of PBIAS is more intuitive, as a larger threshold leads to a larger estimate of water surface and thus more likely overestimation, and vice versa. The PBIAS may indicate poor model performance (Moriasi et al., 2007) in situations of systematic over- or underestimation where NSE and R2 may still have optimistic values. Using PBIAS = 0 alone to determine the optimal threshold might be problematic, as the accumulation of the error term does not reflect its variability at different time steps. In this study, the three performance indices agreed generally well in indicating the optimal thresholds. The determined optimal thresholds reflects a combined consideration on the variability and average magnitude of error (Krause et al., 2005), and they

220

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

Fig. 6. Determination of optimal MNDWI thresholds for 16 reservoirs based on NSE (solid black line), R2 (dashed black line) and PBIAS (dashed blue line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

are not related to reservoir conservation capacity or the length of observation records (Fig. S2). Our determination of MNDWI thresholds is different from previous studies. First, the threshold is evaluated using multi-temporal images over a long term, in contrast to efforts based on a single image. Second, the use of multiple performance indices results in a comprehensive evaluation of the fitness between satellite estimates and local observations. Third, the use of Google Earth Engine allows for the identification of a regional pattern of MNDWI thresholds, which would be difficult to achieve in conventional offline image analysis. 3.2. Estimates of reservoir evaporation Fig. 7 shows average monthly evaporation and net evaporation (accounting for precipitation) for 20 basins. Despite the differences in geographic and climatic conditions of these basins, the results have three patterns in common. First, net evaporation is more variable than evaporation, as precipitation can be quite variable. Second, monthly evaporation shows a unimodal pattern with the peak in June or July, and monthly net evaporation appears to have a bimodal pattern with a local low in April (Fig. 7e and l). This is mainly caused by heavy springtime precipitation in Texas (Longley, 1994). Third, the three-month period from June to August has an absolute loss of freshwater, accounting for 38% of annual total evaporation and 67% of annual total net evaporation; net evaporation is negative from November to February for nine basins, indicating this is an important period of freshwater gain. Across Texas, the total annual evaporation from 200 reservoirs is estimated at 8.0 billion m3 (Fig. S3). It equals to 20% of total conservation capacity (reservoir active storage) and represents a significant portion of the annual water

budget. Compared to state water use in 2015 of 15.0 billion m3 (TWDB, 2017b), reservoir evaporative loss is equivalent to 53% of total water use, 106% of irrigation water use and 156% of municipal water use. Table 1 shows the annual average reservoir area, total evaporation and total net evaporation at the basin level. Basins in the eastern region have stronger evaporation than their western counterparts. Particularly, the reservoirs in Trinity, Sabine and Brazos Basins account for half of the total evaporation from 200 reservoirs. As indicated by net evaporation, evaporative loss is compensated by rainfall on the reservoir surface. This effect is equivalent to 35% of evaporation on average over the 1984–2015 period, but it varies significantly among basins, ranging from over 95% for the Cypress, Sulphur and San Jacinto Basins in the humid region to only 21% for the Rio Grande Basin located in the semiarid region. Rigorous verification of reservoir evaporative loss remains a difficult task in reservoir management. There are no ground measurements available to directly validate the satellite-derived reservoir evaporation estimates. The use of scintillometer, eddy covariance tower, or a combination of both could be a solution to the quantification of energy fluxes from open water surface. However, pilot studies have been conducted only for small water bodies (McGloin et al., 2014; McJannet et al., 2013). There is little evidence on the performance of such technologies over major reservoirs where heat storage, microclimate, water quality, and shoreline landscape are more complex and heterogeneous than small water bodies. We compared the satellite-based estimates to results from a previous study of watershed hydrological modeling. Wurbs and Ayala (2014) conducted the first state-wide study on reservoir evaporation using the Texas Water Availability Modeling (WAM) System. They calculated reservoir area

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

221

Fig. 7. Average monthly evaporation and net evaporation for 20 basins. Blue and green lines are monthly evaporation and net evaporation volumes, respectively. Error bars cover one standard deviation. Red lines indicate zero net evaporation, e.g., the magnitude of monthly evaporative loss is the same as monthly precipitation. The locations of basins are shown in Fig. 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

based on simulated reservoir storage and storage-area curves and applied the TWDB 1°  1° evaporation rates to estimate evaporation volume. Their modeling results are consistent with satellite results of this study at the basin level (Fig. S4), but there are relative overestimations using simulation in the Brazos Basin and relative overestimations using simulation in the Colorado and Trinity Basins. Such discrepancies are associated with the methods of area estimation. First, as in most hydrological simulations, the WAM model was calibrated by comparing simulated flows to observed flows at stations. The simulated reservoir storage and storage-derived reservoir area were not verified. Reservoir storage was simulated using a monthly reservoir water balance and subjected to uncertainties in the model estimates of reservoir seepage and evaporation. Second, the satellitebased monthly estimate was the average of two snapshots that were captured 16 days apart. It didn’t account for the daily

variations between the two dates and was fundamentally different from the monthly water balance used in the hydrological model. Third, the model relied on a static storage-area curve to estimate reservoir area. This curve might need to be updated if sedimentation has significantly changed the elevation-storagearea relationship. 3.3. The significance of reservoir evaporation to regional water resources management The significance of reservoir evaporation is evaluated from three perspectives. First, the average annual reservoir evaporation is compared to the reservoir’s conservation capacity (Fig. 8). When a reservoir’s conservation pool is empty, the reservoir is considered inoperable for water supply. A high ratio of annual evaporation to conservation capacity would indicate that evaporative loss is

222

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

Table 1 Basin annual average reservoir area, total evaporation and total net evaporation. No

Basin Name

Number of Reservoirs

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Brazos Canadian Colorado Cypress Guadalupe-San Antonio Lavaca Neches Neches-Trinity Nueces Nueces-Rio Grande Red Rio Grande Sabine San Jacinto Colorado-Lavaca Sulphur Trinity Trinity-San Jacinto San Jacinto-Brazos San Antonio-Nueces

44 3 31 10 8 1 11 1 3 4 25 8 12 6 2 4 32 1 2 2

Area (km2)

Evaporation (106 m3)

Net Evaporation (106 m3)

mean

sd.

mean

sd.

mean

sd.

7286.1 351.9 4691.2 2118.3 952.4 430.9 6920.2 583.8 1709.4 187.1 6549.4 4814.3 10992.5 1622.8 88.7 1966.7 14938.9 197.2 90.1 150.4

89.0 33.4 111.4 21.5 23.4 8.6 101.0 27.8 86.1 3.0 104.5 318.3 113.1 21.1 3.5 104.4 223.2 2.3 1.7 4.7

915.9 53.2 591.0 224.1 109.0 45.6 744.3 60.3 202.5 25.5 808.2 714.0 1278.3 176.4 9.5 194.6 1792.5 20.5 9.4 16.5

14.8 6.1 17.6 5.5 3.6 2.4 28.7 3.7 13.0 1.0 26.4 59.9 43.9 5.1 0.5 14.3 40.5 0.9 0.3 0.8

412.9 38.3 318.9 10.5 52.6 12.1 101.7 0.6 106.6 15.4 318.6 566.3 214.9 30.0 2.8 3.8 619.2 0.7 0.4 6.8

33.0 6.1 27.4 23.8 8.9 10.2 84.3 12.6 20.9 2.1 77.3 66.8 147.8 19.6 1.4 28.0 91.0 4.3 1.4 2.0

Fig. 8. Comparison of reservoir evaporation to conservation capacity: (a) the ratio of average annual evaporation to conservation capacity, (b) the ratio of average annual net evaporation to conservation capacity, and (c) the difference between (a) and (b).

significant to the water supply function of this reservoir. This ratio was 0.29 on average for 193 reservoirs; seven reservoirs were excluded due to the lack of information of conservation capacity. It exceeds 0.4 for 39 reservoirs (20%), varies between 0.2 and 0.4 for 83 reservoirs (43%), and is below 0.2 for 71 reservoirs (37%). Reservoirs with high ratios are located in the Red, Trinity and Colorado Basins, and coastal regions surrounding the Trinity-San Jacinto Estuary. Replacing evaporation by net evaporation (Fig. 8b), this ratio was 0.10 on average and above 0.4 only for four reservoirs. Second, the average monthly evaporation is compared to the minimum monthly freshwater inflows mandated by the Government of Texas for five major estuaries. The minimum freshwater inflows are derived from an optimization system that consists of a set of multiple linear regression models to characterize the relationships between freshwater inflows, salinity and fishery productivity (Powell et al., 2002). The resulted regime of monthly freshwater inflow is deemed as the minimum environmental flow required to sustain ecosystem health and fishery productivity of the receiving estuary. This regulation has been applied to all major estuaries in Texas to establish goals of freshwater release. This

study used the minimum freshwater inflows of five estuaries: Sabine-Neches, Trinity-San Jacinto, Lavaca-Colorado, Guadalupe, and Nueces (Fig. 1). As shown in Fig. 9, annual reservoir evaporative loss is on average equivalent to 10–47% of the minimum freshwater inflow, but it is more significant in July and August when the evaporation may exceed the minimum freshwater inflow. In the Nueces Estuary (Fig. 9e), the evaporative loss exceeds the environmental flow in every month except in May, and for half of the year, the evaporative loss is at least 150% higher than the minimum environmental flow. The Nueces Basin has made heavy investments in developing alternative water sources, building infrastructures for water diversion and promoting urban water conservation, in order to meet the goal of ensuring the minimum freshwater inflow to the Nueces Estuary. The results of this study indicate that effective suppression of reservoir evaporation, using chemical or physical measures (Assouline et al., 2011), could provide an alternative solution to alleviate the ecological water crisis. Third, reservoir evaporative loss is compared to reservoir storage change per month from 1984 to 2015. Fig. 10 shows the comparison at six reservoirs. The selected reservoirs have long storage records and are representative of a wide range of reservoir size and

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

223

Fig. 9. Comparison of reservoir evaporation and estuarine minimum environmental inflow.

hydrological conditions. For all selected reservoirs except for the smallest reservoir, No. 107, the storage change is considerably more variable than the evaporative loss. Evaporative loss appears to be a significant portion of storage change for small reservoirs (Fig. 10c and g), compared to large reservoirs (Fig. 10d and e). The evaporative loss even exceeded the storage decrease in several months, for example, the period of 2001–2005 at reservoir No. 37 (Fig. 10c) and the period of 1995–1996 at reservoir No. 54 (Fig. 10d). This suggests that the evaporative loss played a decisive role in regulating reservoir storage in those periods. 3.4. The effect of cloud cover threshold on image availability Cloud cover is an important factor for screening satellite images before image analysis. The percentage of cloud cover over the entire image is often used for this purpose. This scene-level cloud cover is calculated using the Automatic Cloud Cover Assessment

(ACCA) algorithm, which runs multiple spectral tests based on the premise that clouds are brighter and colder than most of the Earth’s surface (Irish, 1998). The spatial pattern of this parameter over the study area is shown in Fig. 2d. As this information is available as part of the metadata of Landsat images, it is widely used by researchers to evaluate cloud cover condition. However, this scenelevel cloud cover information can be misleading for studies on separated local features such as reservoirs. A low scene-level cloud cover doesn’t ensure a clear sky over a targeted small feature, and a high scene-level cloud cover doesn’t rule out the possibility that the targeted reservoir is actually cloud free. Screening images using scene-scale cloud cover may exclude some images that have clear or an acceptable cloud cover condition over the targeted feature, resulting a loss of valuable information. Thus it is more effective to evaluate cloud cover condition at the feature (reservoir) scale, e.g., within the adjusted boundary of the targeted reservoir in this work.

224

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

Fig. 10. Comparison of satellite-derived reservoir evaporative loss (red bars) to observed reservoir storage change (blue bars) per month from 1984 to 2015 at eight reservoirs. Evaporative loss is shown as negative values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Data availability under different thresholds of feature-level cloud cover: (a) the total number of records in each year; and (b) the number of records on average in each month. A record is defined as a clipped portion of an image that covers a targeted reservoir with the percentage of cloud-free area higher than a specified threshold (P). As this study only used images from Landsat 5 and 8, there are no images available from 2012.

Fig. 11 shows the annual and monthly variations of featurelevel cloud cover during the study period. In this study, a record is defined as a clipped portion of an image that covers the targeted reservoir with the fraction of cloud-free area higher than a specified threshold (P). Each record can be analyzed using MNDWI to generate an estimate of reservoir area. As shown in Fig. 11a, the number of records varies significantly over years and is particularly low in 1984–1986, 1990, 2005, and 2011. With

a cloud cover threshold of P = 1.0 (e.g., a completely clear sky over the reservoir), only 190 records on average are available annually. Considering the number of reservoirs (n = 200), this would be equivalent to one area estimate per reservoir per year. Lowering the threshold to P = 0.9, the number of records per year increases significantly to 1,222. This would be equivalent to about six area estimates per reservoir per year, which is helpful for analyzing monthly variations. Fig. 11b shows the monthly variations of

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

225

Fig. 12. Misclassified water area in the Choke Canyon Reservoir (No. 22, 29°300 N, 98°210 W). The solid white line indicates the reservoir boundary at conversation capacity, and the dashed white line indicates the 1-km buffered reservoir boundary. Misclassified pixels are shown in red. Subplots (a) and (b) shows the image and result of October 17, 2001, an example of low water level; misclassified water pixels account for 1.04% of the total identified water area. Subplots (c) and (d) shows the image and result of January 13, 2005, an example of high water level; misclassified water pixels account for 0.35% of the total identified water area. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

feature-level cloud cover. In general, there are more records in spring and fall than in winter and summer. For example, with a threshold of 0.9, the average number of records in August (117.1) is 40% higher than in February (83.6). In contrast, using scene-level cloud cover to screen images will reduce the number of reservoirs covered in the images, resulting in a less efficient data utilization (Fig. S5). Also, at the regional level the rainy reason is not necessarily associated with image availability (Fig. S6), as the periods of precipitation events may not always contain the acquisition dates of Landsat images. As cloudy pixels are excluded from the implementation of MNDWI, records with more cloud-covered pixels would result in underestimation of reservoir area. A low threshold of feature-level cloud cover increases the number of records, but at the same time reduces the quality of records. On the other hand, a high threshold results in too few records available for analysis and even no records for some months. In this study, the threshold of 0.9 was selected for the feature-level cloud cover, as this value results in a reasonable balance between record quantity and quality. 3.5. Uncertainties of reservoir evaporation estimates We identified water pixels within an adjusted reservoir boundary. This boundary was established through buffering the polygon of surveyed reservoir conservation capacity by 1000 m. This uniform width of the buffer zone was set carefully to contain the historically-observed maximum water area of all reservoirs in this study. Extracted water pixels within this adjusted boundary inevitably included non-reservoir waterbodies, leading to overestimated reservoir area. However, the impact of these misclassified water pixels is insignificant for two reasons. First, the majority of misclassified water pixels are caused by small ponds and lakes in

the 1-km buffer zone. Second, most rivers are unidentifiable in Landsat imagery that has a spatial resolution of 30 m. Some river sections might be wide enough to be identified as water pixels. If a river section is located within the polygon of conversation capacity (without the buffer zone), it is indeed a part of the reservoir. If it is located within the buffer zone, its length is limited by the width of the buffer zone. In support of our procedure, we estimated the water misclassification error. On average, the misclassified water pixels of small ponds and river sections account for 1% of the total identified water area for the reservoirs in this study. Two examples are presented in Figs. 12 and 13 to show misclassified water pixels in a large reservoir (the Choke Canyon Reservoir) and a small reservoir (the Greenbelt Lake) at both high and low water levels. In the Choke Canyon Reservoir (Fig. 12), misclassified water pixels account for 1.04% and 0.35% on October 17, 2001 (low water level) and January 13, 2005 (high water level), respectively. In the Greenbelt Lake (Fig. 13), misclassified water pixels account for 0.98% and 1.13% on November 27, 2014 (low water level) and September 17, 1997 (high water level), respectively. Uncertainty is also caused by the use of pan evaporation data to estimate reservoir evaporation rates. First, the actual evaporation from a reservoir may differ significantly from pan measurements because of the heat storage in the reservoir (Farnsworth et al., 1982). During the spring and summer when heat is stored in the reservoir, the pan-derived evaporation rate could be much higher than the lake evaporation rate. During the fall and winter when the stored heat is released, the lake could evaporate much faster than the pan. This disparity cannot be fully addressed through applying monthly pan-to-lake coefficients and is recognized as a major concern when using pan observations to estimate reservoir evaporation. A similar challenge also confronts existing satellite-based evaporation products that rely on visible and infrared bands.

226

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

Fig. 13. Misclassified water area in the Greenbelt Lake (No. 37, 35°000 N, 100°550 W). The solid white line indicates the reservoir boundary at conversation capacity, and the dashed white line indicates the 1-km buffered reservoir boundary. Misclassified pixels are shown in red. Subplots (a) and (b) shows the image and result of November 27, 2014, an example of low water level; misclassified water pixels account for 0.98% of the total identified water area. Subplots (c) and (d) shows the image and result of March 17, 1997, an example of high water level; misclassified water pixels account for 1.13% of the total identified water area. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4. Conclusions This study demonstrated the advantage of a GEE-based water resources assessment formwork that can process geospatial big data to assess regional reservoir evaporative loss. This study also confirmed MDNWI as an effective spectral index for water identification and revealed the spatial heterogeneity of MNDWI thresholds at the regional scale. The scale of geospatial analysis in this study is rare for water research community, due to the tremendous time and resources requirements for obtaining, analyzing and storing the involved satellite images and results. Running the geospatial analysis of this study on a single standard computer would have taken two years, but using over 10,000 computers through Google Earth Engine, the analysis was completed in less than 4 h once the script was developed. It is important to note that our method is not developed to estimate evaporation in an area where all reservoirs are totally ungauged. However, this method is capable of generating a satellite-based direct quantification of reservoir area dynamics, completing or extending the current area estimates, expanding the estimates to ungauged reservoirs that are close to gauged reservoirs. Most importantly, the approach provides an efficient means to analyze evaporation rates and patterns over large spatiotemporal scales, where such analysis is infeasible with a single

computer. The combination of cloud computing and satellite remote sensing can reduce uncertainty and significantly improve the efficiency of regional hydrological analysis. The required precipitation and pan evaporation are among the most common measurements at many weather stations. Landsat imagery archive has global coverage, and Google Earth Engine is accessible to the public at no cost. Therefore, the method developed in this study can be readily applied to many regions around the world, but not for areas where all reservoirs are totally ungauged. It would be particularly useful in data-sparse regions where current resources do not provide reliable monitoring of reservoir area, especially where estimation of evaporative loss is critical to the management and expansion of large water storage infrastructure. The results of this study highlight the heterogeneity and significance of reservoir evaporative loss. The amount of evaporative loss reflects the combined effects of various engineering, geological and meteorological factors. Understanding the spatiotemporal patterns of reservoir evaporation provides important information for designing effective evaporation suppression measures, evaluating water infrastructure performance and improving regional water resilience. Results of this study indicate that the magnitude of water loss through evaporation is comparable to water use of major economic sectors and ecosystem water requirements. It can be significant enough to counterbalance the positive effects

H. Zhang et al. / Journal of Hydrology 555 (2017) 213–227

of impounding water and offset the contribution of water conservation and reuse practices. Neglecting the critical impact of reservoir evaporation could result in serious overestimation of water availability and therefore underestimation of the required storage capacity to support the desired water release (Mays, 2011). Acknowledgments We thank the Google Earth Engine team for their generous technical support. We also thank Texas Water Development Board for the help with reservoir precipitation and evaporation data. We are grateful to the three reviewers for their insightful comments. This work was supported by Texas A&M University-Corpus Christi (TAMUCC) under Texas Comprehensive Research and Research Fund (TCRF) and by the National Science Foundation (NSF) under grant GEO/OAD-1342869 to Stanford University. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of either TAMUCC or NSF. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.jhydrol.2017.10. 007. References Assouline, S., Narkis, K., Or, D., 2011. Evaporation suppression from water reservoirs: efficiency considerations of partial covers. Water Resour. Res. 47 (7). https://doi.org/10.1029/2010wr009889. Burke, M., Lobell, D.B., 2017. Satellite-based assessment of yield variation and its determinants in smallholder African systems. Proc. Natl. Acad. Sci. U.S.A. 114 (9), 2189–2194. https://doi.org/10.1073/pnas.1616919114. Du, Y. et al., 2016. Water bodies’ mapping from Sentinel-2 imagery with modified normalized difference water index at 10-m spatial resolution produced by sharpening the SWIR band. Remote Sens. 8 (4), 354. https://doi.org/10.3390/ rs8040354. Duan, Z., Bastiaanssen, W.G.M., 2013. Estimating water volume variations in lakes and reservoirs from four operational satellite altimetry databases and satellite imagery data. Remote Sens. Environ. 134, 403–416. https://doi.org/10.1016/j. rse.2013.03.010. Farnsworth, R.K., Thompson, E.S., Peck, E.L., 1982. Evaporation Atlas for the Contiguous 48 United States. NOAA. 41 pp.. Feyisa, G.L., Meilby, H., Fensholt, R., Proud, S.R., 2014. Automated Water Extraction Index: a new technique for surface water mapping using Landsat imagery. Remote Sens. Environ. 140, 23–35. https://doi.org/10.1016/j.rse.2013.08.029. Gao, H., 2015. Satellite remote sensing of large lakes and reservoirs: from elevation and area to storage. Wiley Interdiscip. Rev. Water 2 (2), 147–157. https://doi. org/10.1002/wat2.1065. Gao, H., Birkett, C., Lettenmaier, D.P., 2012. Global monitoring of large reservoir storage from satellite remote sensing. Water Resour. Res. 48 (9). https://doi.org/ 10.1029/2012wr012063. Hansen, M.C. et al., 2013. High-resolution global maps of 21st-century forest cover change. Science 342 (6160), 850–853. https://doi.org/10.1126/science.1244693. Irish, R.R., 1998. Landsat 7 Automatic Cloud Cover Assessment. NASA, Greenbelt, MD, p. 8. Ji, L., Zhang, L., Wylie, B., 2009. Analysis of dynamic thresholds for the normalized difference water index. Photogramm. Eng. Remote Sens. 75 (11), 1307–1317. Krause, P., Boyle, D.P., Bäse, F., 2005. Comparison of different efficiency criteria for hydrological model assessment. Adv. Geosci. 5, 89–97. https://doi.org/10.5194/ adgeo-5-89-2005. Li, L. et al., 2015. Evaluation of MODIS spectral indices for monitoring hydrological dynamics of a small, seasonally-flooded wetland in Southern Spain. Wetlands 35 (5), 851–864. https://doi.org/10.1007/s13157-015-0676-9.

227

Li, W. et al., 2013. A comparison of land surface water mapping using the normalized difference water index from TM, ETM+ and ALI. Remote Sens. 5 (11), 5530–5549. https://doi.org/10.3390/rs5115530. Longley, W.L., 1994. Freshwater Inflows to Texas Bays and Estuaries: Ecological Relationships and Methods for Determination of Needs. Texas Water Development Board and Texas Parks and Wildlife Department, Austin, TX, p. 386. Lowe, L.D., Webb, J.A., Nathan, R.J., Etchells, T., Malano, H.M., 2009. Evaporation from water supply reservoirs: an assessment of uncertainty. J. Hydrol. 376 (1– 2), 261–274. https://doi.org/10.1016/j.jhydrol.2009.07.037. Martínez Alvarez, V., González-Real, M.M., Baille, A., Maestre Valero, J.F., Gallego Elvira, B., 2008. Regional assessment of evaporation from agricultural irrigation reservoirs in a semiarid climate. Agric. Water Manage. 95 (9), 1056–1066. https://doi.org/10.1016/j.agwat.2008.04.003. Masek, J.G. et al., 2006. A Landsat surface reflectance dataset for North America, 1990–2000. IEEE Geosci. Remote Sens. Lett. 3 (1), 68–72. https://doi.org/ 10.1109/LGRS.2005.857030. Mays, L.W., 2011. Water Resources Engineering. Wiley. McGloin, R. et al., 2014. Quantification of surface energy fluxes from a small water body using scintillometry and eddy covariance. Water Resour. Res. 50 (1), 494– 513. https://doi.org/10.1002/2013wr013899. McJannet, D. et al., 2013. Long-term energy flux measurements over an irrigation water storage using scintillometry. Agric. For. Meteorol. 168, 93–107. https:// doi.org/10.1016/j.agrformet.2012.08.013. Moriasi, D.N. et al., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50 (3), 885–900. Mu, Q., Zhao, M., Running, S.W., 2013. MODIS Global Terrestrial Evapotranspiration (ET) Product – Algorithm Theoretical Basis Document (Collection 5). The University of Montana. Mu, Q., Zhao, M., Running, S.W., Kimball, J.S., McDowell, N.G., 2016. Using MODIS weekly evapotranspiration to monitor drought. In: Gao, W., Chang, N.-B. (Eds.), SPID. Remote Sensing and Modeling of Ecosystems for Sustainability XIII, 9975, UNSP 997502. doi: 10.1117/12.2237749. Müller, M.F., Yoon, J., Gorelick, S.M., Avisse, N., Tilmant, A., 2016. Impact of the Syrian refugee crisis on land use and transboundary freshwater resources. Proc. Natl. Acad. Sci. U.S.A. 113 (52), 14932–14937. https://doi.org/10.1073/ pnas.1614342113. Pekel, J.F., Cottam, A., Gorelick, N., Belward, A.S., 2016. High-resolution mapping of global surface water and its long-term changes. Nature 540 (7633), 418–422. https://doi.org/10.1038/nature20584. Phillips, R.C., Saylor, J.R., Kaye, N.B., Gibert, J.M., 2016. A multi-lake study of seasonal variation in lake surface evaporation using MODIS satellite-derived surface temperature. Limnology 17 (3), 273–289. https://doi.org/10.1007/s10201-0160481-z. Powell, G.L., Matsumoto, J., Brock, D.A., 2002. Methods for determining minimum freshwater inflow needs of texas bays and estuaries. Estuaries 25 (6B), 1262– 1274. Rokni, K., Ahmad, A., Selamat, A., Hazini, S., 2014. Water feature extraction and change detection using multitemporal landsat imagery. Remote Sens. 6 (5), 4173–4189. https://doi.org/10.3390/rs6054173. Rost, S. et al., 2008. Agricultural green and blue water consumption and its influence on the global water system. Water Resour. Res. 44 (9). https://doi.org/10.1029/ 2007WR006331. TWDB, 2016. Precipitation & Lake Evaporation, Austin, TX. Available from: http:// www.twdb.texas.gov/surfacewater/conditions/evaporation/ (accessed 01.12.16). TWDB, 2017a. 2017 State Water Plan – Water for Texas, Austin, TX, 150 pp. TWDB, 2017b. Texas Water Use Estimates – 2015 Summary, Austin, TX. Warren, M.S., et al., 2015. Seeing the Earth in the Cloud: processing one petabyte of satellite imagery in one day. In: 2015 IEEE Applied Imagery Pattern Recognition Workshop (AIPR), pp. 1–12. doi: 10.1109/AIPR.2015.7444536. Wurbs, R.A., Ayala, R.A., 2014. Reservoir evaporation in Texas, USA. J. Hydrol. 510, 1–9. https://doi.org/10.1016/j.jhydrol.2013.12.011. Xu, H., 2006. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens. 27 (14), 3025–3033. https://doi.org/10.1080/01431160600589179. Yang, Y. et al., 2015. Landsat 8 OLI image based terrestrial water extraction from heterogeneous backgrounds using a reflectance homogenization approach. Remote Sens. Environ. 171, 14–32. https://doi.org/10.1016/j.rse.2015.10.005. Zhu, Z., Wang, S., Woodcock, C.E., 2015. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images. Remote Sens. Environ. 159, 269–277. https://doi.org/ 10.1016/j.rse.2014.12.014.