Using MODIS / Terra and Landsat imageries to improve surface water quantification in Sylhet , Bangladesh

Bangladesh has experienced multiple freshwater issues including salinization from monsoonal floods and groundwater over-pumping that induces severe land subsidence. Therefore, using satellite observations to virtually build a monitoring network becomes an efficient and innovative means. We focus on the Sylhet Mymensingh haor area that has the highest annual precipitation and the largest inundation area in northeastern Bangladesh. The modified normalized difference water index is first used to extract water area from MODIS and Landsat-5/-7/-8 optical imageries. A weekly flood chance model is then created from a sequence of images to recover water extent from the cloud-covered images. Using MODIS images for water identification achieves an overall accuracy of 84% in rainy season and 41% in dry season as validated with Sentinel-1A radar images. This model can be further used to refine the Shuttle Radar Topography Mission digital elevation model (DEM). As compared with ICESat laser altimetry, the root-mean-square of the height difference is improved from 1.65 m to 1.16 m after DEM modification. By combining the recovered water area and the refined DEM, surface water volume (WV) is quantified. A comparison with the Gravity Recovery And Climate Experiment (GRACE) gravimetry retrieved equivalent water heights (EWHs) in 2002 2015 is conducted, where the correlation coefficient and root-mean-square of the EWH difference are 91.7% and 0.09 m, respectively. Article history: Received 28 February 2018 Revised 8 November 2018 Accepted 15 November 2018


INTRODUCTION
In populated and low-lying deltaic countries, how to manage water resources is a necessary yet practically difficult issue.It is necessary to know how water cycle progresses and how much water resources are accessible all year round.Traditional methods to monitor water flow in a river or lake are the deployment of water level gauges and in situ wells.However, these methods are not only expensive and time-consuming, but also require sustained labor cost to establish and maintain the infrastructure.In addition, observational quality is not always guaranteed and the data streaming is sometimes interrupted (Andersen et al. 2008;St. Jacques and Sauchyn 2009).
Bangladesh has the 10 th highest population density in the world with a remarkable number of over 1100 people km -2 in 2000 (Neumann et al. 2015).It is located south of the Tibetan Plateau with a low-lying land formed by the Ganges-Brahmaputra-Meghna floodplain (Fig. 1).Hence, its main water resources are acquired from surface flow of these three major rivers.However, although the Indian monsoon and melting snow coming from the Himalayas have brought a considerable amount of surface runoff every year, water flow in river channels is almost impossible to retain by building reservoirs, dams or other water storage facilities due to the flat terrain -most of land elevation is below 10 m countrywide.As a result, Bangladesh is one of the most vulnerable countries that suffer from serious freshwater problems.Periodic and flash floods owing to sudden increase in water volume have caused disastrous damages in many districts, for example, in the Sylhet Mymensingh haor (wetland) area (orange star in Fig. 1) (Islam and Sado 2000).
Due to the difficulties in capturing surface water, Bangladesh relies heavily on groundwater for its municipal and agriculture usages.Unfortunately, the groundwater has been over-pumped in the last few decades and has caused serious land subsidence problems (Hoque et al. 2007).Besides, the Holocene strata in Bangladesh contains mass arsenic poisoning; while people drill water wells to acquire groundwater, the entire aquifer is contaminated by chemicals (Smith et al. 2000).Thus, it is necessary for water management agencies to monitor the migration of surface and subsurface water resources.But the fact is that a water monitoring network is not complete in Bangladesh.The number of water gauges in the entire country is only 85 as of today (Bangladesh Water Development Board 2018, http://www.ffwc.gov.bd/), while lots of idle wells are not used for monitoring groundwater.
To observe basin-scale changes of hydrological regime, many researchers are now seeking possibilities from spaceborne sensors.As reported in the literature, several satellite hydrology applications used a synergy of optical remote sensing, synthetic aperture radar (SAR), pulse-limited radar/SAR altimetry, and satellite gravimetry to monitor three basic parameters: water area (WA), water level (WL), and water volume (WV) (Henry et al. 2006;Hostache et al. 2009;Liu et al. 2016;Tseng et al. 2016a, b).First, the measurement of WA can be achieved by thematically classifying optical or radar satellite images, such as the normalized difference water index (NDWI) (McFeeters 1996) or the modified NDWI (MNDWI) (Xu 2006) for optical (visible to mid-infrared) sensors.For SAR images, one could apply a single-threshold cutoff in the histogram of backscattering intensity for water identification (Hostache et al. 2009;Li and Takeuchi 2016).Second, the WL can be acquired from satellite altimetry or combined satellite images and digital elevation models (Maheu et al. 2003;Berry et al. 2005;Schumann et al. 2008;Hostache et al. 2009;Kuo and Kao 2011;Tseng et al. 2016b).Although satellite altimetry is prone to land contamination in the radar echogram over inland waterbodies, there are many studies aiming at developing customized waveform retracking algorithms to reduce unmodeled effects (Legresy et al. 2005;Frappart et al. 2006).Once both WA and WL parameters are obtained, one can combine them to calculate WV.However, the number of studies focusing on regional surface WV is quite limited.It is because most methods used to measure WA and WL have some restrictions.For example, optical satellite has cloud cover and aerosol contamination issues.In contrast, radar images contain too much noise due to similarity in backscattering over a complex of surface types.For altimetry data, although radar waveforms can be retracked, the uncertainty is still too large to calculate WL for narrow rivers.
Another popular approach to study basin-scale hydrology is the Gravity Recovery And Climate Experiment (GRACE) satellite gravimetry mission.Its main goal is to observe Earth's time-variable gravity field from subtle range changes between a twin-satellite design.The temporal gravity field anomaly (monthly mean minus long-term static background) measured by GRACE is composed of changes of total column water including surface and subsurface contributions (Rodell et al. 2007;Andersen et al. 2008;Feng et al. 2013;Rateb et al. 2017).It reveals a possibility to observe groundwater transition from spaceborne measurements.Once surface water component been well observed or modeled, the groundwater variation and its long-term tendency can be estimated.However, the accuracy of quantifying surface WV strongly influences the estimate of groundwater trend.How to accurately measure surface WV has become an important issue especially in the heavily flooded area.
Considering the advantage and limitation of each satellite sensor discussed above, the goal of this research aims to develop a surface water quantification method by combining multiple satellite data.We first use historical optical satellite imageries including the MODerate-resolution Imaging Spectroradiometer (MODIS) onboard Terra/Aqua and Landsat family to create weekly flood chance models to bypass the cloud-cover problem.This model is a historical statistics of inundation chance used to refer the "relative elevation" in an enclosed basin, where the higher chance means lower elevation and vice versa.This model is later used to fill the cloud gaps left by the MODIS MOD09A1 product.After the cloud-pixel recovery process, European Space Agency's (ESA's) Sentinel-1A radar images is utilized to crosscheck water area predicted by this model.Next, we use a decadal overall flood chance model to modify the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) in Sylhet to quantify surface WV.Finally, the surface WV is compared with GRACE equivalent water height (EWH) product.
There are five sections in this paper.Section 1 introduces the background and objective of this study.The study area is introduced in section 2. Section 3 explains the workflow, satellite datasets, and data processing methods.In section 4. The results of image-derived DEMs are validated with ICESat data and the estimated WV are compared with GRACE.Finally, the conclusion is given in sections 5.

STUDY AREA
Sylhet Mymensingh haor area is located south of the Shilong Plateau of northern Indian shield.Shilong belongs to the Meghalaya state of India, which literally means "The Abode of Clouds".The particular geographic landscape of Shilong makes it the first orographic barrier encountered by the humid southwesterly monsoon winds (Prokop 2014).Hence, annual precipitation is anomalously higher than the regional average.In Fig. 2, the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) monthly release 3B43 version 7 (Huffman et al. 2007) is used to depict the annual total precipitation.We adopt 20-year grid data (0.25° by 0.25°) from 1998 -2017 and take an average from its hourly rate to approximate the annual accumulation.As seen in the contours of Fig. 2, annual accumulation of rainfall reaches more the 5500 mm in Shilong Plateau, which plays as a major cause of periodic flood in Sylhet.In fact, the annual total rainfall in Bangladesh is between 1500 and 5000 mm on average, with an increasing tendency from west to northeast.Figure 2 also displays flat terrain in this region from the Shuttle Radar Topography Mission (SRTM) DEM.Sylhet floodplain is located in the hotspot of high precipitation area.Accordingly, the rise of Brahmaputra River and the overflow of Meghna River during wet season would likely to worsen the flood situation.

Workflow
The data used in this study include MODIS and Landsat optical satellite image, SRTM DEM, and GRACE observations.Besides, we use Sentinel-1A images and ICESat elevation products for water area and DEM validation, respectively.The workflow from data preprocessing to estimating surface WV is shown in Fig. 3.More in details are described in the following subsections.

Optical Remote Sensing for Water Identification
Using optical satellite images, for example MODIS and Landsat, to monitor surface water is a low-cost, efficient, and accurate method for water resource management.Previous studies have demonstrated the usefulness, reliability, and repeatability of these datasets (McFeeters 2013;Liu et al. 2016;Tseng et al. 2016b).In this research, MODIS MOD09A1 product and Landsat-5/-7/-8 images are used to identify surface WA.
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a remote sensing payload onboard the Earth Observing System (EOS) Terra and Aqua satellite, jointly operated by the National Aeronautics and Space Administration (NASA) and U.S. Geological Survey (USGS).Since Terra (EOS-AM) and Aqua (EOS-PM) launched in 1999 and 2002, respectively, these two satellites have provided numerous multispectral images for large scale observation of atmospheric and terrestrial environment.With 36 bands in visible, near/mid infrared, and thermal infrared, MODIS has been proven extremely useful in monitoring aerosol properties (Tanré et al. 1997), vegetation (Huete et al. 2002) and also surface water variation.The spatial resolution of MODIS varies with its spectral bands.For band 1-2 the spatial resolution is 250 m, while band 3-7 has 500 m resolution, and band 8-36 has the lowest resolution at 1 km.The revisit time of Terra and Aqua is 1-2 days, which enables a global coverage in a few days.In this study, the surface reflectance (SR) products are needed to detect water pixels in the wetland.MODIS SR products (MOD09) provide various temporal and spatial resolutions, from daily in 250 -1000 m grid to 8 days with 250 -500 m resolution.Because of the frequent and heavy cloud in Sylhet Mymensingh haor, high-temporal resolution is not useful in this regard.Hence, the 8-day composite products MOD09A1 and MOD09Q1 composed of cloud-free pixels in each 8 days are better to acquire surface information.However, because MOD09Q1 in 250 m contains only band 1 and 2 without shortwave infrared, MOD09A1 in 500 m is finally chosen for this research.All MOD09A1 files in Hierarchical Data Format (HDF) from 2000 -2017 are downloaded from NASA Land Processes Distributed Active Archive Center (LP DAAC) data pool.Sylhet wetland is located in the sinusoidal grid number h26v6.After gathering all available images, the sinusoidal projection of raster data is converted into WGS84 coordinate by using the MODIS Reprojection Tool (MRT).
In addition to the coarse grid information given by MODIS, we adopt more fine-granule Landsat images to provide details in surface inundation.Landsat is a well-known series of Earth resource exploring satellites.Since Landsat-1 launched in 1972, this series of NASA/USGS joint mission has provided the longest and continuous land observations.As of today, Landsat-7/-8 are still working to provide images for environmental analysis seamlessly around the globe.In this research, the image taken from the Thematic Mapper (TM) onboard Landsat-5 in 2000 -2013, the Enhanced TM Plus (ETM+) onboard Landsat-7 in 2000 -2017, and the Operational Land Imager (OLI) of Landsat-8 in 2013 -2017 are used as model training data.Each of these three satellites has a 30 m spatial resolution and 16-day repeat cycle.Higher spatial resolution can help on extracting more detailed surface water, especially for wetland channels.Unfortunately, due to the Scan Line Corrector (SLC) failure in 2003, Landsat-7 images observed thereafter contain null strips that makes some information lost.But the images still provide partial information for surface water extraction.We simply apply a linear interpolation to remove gaps at the step after the binarized water classification.
In 2017, USGS reprocessed Landsat images and offered the "Collection 1 Level-1/-2" data.The Level-2 surface reflectance product has been atmospheric-corrected by either the 6S radiative transfer model (Vermote et al. 1997) for Landsat-4/-5/-7 or by the internal algorithm for Landsat-8 (Vermote et al. 2016).Landsat Collection 1 Level-2 GeoTiff images are accessible in the USGS EarthExplorer website.The Digital Numbers (DNs) in raster format are converted into surface reflectance by a simple rescaling factor.Following that, we adopt MNDWI (Xu 2006) to extract surface water owing to its efficiency in delineating water area from surrounding sandbar and manmade buildings.MNDWI was evolved from NDWI (McFeeters 1996), which adopts IR as weak reflectance band for water.Based on a few studies, MNDWI outperforms NDWI especially in vegetated areas (Singh et al. 2015) and had been used to map flood area in Cambodia (Ledien et al. 2017) and coastal Bangladesh (Ghosh et al. 2015).MNDWI uses the band contrast of green band and middle infrared band because water body shows a strong reflectance in green and higher absorption in middle infrared.These two bands are thus sensitive to check water existence.The formula of MNDWI is shown in Eq. ( 1) (Xu 2006).

MNDWI Green MIR
Green MIR where green band is MODIS band 4, Landsat-5/-7 TM/ ETM+ band 2 or Landsat-8 OLI band 3. MIR band is MODIS band 6, Landsat-5/-7 TM/ETM+ band 5, or Landsat-8 OLI band 6.A zoom-in view of MNDWI map of Landsat near Sylhet Mymensingh haor area (orange star in Fig. 1) is shown in Fig. 4, where the warm color means the stronger signal potentially contributed by water surfaces.A threshold is set manually by confirming with original truecolor images.The threshold is set as 0.2 for MODIS and 0 for Landsat series.

Seasonal Inundation Model from Water Area Time Series
As mentioned, one major restriction hampering satellite optical imageries to continuously monitor surface water in Sylhet is the persistently cloud-cover problem.A method is proposed here to potentially recover cloud-covered areas by using the flood chance model of each week.We use 771 MODIS MOD09A1 images between 2000 and 2017 (18 years) and 487 Landsat-5/-7/-8 TM/ETM+/OLI images in Sylhet area, which has seasonal floods every year to build weekly flood chance models.
All images in the same week of year are stacked or summed up and then divided by the number of images to get the flood chance of each pixel on a weekly basis, as shown in Eq. ( 2).
( , ) ( , ) ( , ) % ( , ,..., ) where P m (i, j) is the probability (0 -100%) of flooding in the m th week in row i and column j. ( , ) ( , ) max where T k in % is the flood chance threshold of the k th satellite optical image.( , ) S i j m k is the predicted map of water area from a cloud-cover image belonged to the m th week.
The abovementioned calculations are programmed and executed in MATLAB.

Sentinel-1A SAR Image for Validation
To verify the accuracy of the recovered result, three Sentinel-1A images are collected as validation data.Sentinel-1A is the first Sentinel-1 constellation launched in 2014 as part of European Union's Copernicus program, which was earlier named the Global Monitoring for Environment and Security (GMES) mission.The peer satellite Sentinel-1B was launched in 2016.The goal of GMES is to monitor the changes of ocean and land in more details without weather restrictions.With a C-band SAR sensor onboard, Sentinel-1 satellites provide high spatial (~10 m) and moderate temporal resolution in 12 days (6 days for Sentinel-1A and 1B constellation).Sentinel-1 has different modes in multiple resolution and swath.In this research, we adopt the Interferometric Wide (IW) mode VV/VH polarization Level-1 Ground Range Detected (GRD) product as validation data.The data is downloaded from ESA Copernicus data portal and processed by the Sentinel Application Platform (SNAP) complementary software.
The data preprocess includes four steps: radiometric correction, geometric correction, speckle filtering, and image subset.The radiometric correction aims to remove the effects of slant range, incidence angle and illumination pattern of antenna.After correcting radiometric noise, the DNs value of Sentinel-1A image can be converted into intensity as expressed in Eqs. ( 5) -( 7).
( ) log sin I 10 where i v is the backscatter coefficient in dB, i b is the radar brightness in dB, and I i is the incidence angle in degree.All these three parameters are orientationally dependent in the i th pixel.i b is derived by Eq. ( 6), where A3 is offset and A2 is the gain value.In Eq. ( 7), R e is the radius of Earth in km, H is the orbital height in km, i z is the spherical angle in degree, and R i is the slant range in km.
Next, the geometric correction step aims to correct surface deformation caused by the slanting photography of radar images.We apply SRTM 1 arcsec DEM to correct the distortion.The third step is speckle filtering, which removes the speckle noise caused by the interference of radar signal.Finally, the Lee filter (Lee 1980) is used to keep features of image after smoothing and also to eliminate the remaining noise.
After pre-processing, a histogram can be formed by the intensity values of an image, as shown in Fig. 6.In this histogram, a double-peak curve formed by the water and non-water components behaves distinguishable scattering characteristics.Since water-like features with relative smooth surfaces show relatively low backscatter coefficient as compared to other objects, most of the pixels in the left peak are formed by water pixels.Here, we simply use a two-degree Gaussian function to fit the histogram and set a threshold to extract water area, as expressed in Eq. ( 8) (Hostache et al. 2009).

Threshold mean
Three examination results for the recovered inundation areas in wet/dry seasons using MODIS and Landsat are given in Figs.7 -9.Figure 7a is the original cloud-covered MODIS 8-day composite (2015/08/13) after MNDWI thresholding, where red patches indicate water area, blue means non-water, and white represents cloud pixels.Figure 7b is the recovered result using flood chance model of that week.A Sentinel-1A image near the same time (2015/09/09) is obtained to delineate water area as shown in Fig. 7c.
By comparing Figs.7b and c, we notice that most water patches (red) are successfully predicted from the help of weekly model, even if the original image is heavily covered by cloud at ~50%.The overall accuracy of recovered water area is 84% as compared against Sentinel-1 observation, where the recovered and observed water area is 2.17 × 10 3 km 2 and 2.39 × 10 3 km 2 , respectively.
Another cloud-covered MODIS image (2016/12/14) is demonstrated in Fig. 8a for dry season.By using flood chance model for WA recovery, Fig. 8b removes cloud patches and suggests potential WA pixels.A Sentinel-1 image taken on 2016/12/15 (Fig. 8c) is used to verify water pixels.The overall accuracy of recovered water area is 41% as compared against Sentinel-1 observation, where the recovered and observed water area is 6.11 × 10 2 km 2 and 5.97 × 10 2 km 2 , respectively.A low accuracy appearing in this case is because the coarse resolution of MODIS data is not sensitive to narrow wetland channels.Once in dry season when the water bodies are small and sparse, the predicted water patches are overestimated owing to spectral mixing.It should be noted that the performance of flood chance model is determined  by the resolution of input data.Because MODIS in 500 m is wider than many narrow rivers or small water bodies, once a pixel in the input image classified as water has a pixel size large than the actual waterbody, the resampled 30 m grid will cover "high elevation" or "low flood chance" pixels in the flood chance model.Consequently, the threshold T k in Eq. ( 3) will be lowered and many high elevation pixels will be classified as water.However, it would less harmful because the absolute WA error is also small in dry season.The water volume quantification would be less affected by the WA estimating error.
In order to proof the usefulness of flood chance model, another cloud-covered Landsat image (2016/03/19) is tested as shown in Fig. 9a.By using the flood chance model for this heavily cloud-covered (85%) image, Fig. 9b clears cloud patches and suggests potential WA pixels.A Sentinel-1 image taken on 2016/03/07 (Fig. 9c) is used to verify water pixels.The overall accuracy of estimated water area is 96% as compared against Sentinel-1 observation, where the estimated and observed water area is 2.42 × 10 2 km 2 and 2.18 × 10 2 km 2 , respectively.This case successfully demonstrates the possibility of WA recovery from a retrospective analysis of surface conditions.

GRACE Temporal Gravity Field Terrestrial Water Storage
GRACE is a joint mission of NASA and the German Aerospace Center (DLR) since 2002.This satellite is composed of two identical spacecrafts separated along track by 220 km.By using the K-Band Ranging System (KBR) onboard each satellite, the distance change as small as 10 μm over 220 km is measured.The distance change due to orbit perturbation is caused by the mass inhomogeneity of the Earth, where the perturbed orbit can be converted into the temporal gravity field anomaly.It can be further used to study the terrestrial water storage change in terms of equivalent water height (EWH) change.The data set represents a long-term change of the Earth's gravity field or its mass change including, for example, seismic deformation, geodynamics, ice ablation, precipitation, hydrologic cycles, and mass change of the ocean.By removing the long-term observed or inverted gravity field from each individual monthly solution, the result is the temporal gravity field anomaly sampled every month.
The GRACE Level 2 data products are produced by official data centers including Center for Space Research (CSR) at The University of Texas at Austin, NASA/Jet Propulsion Laboratory, and German Research Centre for Geosciences (GFZ).We choose CSR Release (RL) 05 monthly EWH product to obtain temporal gravity field changes.This product is originated from CSR release (RL) 05 monthly 60 degree/order spherical harmonics product that has an approximate 330 km spatial resolution.Since GRACE ob-servations have errors that result in spatial "stripe" errors in north-south direction, all GRACE products need to be spatially filtered or post-processed to mitigate these errors.CSR RL 05 product is smoothed by the Decorrelation Filter (DDK).DDK 1 filter is a non-isotropic smoothing method proposed by Kusche (2007) to reduce the non-isotropic noise in GRACE data.Based on this method, Kusche et al. (2009) modified the smoothing method and named it as DDK 5.As an improvement over DDK 1 filter, DDK 5 can be applied to higher resolution with far less unsolved coefficients in the decorrelation process (e.g., Duan et al. 2009;Guo et al. 2010).The CSR RL 05 data product can be acquired from the GRACE Plotter website provided by the Centre National d'Etudes Spatiales/Groupe de Recherches de Géodésie Spatiale (CNES/GRGS) of France.The unit of temporal gravity field change has been converted into EWH (Wahr et al. 1998).EWH is used to represent mass change of a unit area as totally covered by water.To compare with the GRACE derived EWH, we divide estimated surface WV by a 330 km × 330 km grid corresponding to the GRACE spatial resolution.

Image-Derived DEM and WV Quantification
After recovering WA from all optical satellite images, WL needs to be determined for the quantification of surface WV.We utilize SRTM v3 global 30 m product released in 2015 as a height reference.SRTM provides orthometric height based on EGM96 geoid and covers terrestrial area from 56°S -60°N.This model is also available in the USGS EarthExplorer data portal.SRTM has been used in many studies, however, 90% random error is about 5 m in Bangladesh (Rodríguez et al. 2005), while most terrains in Bangladesh is less than 10 m.It implies that the relative error of SRTM in Bangladesh may still be high.Here we use our developed method to demonstrate whether the referenced SRTM could be improved by removing possible errors or outliers.
As mentioned, the flood chance model can represent terrain information because the boundary of an enclosed waterbody would have a unique orthometric elevation.Therefore, we assume that the pixels with the same flood chance have the same elevation.This method is similar to the hypsometry approach widely adopted in hydrologic research to recover WL of waterbodies.With this presumption, the overall flood chance model is first used to define multiple independent wetland watersheds (Fig. 10a).In this step, pixels with flooding chances less than 5% are ignored to avoid subtle noise in small waterbodies, such as agriculture zones.Next, we convert raster data into multiple polygons, where each polygon means an independent watershed.Each independent watershed is used to create a mask for the co-registered SRTM and flood chance model.For each watershed, we set the highest and lowest elevation by averaging SRTM pixels with the same boundary flood chance (lowest and highest flood chance).The elevation of remaining pixels in the same wetland are linearly interpolated based on flood chance order.The computation follows a similar work (Tseng et al. 2017) applied for coastal DEM reconstruction: where DEM sat (i, j) is the satellite image-derived DEM, P fs (i, j) is the feature-scaled probability of P m (i, j) in Eq. ( 2).H h is an average of the highest 2% (to avoid extremely values) of SRTM elevation within the co-registered water patch.H l is an average of the lowest 2% of SRTM elevation within the co-registered water patch.
To verify the image-derived DEM, ICESat surface elevation (GLA14, release-33) is used for comparison.ICESat is a laser altimetry mission operated in campaign mode during 2003 -2009.The elevation of the illuminated surface can be precisely retrieved from the photon counting technique (Abshire et al. 2005;Urban et al. 2008).Here, we choose 6 out of 18 operating periods in dry seasons (February to April) for the validation of surface height.The original SRTM is displayed in Fig. 10b, in which a few ICESat ground tracks are overlaid as black lines.A comparison of image driven DEM and SRTM is shown in Figs.10c -j.Figures 10c, e, and g are the image-derived DEM based on MODIS, Landsat, and a combination of two, respectively.Figures 10d, f, and h are the scatter plots of derived DEM (red dots) and SRTM (black dots) against ICESat measurements in geoidal heights.For all these comparisons, the SRTM and satellite images have been resampled to 30 m and co-registered.Meanwhile, all DEM results are linearly interpolated from the longitude and latitude of ICESat ground points.The statistical results of comparison are shown in Table 1.In Figs.10c -d, we find that the MODIS-derived DEM appears a positive bias as compared with ICESat.However, it has a smoother surface as seen in the continuous distribution of red dots in Fig. 10d.In contrast, SRTM with integer number of elevations appears a vertical-strip pattern.The root-mean-square of the difference (RMSD) and correlation coefficient (CC) is 1.65 m/0.72 for SRTM and 1.48 m/0.75 for MODIS.For Landsat-derived DEM as shown in Figs.10e -f, the bias is reduced and the uncertainty level is a little smaller than that of SRTM.The RMSD/CC is 1.37 m/0.75 for Landsat.When we merge two datasets altogether to calculate the inundation chance and create a merged DEM, the result in Figs.10g -h shows an even better result, with RMSD/CC equal to 1.16 m/0.79.It is noticed that the noisy SRTM has been smoothed and improved by using this method, regardless of the image resolution in this case.
After the modification of SRTM, extracted water area in MODIS images can be co-registered with the merged DEM to compute water volume.Here, the water level is assigned to each watershed by finding the highest elevation within the co-registered water patch, and column water is integrated as the volume between DEM and water surface.

Comparison and Validation Using GRACE EWH
The comparison of estimated surface water volume and GRACE data in Sylhet haor area is shown in Fig. 11.Here, the MODIS product is linearly interpolated to synchronize the timing with GRACE data in the middle of each month.This plot demonstrates the capability of the developed method to capture the transportation of surface water storage.The amplitude of annual signal is similar in two independent approaches for WV quantification.The correlation coefficient in between is 0.91, where the root-mean-square of the difference (RMSD) is 0.09 m in EWH during overlapping years (2002 -2014).Also, two major floods happened in 2004 and 2007 are observed in the estimated WV.They were barely seen in the studies using remote sensing and altimetry satellites owing to cloud cover or the lack of ground tracks, respectively.Our work is one of a few in the literature that can estimate temporal water volume changes solely from the optical remote sensing imageries and SRTM.

Discussion About Image Sources for Deriving DEM
In the current form of our workflow, there are only optical imageries considered as the sources for DEM production, because a sufficient amount of water areas in different levels needs to be collected to form an inundation chance map and to further generate DEM.However, one major disadvantage is that we need to fill cloud gaps in each image on the stage of water identification.This step may introduce commission/omission errors in predicted water patches due to DEM uncertainty.Thanks to the open data policy of Sentinel-1 (Torres et al. 2012), this problem could be mitigated in the future while using frequent Sentinel-1A/B constellation with 6-day revisit, or even shorter if both ascending and descending passes are included.SAR images will provide faster inundation assessment (Borah et al. 2018) with high accuracy (Huang et al. 2018).Sentinel-1 SAR will enormously increase the efficiency in water identification (Palmer and Ruhi 2018).In other words, it requires short time span to collect sufficient number of images for weekly flood chance model, so that the dependence on optical images could be reduced.Also, more images-recognized water area (WA) means that we can estimate water level and water volume in a frequent manner.In our future work, we will use weekly SAR images, or at a higher sampling rate, to investigate flood-induced rapid surface water migrations.

CONCLUSIONS
In this paper, we propose an effective method to improve the quantification of surface water by integrating historical optical satellite imageries and convert the water appearance chances into relative orthometric heights.The study mainly focuses on Sylhet floodplain which has the highest annual precipitation in Bangladesh.During monsoon season, the recovered water area shows 84% overall accuracy as compared with Sentinel-1A classification result.The accuracy may be less than 50% in dry season if the study area contains many narrow tributaries that is hardly detected by MODIS MOD09A1.But the accuracy improves as an increasing resolution of input sources, such as Landsat 30 m imageries.The method of flood chance model is useful considering its ability of extending surface WA measurements, especially in the heavily cloud-covered area in a long rainy season.The model accuracy can be improved if we keep adding satellite images with higher resolution sources, for example Landsat-7/-8, Sentinel-2 and even SAR satellites.
Furthermore, the flood chance model also provides a means to reduce the uncertainty in SRTM DEM, and actually improve the model.The RMSD improves from 1.65 m in SRTM to 1.16 m in the merged (MODIS + Landsat) DEM as compared with ICESat altimetry results.The validation indicates that the method developed in this study smooths out the stepwise SRTM in integer meters.With the recovered water area and reconstructed DEM, a continuous quantification of surface water volume can be derived.The signals of two extreme flood events in 2004 and 2007 are observed in the curve of estimated surface WV.From the comparison with GRACE EWH, the water storage changes in Sylhet Mymensingh haor can be modelled with a good spatiotemporal resolution.Assuming other mass changes, such as tectonics and groundwater variability are ignorable in this area, a majority of mass changes caused by surface water is well explained by using the introduced workflow.
In conclusion, this research develops a monitoring system which relies only on the low-cost data to the end users.It effectively monitors surface water volume on a basin scale.Owing to the accessibility of abundant satellite images nowadays, this monitoring system is transportable to other study areas of the world.Moreover, the open data of Sentinel-1 may also serve as an important image source in the flood chance model as the number of images accumulated.A better estimation of surface WV may further be used for investigating subsurface water changes.With more and more satellite missions being launched, it is expected that this monitoring approach could become an efficient tool for hydrological studies and water resources management.Ministry of Science and Technology, Taiwan under project number 104-2119-M-008-005-MY2 and 104-2221-E-008-099-MY3 and in part, by the National Central University New Faculty Research Award.We thank NASA, USGS, and ESA for providing satellite imageries used in this study.

Fig. 1 .
Fig. 1.Political boundaries of Bangladesh and three major rivers.Orange star is the center of our study area.Orange box is a reference boundary approximating Gravity Recovery And Climate Experiment (GRACE) gravimetry data derived equivalent water height (EWH) change at 330 km spatial resolution.Red dots indicate major urban areas.

Fig. 2 .
Fig. 2. SRTM elevation of Bangladesh and neighboring countries.Contours are the average annual precipitation in mm calculated from TRMM 3B43 v7.
Fig. 4. A sample of MNDWI map from a single Landsat image.

Fig. 10 .
Fig. 10.(a) Isolated wetland basins in different colors.(b) Original SRTM 30 m DEM in Sylhet area, with ICESat tracks (black line) in dry season for validation.(c) MODIS-derived DEM.(d) Comparison between MODIS-derived DEM and ICESat estimates.(e) Landsat-derived DEM.(f) Comparison between Landsat-derived DEM and ICESat estimates.(g) Merged (MODIS + Landsat) DEM.(h) Comparison between merged DEM and ICESat estimates.(i) A sample of river channel in SRTM.(j) A sample of river channel in the merged DEM.Panel (b), (c), (e), (g), (i), (j) share the same color bar as in (b).

Fig. 11 .
Fig. 11.Time series of estimated surface water and GRACE EWH in Sylhet.Blue line is the estimated surface water EWH from merged DEM and satellite images.Green line is GRACE-observed EWH.