**Geo-Spatial Analysis in Hydrology**

Editors

**Qiming Zhou Jianfeng Li**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Qiming Zhou Hong Kong Baptist University Hong Kong

Jianfeng Li Hong Kong Baptist University Hong Kong

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *ISPRS International Journal of Geo-Information* (ISSN 2220-9964) (available at: https://www.mdpi. com/journal/ijgi/special issues/spatial hydrology).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-980-5 (Hbk) ISBN 978-3-03936-981-2 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

**Qiming Zhou** (Professor) is a Professor of Geography and Associate Dean of Faculty of Social Science (Research) and Director of Centre for Geo-Computation Studies at Hong Kong Baptist University. His research interests cover a broad area of geo-spatial information science, particularly in geo-computation and remote sensing applications. He has been actively engaged in research such as digital terrain analysis, landuse and land cover change detection, monitoring and modeling, spatio-temporal pattern analysis in multi-temporal remote sensing, and GIS and remote sensing applications to urban, environmental, and natural resource management. Prof Zhou has published widely in fields of GIS, terrain analysis, remote sensing change detection, and spatio-temporal modelling. He currently holds visiting professorships at four universities/research institutions, and serves as Asia and Australasia Editor of Transactions in GIS and Editorial Board members of six academic journals.

**Jianfeng Li** (Assistant Professor), Dr., is an Assistant Professor in the Department of Geography, and the deputy director of the Centre for Geo-computation Studies of Hong Kong Baptist University. His major research interests include hydroclimatology, environmental change, and water hazards, focusing on climate change impacts on hydrological processes and the environment especially hydrological and climatic extremes. He is an Associate Editor of Hydrological Processes, an international journal in hydrology, and has been active in serving in various professional communities. His studies have been published in top-tier journals, including Nature Climate Change, PNAS, Journal of Hydrometeorology, and Journal of Geophysical Research: Atmospheres.

### *Editorial* **Geo-Spatial Analysis in Hydrology**

#### **Qiming Zhou \* and Jianfeng Li**

Department of Geography, Hong Kong Baptist University, HKSAR, Hong Kong, China; jianfengli@hkbu.edu.hk **\*** Correspondence: qiming@hkbu.edu.hk

Received: 9 July 2020; Accepted: 9 July 2020; Published: 11 July 2020

**Abstract:** With the increasing demand for accurate and reliable hydrological information, geo-spatial analysis plays a more and more important role in hydrological studies. The development of the geo-spatial technique advances our understanding of the complex and spatially heterogeneous hydrological systems. Meanwhile, how to efficiently and effectively process and analyze multi-source geo-spatial data has become more challenging in the fields of hydrology. In this editorial, we first review the development and application of geo-spatial analysis in three major topics in hydrological studies, namely the scaling issue, extraction of basin characteristics, and hydrological modelling. We hence introduce the articles of the Special Issue. These studies present the latest results of geo-spatial analysis in different topics in hydrology, and improve geo-spatial analytic methods for better accuracy and reliability.

**Keywords:** hydrology; geo-spatial analysis; scaling issue; basin characteristic extraction; hydrological modelling

#### **1. Introduction**

With the rapid development of geo-spatial technology and the increasing demand for accurate and reliable hydrological information, geo-spatial analysis plays a more and more important role in hydrological studies, such as hydrological monitoring, water resources assessment, and water-related decision making. The applications of geographical data acquiring, storing and analytic technique in hydrology provide more detailed (e.g., finer resolution) and more extensive (e.g., larger spatial and temporal coverage) information of the distribution, movement, and dynamic of water components, such as precipitation, surface runoff, soil moisture, and groundwater. These geo-spatial data and techniques advance our understanding of the complex and spatially heterogenous hydrological systems especially in ungauged regions such as arid and semi-arid areas. Meanwhile, how to efficiently and effectively mine, process and analyze multi-source geo-spatial data has become more challenging in the fields of hydrology.

Geo-spatial analytic methods have been widely adopted in hydrological studies. The applications of geo-spatial analysis in three fundamental issues in hydrology are briefly reviewed in this editorial, namely the scaling issue, extraction of basin characteristics, and hydrological modelling.

(1) Scaling issue of hydrometeorological variables

Traditional in situ measurements of hydrometeorological variables, such as precipitation, temperature, streamflow, and soil moisture, are at the site scale, while some hydrological studies require basin-scale analysis [1]. Therefore, how to solve the scaling mismatch issue has been one of the most fundamental issues in hydrology [2]. Various spatial interpolation techniques have been developed to interpolate point-scale observations to the areal scale, such as Thiessen Polygons, Inverse Distance Weighting, and Ordinary Kriging [3]. Hydrological variables such as precipitation are highly heterogeneous in space and time, and their spatial distributions are strongly affected by topographic factors. In order to improve the skill of the interpolation methods, especially in the regions with complex terrain, topographic factors are considered as secondary predictors in certain

interpolation methods, such as Thin Plate Smoothing Splines [4] and Kriging with External Drift [5]. Besides spatial interpolation techniques, Areal Reduction Factor, a ratio between the areal average value and the point-scale value, is another popularly used strategy to deal with the scaling issue of precipitation in hydrological risk analysis and design [6]. Moreover, different types of statistical downscaling methods have been developed to downscale future projections from coarse Global Climate Models to a fine resolution [7,8]. Statistical downscaling methods are based on the statistical relationship between values at a coarse resolution and a fine resolution in a selected training period.

(2) Extraction of basin characteristics

Basin characteristics, such as elevation, slope, basin area, river channel, and land cover land use, are the most basic and essential information for hydrological studies such as hydrological modelling [9]. Geo-spatial methods have been developed to extract these key characteristics from grid-based Digital Elevation Models (DEM) [10]. For instance, to acquire terrain parameters for hydrological studies at various spatial scales, DEM may need to be generalized or reconstructed using geo-spatial techniques such as triangulated irregular network (TIN). However, the reliability of basin characteristics extracted depends on DEM accuracy, data structure, and derivation algorithm [11,12]. For example, flow routing extracted from DEM can be affected by the choice of extraction algorithms [13]. A number of studies have been conducted to improve the accuracy and reliability of basin characteristics extracted from geo-spatial data. Zhou and Chen [9] proposed a compound method to integrate the point-additive and feature-point methods to construct a TIN. The compound method is capable of keeping the major terrain features while significantly reducing the elevation data points. To overcome the complication of flow divergence/convergence in traditional raster-based methods for estimating surface flow paths, Zhou et al. [14] developed a triangular facet network algorithm. They found that the facet-based algorithm outperforms the traditional methods with better representation and more consistent outcomes.

#### (3) Hydrological modelling

The advancement of geo-spatial techniques and the increasing availability of multi-source geographical data facilitate the development of hydrological models. Remote sensing provides a variety of hydrometeorological (e.g., precipitation, evapotranspiration, and soil moisture) and land surface data (e.g., DEM and land use land cover) with extensive spatial coverage for hydrological models as inputs or as reference for model calibration and validation [15,16]. A data assimilation technique has been developed to couple instrumental observations (more accurate but lacking spatial representativeness) and hydrological models (better spatial representativeness but limited by model errors and uncertainties), such as the popularly used Global Land Data Assimilation System (GLDAS) [17]. The resolutions of hydrological models have become higher and more physical processes have been incorporated in the state-of-the-art models, making hydrological modelling more computationally demanding. More recent studies focused on improving the computation efficiency of hydrological modelling, such as those using the parallel computing strategy. Spatial domain decomposition is a parallel strategy that partitions a basin into a number of sub-basins and conducts hydrological simulations in different sub-basins among multiple processors [18]. Zhang and Zhou [19] developed a particle-set strategy to parallelize the flow-path network model to achieve higher performance in the simulation of flow dynamics. Compared to previous partition-based strategy, the strategy developed by Zhang and Zhou [19] focused on dynamic water flows instead of statistic basin units.

In addition to the topics introduced above, there are many other hydrological topics in which geo-spatial analysis plays a crucial role, such as water quantity and quality evaluation for sustainable management [20,21]. Given the increasing importance of geo-spatial analysis in hydrological studies, the Special Issue aimed to seek studies in a wide range of topics related to the development and application of geo-spatial analysis in hydrology.

#### **2. Content of the Special Issue**

This Special Issue was set up for exchange of the latest ideas, methods, and results of hydrological studies using geo-spatial analytic technique. The articles of the Special Issue cover a variety of topics in hydrology, including flood modelling, extraction of basin characteristics, and the monitoring of water quantity and quality for sustainable management using geo-spatial methods.

Floods are one of the costliest and deadliest natural hazards in the world. Flood simulation and visualization with acceptable accuracy and reliability are crucial for flood monitoring, forecast, and management. Effective and efficient policy and decision making for flood prevention and control require not only high-quality and fine-resolution hydrological data, but also timely processing and visualization of geo-spatial information. Geo-spatial technique is the key for flood simulation and decision support system to store, process, and visualize hydrological data and flood risks. Wu et al. (2019) integrated one- (1D) and two-dimensional (2D) hydrodynamic models with a spatio-temporal Geographic Information System (GIS) to dynamically simulate flood risks. In this decision support system, a three-dimensional (3D) model of the study area and hydraulic engineering facilities can be quickly established using the photography-based 3D modeling technology. Based on this framework, a multi-source spatio-temporal data platform for flood risk simulations was developed for the Xiashan Reservoir in the Weihe River as a case study. The model assessment results in Wu et al. [22] showed that the integration of spatio-temporal GIS and hydrodynamic models can improve the efficiency of flood risk simulations for decision support, such as dam-break flood simulations and dynamic visual simulations. Munir et al. [23] integrated a hydrological model and a hydraulic model based on remote sensing and GIS-derived estimates to simulate torrential streamflow response to flash flood events in Pakistan. The study made use of the integration of GIS and remote sensing technique to derive different hydrological parameters for the models at a pixel level. The integration of the models was found to be able to simulate flash flood conditions and extents more accurately. These two studies indicated that the integration of models with geo-spatial analytic methods can improve the efficiency and accuracy of flood simulations for flood management.

The second group of studies focused on using geo-spatial technique to improve the accuracy of extraction of basin characteristics. Li et al. [24] proposed a new method of watershed delineation for flat terrain based on Sentinel-2A images with the Canny algorithm on Google Earth Engine (GEE). Using the traditional DEM for water delineation in local-scale plains may not obtain realistic drainage networks primarily because of large depressions and subtle elevation differences. In their study, Sentinel-2A images were first used to identify water bodies such as rivers, lakes and reservoirs to compose the drainage network. Afterward, catchments were delineated based on the flow direction of these water bodies. The proposed method was applied in the Taihu Basin of the Yangtze River basin as an example. The catchment characteristics extracted from the Sentinel-2A images were compared with those based on DEM. Their results showed that the catchment delineation based on the Sentinel-2A images can more precisely represent drainage networks and catchments especially in the plain areas. In another paper of this Special Issue, a level of confidence approach was proposed to quantify the confidence and improve the accuracy of the bathymetry [25]. The quantification of the confidence of satellite-derived bathymetry is challenging because of the lack of in situ data for validation. This proposed approach considers multiple satellite-derived bathymetry techniques, i.e., empirical, classification, and photogrammetric (including automatic and manual). They found that the proposed approach increases the overall accuracy of the satellite-derived bathymetry. Furthermore, certain levels of uncertainties in bathymetry were removed in the proposed method.

The evaluation of water quality and crop water budgeting are two important issues in water resources management for achieving sustainable development of society and ecosystems. Remotely sensed data and geo-spatial methods have been widely used by researchers and policy makers to monitor water quality and estimate crop water deficit. However, the commonly used moderate resolution sensors hardly fulfill the monitoring requirements for small-sized water bodies. Avdan et al. [26] introduced the high-resolution RapidEye satellite data to assess water quality

parameters in a small dam, such as electrical conductivity, total dissolved soils, water turbidity, suspended particular matter, and chlorophyll-a. Their results showed that almost all of the water quality parameters have correlations with Rapid Eye reflection higher than 0.80. In Javed et al. [27], a remote sensing technique was used to estimate crop water deficit. Crop classification was determined by NDVI using crop cycles based on reflectance curves. The crop water deficit was defined as the difference between potential and actual evapotranspiration derived from the reflectance-based crop coefficients. Their results showed strong correlations between evapotranspiration, temperature, and rainfall. Crops in summer suffered a higher water deficit than in winter, because of the higher evapotranspiration demand due to higher temperature. Both studies demonstrated that remote sensing and geo-spatial methods are important tools for small-scale water quality assessments and agricultural water consumption evaluations for sustainable water management.

#### **3. Closing Remarks**

The articles in this Special Issue presented the latest results of geo-spatial analysis in different topics in hydrology, and developed new methods to improve the accuracy and reliability of the results of geo-spatial analysis. These studies showed that the integration of traditional hydrological/hydraulic models with geo-spatial techniques can improve the efficiency and performance of flood simulations for decision making. The needs of high-resolution hydrological data for scientific studies and practical operations at the local scale have been increasing. The results of these studies demonstrated the reliability of using geo-spatial techniques to acquire high-resolution hydrological parameters at a small spatial scale. In summary, geo-spatial analysis has been an essential and important element in a wide range of research topics in hydrology.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Dynamic 3D Simulation of Flood Risk Based on the Integration of Spatio-Temporal GIS and Hydrodynamic Models**

#### **Yongxing Wu 1,**† **, Fei Peng 1,**† **, Yang Peng 2, Xiaoyang Kong 1, Heming Liang <sup>1</sup> and Qi Li 1,\***


Received: 26 September 2019; Accepted: 13 November 2019; Published: 18 November 2019

**Abstract:** Dynamic visual simulation of flood risk is crucial for scientific and intelligent emergency management of flood disasters, in which data quality, availability, visualization, and interoperability are important. Here, a seamless integration of a spatio-temporal Geographic Information System (GIS) with one-dimensional (1D) and two-dimensional (2D) hydrodynamic models is achieved for data flow, calculation processes, operation flow, and system functions. Oblique photography-based three-dimensional (3D) modeling technology is used to quickly build a 3D model of the study area (including the hydraulic engineering facilities). A multisource spatio-temporal data platform for dynamically simulating flood risk was built based on the digital earth platform. Using the spatio-temporal computation framework, a dynamic visual simulation and decision support system for flood risk management was developed for the Xiashan Reservoir. The integration method proposed here was verified using flood simulation calculations, dynamic visual simulations, and downstream river channel and dam-break flood simulations. The results show that the proposed methods greatly improve the efficiency of flood risk simulation and decision support. The methods and system put forward in this study can be applied to flood risk simulations and practical management.

**Keywords:** spatio-temporal GIS; hydrodynamic model; spatio-temporal computation framework; flood risk; 3D simulation

#### **1. Introduction**

Water is important for human survival, but it is also the source of many disasters. According to the 2018 disaster report for China, flood disasters are one of the main natural disasters for the country [1]. It is crucial for ecologically-based development to understand how to scientifically prevent flood disasters and protect and properly use water resources. Hydrodynamic models are core quantitative calculations in emergency flood disaster management [2–4] and can accurately simulate the instantaneous dynamic evolution and medium- to long-term development of a flood [5]. Since hydrodynamic models are very complicated, usually only domain experts can understand the models. Therefore, three-dimensional (3D) visualization of the computation results is important in practice. In recent years, texture mapping technology, animation technology, and simulation technology have been used to build 3D virtual environments based on data such as LiDAR point clouds, satellite images, contour lines, and topographic maps. By using digital hydrodynamic models to calculate the water depth and flow velocity, flood evolution and inundated areas could be simulated and displayed in a 3D virtual environment, thus enhancing decision-makers', the public's, and other stakeholders' awareness of floods. It will also promote scientific management of flood disasters [6–9]. Although the effect of visualization (such as rainfall and rollback waves) in these studies is reasonably well presented [8], most of them use simulation software or animation technology to show the effects of floods and inundated areas at different times [9]. The flexibility of 3D presentations and spatial analysis functions should be strengthened. For example, queries of water depth at an arbitrary point and displays of flowrate and direction are not supported in a 3D virtual environment. On the one hand, 3D virtual environments have not been seamlessly integrated with hydrodynamic models, which read the computation results of hydrodynamic models (such as text files). On the other hand, terrain topology is established by using contours of topographic maps. This method takes longer to process, and its precision depends on the scale and production time of the map [7,8]. Along with oblique photogrammetry, the seamless integration of spatio-temporal GIS and hydrodynamic models provide a way to solve the problem.

Based on a unified spatio-temporal reference, spatio-temporal datasets comprise geographic elements or phenomena related to a location [10]. Spatio-temporal data have a spatial dimension (S), an attribute dimension (D), and a time dimension (T) [11]. Space–time data reflect quantitative and qualitative characteristics, the spatial structure and the spatial relationships among the geographical elements or phenomena, and their changes through time [12,13]. Spatio-temporal data are the basis for how humans understand the geographical world. A spatio-temporal Geographic Information System (GIS) is able to acquire, store, analyze, and visualize spatio-temporal data. It represents the position and spatial form of geographic objects or phenomena, and their changes over time. Compared with two-dimensional (2D) GIS and three-dimensional (3D) GIS, spatio-temporal GIS is more powerful for visualization and spatial analysis and can more accurately reflect how objects change [14–16]. Over time, an increasing amount of meteorological data have become available. The costs of high-resolution satellite images and oblique photogrammetric surveys have been reduced, and their reliability has improved. The new generation of information technologies (such as the Internet of Things (IoT), big data, and cloud computing) have rapidly developed and are increasingly applied in the water industry [17]. As a result, the number of data types has increased, and the volume of data has grown rapidly. These developments have enabled real-time simulations of flood risk and intelligent emergency management. In addition, there is an increased demand for data management, visualization, spatial analysis, and business integration. These developments have also led to changes from static 3D data to dynamic 3D data and time series data, from static to dynamic and continuous visualizations, from spatial analysis to real-time spatio-temporal simulations, and from decision support aids to operational running. Spatio-temporal GIS can better satisfy these changes, can better manage spatio-temporal data from flood disasters, and can be used to reveal patterns of spatio-temporal changes in incidents (i.e., floods) [18–21].

Hydrodynamic models are very complicated and involve large amounts of data [6–9]. The application value of a hydrodynamic model is affected by the following factors: (a) The availability, timeliness, and resolution of basic data pertaining to the river channel and its surroundings; (b) the effect of the visualization of the computation results; and (c) the extent to which the hydrodynamic models and the business system are integrated. The disaster-inducing factors of flood risk include the flood volume, the inundated area, the inundated depth, the inundation duration, the flood flow velocity, and the flow direction. Changes in the flood flow velocity, the flow, and the water level through time are important for decision making in flood disaster emergency management [22]. Spatio-temporal GIS methods allow three-dimensional visualizations of the locations and conditions of hydraulic engineering facilities, hydrologicalinformation, meteorological data, flood factors, and the results ofmodel computations. This type of GIS can visually and dynamically show spatio-temporal changes in floods and the spatial distribution of affected persons, facilities, and ecological environments [23,24]. It provides dynamic, quantitative, refined, and real-time information for decision making by flood evolution simulations, condition evaluation of hydraulic engineering facilities, disaster evaluation, and emergencymanagement. Additionally, it can provide powerful data and platform support for early disaster warning systems and monitoring and evaluation

of flood disasters [25–27]. Spatio-temporal GIS uses high-precision and high-resolution spatio-temporal data and dynamic visualization effects. It has greatly enhanced flood simulations, flood disaster forecasts, and early warning and emergency management of flood disasters. It provides a valid solution to data quality and visualization problems in hydrodynamic model computations [24]. The seamless integration of spatio-temporal GIS and hydrodynamic models will considerably solve the three problems that affect the application value of hydrodynamic models.

This study presents a technical method for quickly establishing a spatio-temporal GIS framework for a reservoir, river, and the surrounding areas using an oblique photogrammetric survey and the digital earth platform. This method allows the seamless integration of the spatio-temporal GIS framework with the hydrodynamic models. A platform for the dynamic simulation of flood risk was established; this platform can integrate data flow, business flow, computation resources, and visual decision making. The platform provides basic data, such as riverbed section locations and a regional digital elevation model (DEM), which are required for the hydrodynamic model computations. The resulting values, including the flow velocity, flow direction, water depth, and inundated range, are displayed in an integrated, dynamic, and three-dimensional way. The platform can also show the variations in the flood factors in three dimensions, relative to changes in the input parameters such as the rainfall and reservoir flood discharge. This platform provides quantitative, scientific, efficient, and visual decision support information for flood simulations, early flood disaster warning systems, and water resources management. A flood-dispatching dynamic and visual simulation platform for the Xiashan Reservoir and the Weihe River was established. The platform provides web-end functions, including data management, data query, model computation, flood simulation, and visual decision making. The methods described here are of great significance for the emergency management of flood disasters and for the intelligent dispatching of water resources in the era of IoT, big data, and cloud computing [17].

#### **2. Materials and Methods**

#### *2.1. Method for Rapidly Building a Spatio-temporal GIS Platform for Flood Risk Simulation*

#### 2.1.1. Data Acquisition Through Oblique Photography

A multi-lens oblique photography system was used to photograph the ground, based on the position of exposure points, to obtain multi-angle ground images with multiple degrees of overlap. With its enhanced performance and more convenient operations of graphics processing units (GPU), cloud computing, unmanned aerial vehicles, and digital cameras, oblique photography can rapidly generate a high-precision and high-resolution three-dimensional model of the physical world. This three-dimensional model can be resolved at the centimeter scale, which can truly and objectively represent the land surface configuration [28,29] (Table 1).


**Table 1.** Key quality requirements for high-definition image data and landform data.

Abbreviations: DOM, digital orthophoto map; DEM, digital elevation model.

The data from oblique photography were processed using smart3D. Photos captured from different angles were used as data sources for smart3D to read the information, such as photo locations and control points. Smart3d output 3D models of terrain and buildings with real textures without manual intervention, which could accurately show the geometrical morphology and detailed composition of ground objects. The processing flow is shown in Figure 1. A triangulation network (TIN) was established using a dense point cloud generated by aerial triangulation and the dense matching of images [30]. The TIN forms an untextured model that can reflect the 3D spatial form of objects. Calculating the corresponding texture from the images using the software and mapping the texture on the corresponding untextured model can form a real 3D scene that reflects the spatial relationships and surface features of objects. Then, data, such as digital orthophoto maps (DOMs), digital elevation models (DEMs), DSMs, 3D models (3DMs), DLG (Digital Line Graphic) and digital object models (DOBs), can be generated as needed.

**Figure 1.** Schematic diagram of the complete data processing flow for oblique photography-based three-dimensional (3D) modeling.

#### 2.1.2. Construction of 3D Models for Hydraulic Facilities

In this study, three-dimensional models of hydraulic facilities were constructed for separate queries, spatial analyses, connecting time series data of sensors, and so on. The model precision standards refer to those listed in Table 2.


To create the three-dimensional models, the contour lines of the hydraulic engineering facilities were imported into the modeling software. The outline structures of the reservoir dam and its appurtenance were obtained through a general survey. Editing of all models was completed in 3ds Max.

#### 2.1.3. Building of a 3D Scene for Flood Risk Simulation

Using LOD (level of detail) technology and by integrating spatial data at different scales, the real 3D scene was built to provide data, a computing platform, and visualization support for the dynamic simulation of flood risk in the spatio-temporal GIS framework [26,27] (the building process is shown in Figure 2). The large-scale data include a 1:10,000 topographic map and a 2.5-m resolution image of the study area and the surrounding areas within a certain range. The medium-scale data are the oblique photography-based 3D modeling data within the study area, including the DEM, DOM, and 3DM. The small-scale data include the 3D single-body models of the hydraulic engineering facilities and the 3D models of the equipment and facilities in buildings, such as the pump stations. The 3D scene generated from oblique photography (osgb format) was directly used in the 3D models of key areas. The textured 3D terrains were generated by superimposing the satellite image on DEM, which were used as the 3D models of the rest areas. These two 3D scenes and single 3D models of water conservancy facilities were imported into spatio-temporal GIS. The whole study area could then be displayed in a three-dimensional way through the spatio-temporal GIS.

**Figure 2.** Process of 3D scene construction for flood risk simulation.

#### *2.2. Hydrodynamic Models*

For simulating spatio-temporal variations in the river channel flow upstream and downstream of the reservoir and early warnings of flooding in the watershed, we selected the 1D and 2D hydrodynamic models to perform quantitative simulations of the flooding process in the downstream river channel, based on current and historical measured hydrological variables (including rainfall, discharge, and water level) [31,32]. The data required by the model and the output results are provided in Table 3.

**Table 3.** Relevant parameters for the hydrodynamic model system.


2.2.1. One Dimensional (1D) Unsteady Flow Model of the River Network

One-dimensional unsteady flow movement in a single river channel is often described using the Saint-Venant equations [31,33]—namely, the continuity equation of flow (Equation (1)) and the dynamic equation of flow (Equation (2))—as follows:

$$\frac{\partial \underline{Q}}{\partial \mathbf{x}} + B \frac{\partial \underline{Z}}{\partial \mathbf{t}} = q\_1 \tag{1}$$

*ISPRS Int. J. Geo-Inf.* **2019**, *8*, 520

$$\frac{\partial \mathcal{Q}}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (a\_1 \frac{\mathbf{Q}^2}{A}) + gA \frac{\partial \mathcal{Z}}{\partial \mathbf{x}} + g \frac{n^2 \mathcal{Q} |\mathbf{Q}|}{A R^{4/3}} = q\_1 (u\_1 - V) \tag{2}$$

where *t* is time (s); x is the flow path (m); *Q* is the discharge; *A* is the cross-sectional area of flow; *B* is river width; *Z* is the water level; R is the hydraulic radius; *n* is the Manning roughness coefficients; *V* is the average flow velocity of the section; *ql* and *ul* are the lateral inflow per unit length of the river reach and the component of the lateral inflow in the x direction; α<sup>1</sup> is the momentum correction factor, where α<sup>1</sup> = - *<sup>A</sup> <sup>u</sup>*2*dA* /(*Q*2/*A*), and g is the gravitational acceleration.

Flow and energy in single river reaches are exchanged at the bifurcation point. Therefore, movement of flow at the bifurcation point must conserve mass and energy. In other words, a balance between the net flux at the bifurcation point and the net change of the actual water volume at the bifurcation point (Equation (3)) must be maintained, as described in the equation as follows:

$$
\sum Q\_i = \frac{\partial \Omega}{\partial t} \tag{3}
$$

where *Qi* is the influx at the bifurcation point after passing through cross-section *I* and the influx at the bifurcation point (node) is positive, while the outflux at the bifurcation point (node) is negative. Ω is the water storage capacity of the bifurcation point. Considering the velocity head, resistance loss, and other factors at the end point of river reaches at the bifurcation point, the water levels at the end points of all river reaches at the bifurcation point shall meet the following requirements (Equation (4)):

$$\begin{aligned} Z\_i = Z\_j = \dots \dots = \overline{Z} \\ \Delta Z\_i = \Delta Z\_j = \dots \dots = \Delta \overline{Z} \end{aligned} \tag{4}$$

If the section of the river reaches connected by the bifurcation point is very close to the bifurcation point, resistance loss at the bifurcation point will be negligible. Therefore, it can be assumed that the water level at the end points of all river reaches at the bifurcation point is the same.

Equations (1) and (2) can be solved using numerical techniques [34]. The solution to these equations comprises estimates of *Q* and *Z* for every cross-section at each time step.

#### 2.2.2. Two-Dimensional (2D) Flow Routing Model of the Downstream River Channel

Generally, the 2D flow equation for the plane flood routing in the flood detention area can be described using the shallow water equation [35,36]. The continuity equation of flow (Equation (5)) is as follows:

$$\frac{\partial \mathcal{Z}}{\partial t} + \frac{1}{\mathcal{C}\_{\mathcal{E}} \mathcal{C}\_{\eta}} \frac{\partial}{\partial \mathcal{E}} (hv \mathcal{C}\_{\eta}) + \frac{1}{\mathcal{C}\_{\mathcal{E}} \mathcal{C}\_{\eta}} \frac{\partial}{\partial \eta} (hv \mathcal{C}\_{\mathcal{E}}) = 0 \tag{5}$$

The dynamic equations of flow (Equations (6) and (7)) are as follows:

$$\begin{bmatrix} \frac{\partial \boldsymbol{u}}{\partial t} + \frac{1}{\boldsymbol{C}\_{\xi} \boldsymbol{C}\_{\eta}} \Big[ \frac{\partial}{\partial \boldsymbol{\xi}} (\mathbf{C}\_{\eta} \boldsymbol{u}^{2}) + \frac{\partial}{\partial \eta} (\mathbf{C}\_{\xi} \boldsymbol{v} \boldsymbol{u}) + \boldsymbol{v} \boldsymbol{u} \frac{\partial \mathbf{C}\_{\xi}}{\partial \eta} - \boldsymbol{v}^{2} \frac{\partial \mathbf{C}\_{\eta}}{\partial \boldsymbol{\xi}} \Big] = -\boldsymbol{g}\_{\boldsymbol{C}\_{\xi}} \frac{1}{\partial \boldsymbol{\xi}} + \boldsymbol{f} \boldsymbol{v} \\\ -\frac{\boldsymbol{u} \cdot \nabla \boldsymbol{u}^{2} + \boldsymbol{v}^{2} \boldsymbol{n}^{2}}{h^{4/3}} \boldsymbol{g} + \frac{1}{\boldsymbol{C}\_{\xi} \boldsymbol{C}\_{\eta}} \Big[ \frac{\partial}{\partial \boldsymbol{\xi}} (\mathbf{C}\_{\eta} \boldsymbol{\sigma}\_{\xi \xi}) + \frac{\partial}{\partial \eta} (\mathbf{C}\_{\xi} \boldsymbol{\sigma}\_{\eta \xi}) + \boldsymbol{\sigma}\_{\xi \eta} \frac{\partial \mathbf{C}\_{\xi}}{\partial \eta} - \boldsymbol{\sigma}\_{\eta \eta} \frac{\partial \mathbf{C}\_{\eta}}{\partial \boldsymbol{\xi}} \Big] \end{bmatrix} \tag{6}$$

$$\begin{bmatrix} \frac{\partial \boldsymbol{\upsilon}}{\partial t} + \frac{1}{\mathbb{C}\_{\boldsymbol{\zeta}} \mathbb{C}\_{\boldsymbol{\eta}}} \Big[ \frac{\partial}{\partial \boldsymbol{\zeta}} (\mathbb{C}\_{\boldsymbol{\eta}} \boldsymbol{\upsilon} \boldsymbol{u}) + \frac{\partial}{\partial \boldsymbol{\eta}} (\mathbb{C}\_{\boldsymbol{\zeta}} \boldsymbol{\upsilon}^{2}) + \boldsymbol{\mu} \boldsymbol{\upsilon} \frac{\partial \mathbb{C}\_{\boldsymbol{\eta}}}{\partial \boldsymbol{\zeta}} - \boldsymbol{u}^{2} \frac{\partial \mathbb{C}\_{\boldsymbol{\zeta}}}{\partial \boldsymbol{\eta}} \Big] = -\boldsymbol{g} \frac{1}{\mathbb{C}\_{\boldsymbol{\eta}}} \frac{\partial \mathbb{Z}}{\partial \boldsymbol{\eta}} - \boldsymbol{f} \boldsymbol{u} \\\ -\frac{\boldsymbol{\upsilon} \sqrt{\boldsymbol{u}^{2} + \boldsymbol{\upsilon}^{2} \boldsymbol{u}^{2}}}{\boldsymbol{h}^{4/3}} \boldsymbol{g} + \frac{1}{\mathbb{C}\_{\boldsymbol{\zeta}} \mathbb{C}\_{\boldsymbol{\eta}}} \Big[ \frac{\partial}{\partial \boldsymbol{\zeta}} (\mathbb{C}\_{\boldsymbol{\eta}} \boldsymbol{\sigma}\_{\boldsymbol{\zeta} \boldsymbol{\eta}}) + \frac{\partial}{\partial \boldsymbol{\eta}} (\mathbb{C}\_{\boldsymbol{\zeta}} \boldsymbol{\sigma}\_{\boldsymbol{\eta} \boldsymbol{\eta}}) + \boldsymbol{\sigma}\_{\boldsymbol{\eta} \boldsymbol{\xi}} \frac{\partial \mathbb{C}\_{\boldsymbol{\eta}}}{\partial \boldsymbol{\xi}} - \boldsymbol{\sigma}\_{\boldsymbol{\xi} \boldsymbol{\xi}} \frac{\partial \mathbb{C}\_{\boldsymbol{\zeta}}}{\partial \boldsymbol{\eta}} \Big] \end{bmatrix} \tag{7}$$

where ξ and η represent two orthogonal curvilinear coordinates in the orthogonal curvilinear coordinate system; *u* and *v* represent the flow velocities in the ξ and η directions, respectively; *h* represents the water depth; *Z* represents the water level; *f* represents the Coriolis coefficient and can be calculated as *f* = Ω*<sup>d</sup>* sin φ where Ω*<sup>d</sup>* is the rotational angular velocity of the earth and φ is the latitude; *n* is the roughness coefficient; *vt* represents the coefficient of turbulent viscosity and is calculated using *Vt* = α*U* ∗ *h*, where α is a constant (α = 0.25 ∼ 1.0) and *U*\* is the frictional velocity; *C*<sup>ξ</sup> and *C*<sup>η</sup> represent the Lamé coefficients in the orthogonal curvilinear coordinate system, expressions of which are *C*<sup>ξ</sup> = *x*ξ <sup>2</sup> + *y*<sup>ξ</sup> <sup>2</sup> and *C*<sup>η</sup> = *x*η<sup>2</sup> + *y*η2, respectively; and σξξ, σξη, σηξ and σηη represent the turbulence stresses, which are calculated as follows (Equations (8)–(10)):

$$
\sigma\_{\xi\xi} = 2v\_l \left[ \frac{1}{\mathcal{C}\_{\xi}} \frac{\partial u}{\partial \xi} + \frac{v}{\mathcal{C}\_{\xi}\mathcal{C}\_{\eta}} \frac{\partial \mathcal{C}\_{\xi}}{\partial \eta} \right] \tag{8}
$$

$$\sigma\_{l\parallel\eta} = 2v\_l \left[ \frac{1}{\mathcal{C}\_{\eta}} \frac{\partial v}{\partial \eta} + \frac{u}{\mathcal{C}\_{\xi} \mathcal{C}\_{\eta}} \frac{\partial \mathcal{C}\_{\eta}}{\partial \xi} \right] \tag{9}$$

$$\sigma\_{\xi\eta} = \sigma\_{\eta\xi} = 2v\_l \left[ \frac{\mathbb{C}\_\eta}{\mathbb{C}\_\xi} \frac{\partial}{\partial \xi} (\frac{\upsilon}{\mathbb{C}\_\eta}) + \frac{\mathbb{C}\_\xi}{\mathbb{C}\_\eta} \frac{\partial}{\partial \eta} (\frac{u}{\mathbb{C}\_\xi}) \right] \tag{10}$$

Depending on numerical discretization strategies, the finite volume methods (FVM) were introduced in our study to obtain the numerical solutions [5,9]. To avoid the formation of a jagged velocity field and pressure field, the staggered mesh method was introduced. In addition, the SIMPLEC algorithm and under-relaxation technology were used to complete the correction of the velocity equation and the depth equation and further accelerate the convergence of the correction equation [5,37].

In this study, a method based on the empirical formulas [38,39] was introduced to estimate the Manning roughness coefficients, and the cross-sectional flow resistance was corrected [2,40].

#### *2.3. Spatio-temporal Computation Framework*

Floods are a common process on the Earth's surface and have apparent spatio-temporal characteristics. Multiple factors are often involved in spatial analyses of river water and dynamic simulations of flood risk. Spatio-temporal big data have laid an important foundation for flood monitoring, dynamic simulations, and risk assessments [36].

The spatio-temporal GIS framework is an important aspect of spatial information science. Hydrodynamic models are an important component of the study of the hydraulics. The seamless integration of spatio-temporal GIS and hydrodynamic models involves a crossover study of the two subjects [36]. Due to the complexity of hydrodynamic models, model researchers often emphasize studies on model parameters, model applicability, and model accuracy and develop separate software tools for model calculations. The preparation of data, the calculation process, and the output results of models create their own system, without in-depth docking with application scenarios. Most model outputs are 2D tables with poor visualization effects [2]. Spatio-temporal GIS has strong capabilities in data acquisition, management, storage, analysis, calculation, and multidimensional visualization. Calculations of hydrodynamic models are embedded into the spatio-temporal GIS, and seamless integration is performed on the data flow, calculations, outputs and storage, and visualization [15]. Spatio-temporal GIS provides DEM, rainfall, and other types of data inputs about the river landforms and flood detention areas that are used to calculate the hydrodynamic model, as well as calculation resources and storage of the model calculation results. Multidimensional dynamic visualizations are performed for the model results using the visualization capability of the spatio-temporal GIS, thus realizing the seamless integration of the spatio-temporal GIS with the hydrodynamic model, as shown in Figure 3.

**Figure 3.** Conceptual model of the spatio-temporal computation framework. GIS, Geographic Information System.

The logical framework for the seamless integration of the spatio-temporal GIS with the hydrodynamic model is shown in Figure 4. It is based on service-oriented architecture (SOA) and is divided into the hardware layer, the data layer, the platform layer, the service layer, and the application layer, from bottom to top. The data layer includes static 3D scene data, IoT sensed and socioeconomic data, and spatio-temporal flood risk data. The static 3D scene data mainly include DOMs acquired by oblique photography-based 3D modeling, 1:2000 DEMs, 3DM data, river body data, and 3D model data of hydraulic engineering facilities. The IoT perception data include time series data acquired by sensing equipment at the hydrometrical station and the precipitation station. The spatio-temporal flood risk data include the flow velocity, flow direction, discharge, water depth, and inundated area at different sections, locations, and times, which are the outputs of the hydrodynamic model. The core of the platform layer is a spatio-temporal GIS, which is based on a digital earth platform developed by the authors. The platform layer connects the preceding and the following procedures, performs data management for the following procedures, and provides service integration, calculations, and data for the preceding procedures [41,42]. The service layer encapsulates core modules or functions for services, which allow the application and operational system developers to focus on the business process and user demands and simplifies the development workload in the application layer. The service and application layers are loosely coupled to enhance system flexibility and scalability [43,44]. The application layer performs a dynamic visual simulation of the flood risk for decision support, acting as an interface for the conversion of data to information and knowledge.

**Figure 4.** Logical framework for the integration of spatio-temporal GIS with hydrodynamic models.

#### **3. Case Study**

#### *3.1. The Study Area*

As the largest reservoir in the Shandong Province, the Xiashan Reservoir is in the middle reach of the Weihe River (Figure 5). The total watershed area of the upstream and downstream reaches is up to 107.7 km2. It has a total capacity of 1.405 billion m3. The reservoir area is large and includes 4 counties and cities, 11 townships, and 97 immigrant villages [45]. The reservoir has been in service since September 1960. It is a large water conservancy project that integrates flood control, irrigation, power generation, aquaculture, urban and industrial water supply, and other comprehensive utilizations. The designed irrigation area of the reservoir is 1020 km2 and the effective irrigation area is 693 km2. The main dam of the Xiashan Reservoir is 2680 m long. The average annual rainfall in the reservoir area is 615.3 mm, approximately 80% of which is concentrated in the period from June to September. The study area includes the Xiashan Reservoir, 70 km from the upstream and downstream river channel of the Weihe River, and the flood detention area downstream of the reservoir, with a total area of 175 km2.

#### *3.2. Construction of a Spatio-Temporal Database*

Spatio-temporal data include basic spatial data, perception data, socioeconomic data, and spatiotemporal data of flood risk factors. The 3D spatial data of the reservoir area and the 70-km river channel upstream and downstream stretch of the Weihe River were acquired by oblique photography-based 3D modeling, including a DOM with a resolution of 0.1 m, a 3DM, and a 1:2000 DEM. The perception data are mainly the time series data at the precipitation station and the hydrometrical station. The demographic and socioeconomic data mainly include data on road traffic, demographic data from administrative regions at all levels, the locations of administrative regions, and the spatial positions and associated attributes of schools, hospitals, public security organizations, factories, and other institutes. The spatio-temporal data contents are shown in Table 4.

#### *3.3. Dynamic Visual Simulation System for Flood Risk*

By referring to Section 2.3 ("Spatio-temporal Computation Framework") (Figure 4) and using digital earth as the engine, a dynamic visual simulation system for the flood risk at the Xiashan Reservoir was established to complete the seamless integration of the spatio-temporal GIS with the hydrodynamic model. The system functions, as shown in Figure 6, mainly include such modules as data management, data query, model calculation, flood process simulation, and visual decision making. The spatio-temporal GIS module provides functions such as spatio-temporal data processing, visualization, spatio-temporal analysis, and spatial data management. The "Calculation of Hydrodynamic Model" module supports computing environment settings such as model selection, parameter input, and flood type selection. The "Simulation and Visualization of Hydrodynamic Model" module simulates the evolution process of floods and performs 3D visualization for the inundated area, water depth, flow velocity, flow direction, and so on. The "Query and Statistics" module can analyze and visualize the perception data from the precipitation and hydrometrical stations by time frame. The system can access different types of sensor data and other business system data through the Integration Interface. Based on the system, the hydrodynamic data can be input, calculated, and visually embedded into the flood process simulation and the flood control and emergency management services to improve the timeliness, scientific reliability, and visualization effect of the flood risk simulations [46].

**Figure 5.** Location of the Xiashan Reservoir: (**a**) The location of the Xiashan Reservoir in China; (**b**) the Xiashan Reservoir Dam; (**c**) the layout of the Xiashan Reservoir.

The system adopts a Brower/Server architecture and releases a GIS service as a Web Service [47]. It integrates the data from the rainfall monitoring system, hydrologic monitoring system, video monitoring system, meteorological system, land geological disaster system, and other operational systems into a unified spatio-temporal database. Under the unified spatio-temporal calculation framework, flood risk simulation, hydrological data query statistics, and other operations are packaged to complete the seamless integration of the data flow, operation flow, and system functions. The system architecture is shown in Figure 7.




**Table 4.** *Cont.*

**Figure 6.** Functional modules of the dynamic visual simulation system for the flood risk of the Xiashan Reservoir.

**Figure 7.** System functions and service architecture based on service-oriented architecture (SOA).

#### **4. Results**

#### *4.1. Three-Dimensional Visualization of the Study Area*

According to the methods in Section 2.1, a 3D visual simulation of the watershed, reservoir area, river, and hydraulic engineering facilities was realized by organizing and visualizing the spatial data in a spatio-temporal GIS framework. The 3D visualization of the study area is shown in Figure 8.

The system can quickly retrieve information about the ecological environmental system of the reservoir, river channel, hydraulic engineering facilities, and the peripheral area. Relying on the strong spatial analyses and visualization abilities of spatio-temporal GIS, the integration of the dynamic evaluation, visualization, and decision support for flood risk is greatly improved.

(**a**)

**Figure 8.** 3D visualization of the study area: (**a**) A bird's-eye view of the Xiashan reservoir and Weihe river; (**b**) A 3D view of the indoor facilities; (**c**) A 3D view of the administration agency.

(**b**) (**c**)

#### *4.2. Visual Simulation of 1D and 2D Flood Routing in the River Channel*

By using the 1D and 2D hydrodynamic models embedded in the system, values for the flooding in the river channels upstream and downstream of the Xiashan Reservoir were calculated, and the 3D routing simulation was eventually realized. The system automatically calculated the flood process in the river channels upstream and downstream of the Xiashan Reservoir according to the user-supplied input parameters and the system data. The calculation results include flood flow hydrographs of the channel upstream of the reservoir and its statistical flood characteristics (including time, water volume, and the highest water level of the river channel) (Figure 9).

**Figure 9.** Calculation results of the 2D flow mathematical model of the river channel downstream of the Xiashan Reservoir (river channel).

We dynamically showed the inundated area based on the flood and water level at any point on the river channel for a specific time period using time as the axis. A 3D visual simulation analysis of the evolution of the in-channel water level was produced by using the visualizing function of the system (Figure 10).

**Figure 10.** 3D simulation of two-dimensional (2D) flow in the downstream river channel. (**a**,**b**) The area of flood submergence and the color from red to blue represents the water depth from deep to shallow, respectively; (**c**) the dynamic monitoring and display of the water level in front of the reservoir dam during flood discharge; (**d**) the dynamic visualization simulation of one-dimensional (1D) and 2D hydrodynamic models in the river channels upstream and downstream of the reservoir, where the length of the arrows represents the flow velocity, the color of the arrows from red to green represents the water depth from deep to shallow, respectively, and the direction of the arrows represents the flow direction.

#### *4.3. Visual Flood Control Dispatching for the Xiashan Reservoir*

Based on the high-accuracy 3D scene data in the system and with time as the system axis, the model calculation results were dynamically simulated and shown, including the change in the water level in front of the reservoir dam over time for a given flood or runoff event. The degree of backwater in the reservoir area and the inundation level and its dynamic routing were simulated under the designated dispatching scheme (Figure 11).

**Figure 11.** 3D demonstration of the flood control dispatching process of the Xiashan Reservoir. (The arrow direction indicates the flow direction. The longer the arrow, the higher the flowrate. Color is used to show the depth of the water. From green, blue, yellow to red, the darker the color, the greater the water depth.)

This function module can perform reservoir flood control calculations based on a given (or forecasted) flood or runoff event and presents the changes in the reservoir water levels and discharged volumes. It also calculates the amount of backwater in the reservoir area, the inundation level, the discharged flow, the flow velocity, and the water level change under the dispatching scheme, in combination with the 1D hydrodynamic calculation model and the discharge calculation model of the spillway structure of the dam. The integrated calculation and dynamic visual simulation processes can improve the decision-making efficiency of flood control dispatching consultations and flood dispatching levels and minimize the impacts of flood risk [48].

#### *4.4. Visual Simulation of Dam Break*

Once the reservoir dam breaks, a destructive flood will form, causing severe casualties and loss of social properties on the downstream side, which would significantly influence social and economic development in the downstream area of the reservoir [32]. The dam-break flood analysis mainly involves calculating the flood process at the dam site and in the downstream area, including the flow and water-level hydrographs at the dam site and the discharge, water level, flow velocity, and peak arrival time in the downstream flood routing.

In addition, we also calculated the downstream flood fields at the moment of main dam break at T = 10 s, T = 0.5 h, T = 1.0 h, and T = 1.5 h (Figure 12). The figure shows that because the lateral length of the main dam is much greater than the width of the spillway, the flow velocity at the moment of main dam breakage (t = 10 s) is lower than it is at the moment of spillway breakage, with a maximum flow velocity of 16.77 m/s. The flood at the moment of main dam breakage covers a wider downstream inundated area compared with that at the moment of spillway breakage.

**Figure 12.** Calculation results of the main dam break.

Based on the DEM data from the middle and lower reaches of the system and on the 1D and 2D river channel hydrodynamic models embedded in the system, we dynamically simulated the downstream flood process of the reservoir after a dam break (Figure 13), including calculating the inundated area at different times and other flood characteristics, such as the water depth, flow velocity, and flow direction at any point.

**Figure 13.** 3D simulation of dam-break flow. (There are two rivers in this figure. After the dam broke, the area between the two rivers was also inundated.)

#### *4.5. Verification of Hydrodynamic Model and Sensitivity Analysis*

According to the water level and discharge data of seven typical available sections, we verified the 1D and 2D hydrodynamic models. Figure 14 shows a comparison between the calculated values of the water levels during floods at different design frequencies and the measured data from seven typical sections.

**Figure 14.** Preliminary verification of the mathematical model of flow in the river channel downstream of the Xiashan Reservoir.

As seen in Figure 14, under the designed discharge, the calculated value of the model is consistent with the measured data, indicating that the established mathematical model can largely simulate the dynamic pattern of water flow of the river channel downstream of the Xiashan Reservoir.

#### **5. Discussion**

Hydrodynamic models are fundamental to flood quantitative simulations. The accuracy and reliability of these computations are closely related to the availability, timeliness, and scale or resolution of data [5]. Hydrodynamic models are highly specialized; therefore, visualization and intuitive representations of their results are very important for scientific decision making. The integration of hydrodynamic model computations and emergency management is one of the important factors for improving the efficiency of flood emergency administration [4]. This paper offers methods for quickly establishing a spatio-temporal GIS framework for flood risk simulations and seamlessly integrating spatio-temporal GIS and hydrodynamic models. These methods can effectively solve the problems of data availability, visualization, and integration of hydrodynamic models, which will improve flood risk simulations and emergency administration.

The method proposed in this paper is used to quickly construct flood risk simulation spatio-temporal GIS, which provides up-to-date data for hydrodynamic model computations. Basic geo-information data, such as topography, are the basis of hydrodynamic model computations [5]. In practice, hydrodynamic model computations are often affected by the poor availability of basic data, poor timeliness, and low resolution of data. The following factors will lead to the failure of hydrodynamic model computations or the poor reliability and accuracy of computation results: (a) Difficulty in obtaining basic data for various reasons; (b) poor timeliness of basic data; and (c) low-resolution basic data. In this paper, tilt photogrammetry, three-dimensional modeling, and a digital earth platform are used to rapidly construct a spatio-temporal GIS framework, which provides up-to-date and high-resolution basic data for hydrodynamic model computations. This framework also includes increased storage and multi-dimension visualization.

The integrated framework realizes seamless integration of data, business information, and system functions contained within spatio-temporal GIS and the hydrodynamic model. It promotes the application of flood risk simulation results in emergency administration, which will improve the efficiency of decision making. In most realistic applications, the data acquisition system, the hydrodynamic model computation system, the visualization system, and the business management system are individual systems [18]. Therefore, interoperability among these systems is usually not very good, which affects the efficiency of emergency administration of flood disasters. The integration framework described in this paper integrates data collection, model computation, dynamic visualization, and emergency administration, which enables data flow and business processes to flow through a unified platform. In practice, this will shorten decision-making time in emergency management. Based on this framework, a dynamic visual simulation platform for the flood risk at Xiashan Reservoir was established, which could be used to conduct three-dimensional dynamic simulations of flood evolution The platform intuitively displays the variations in the inundation range over time and the changes in the flood velocity, flow direction, and flow discharge over time three-dimensionally, thus providing an easier way to understand flood evolution. The platform offers more scientific decision-support information for flood emergency management. Furthermore, the platform provides rich interfaces for the integration of monitoring systems, such as rainfall stations, hydrological stations, and so on. As a result, time series monitoring data can be connected to the platform in real time and can serve as input parameters for hydrodynamic models. In several cases, the static flood risk map will be created for specific rainfall and flood flow scenarios, such as a 100 mm rainfall event or a 50-year flood recurrence [49]. The map shows the inundation range, flood depth, etc. The methods proposed in this paper can be used to simulate flood processes under any rainfall or flow conditions, which can reflect changes in the flood inundation range as rainfall varies. Meanwhile, the simulation will be dynamically visualized in a three-dimensional way.

Flood risk simulations and emergency administration are complex systems, and their effectiveness and accuracy in realistic applications are influenced by many factors. Four research goals remain:

First, in common modeling practices, a physical model will be preferably constructed to assure the validity of the model outcomes [6]. This study focused on the construction of spatio-temporal GIS by using oblique photography and the seamless integration of spatio-temporal GIS and hydrodynamic models. Therefore, the physical models of Xiashan Reservoir and river channels were not constructed. The historical observation data from seven typical sections were used to verify the model (Section 4.5).

Second, the traditional hydrodynamic equations were used in the study. There are simplified equations in terms of the development of the hydrodynamic model [5]. The integration framework proposed in this study can integrate multiple hydrodynamic models. Therefore, hydrodynamic models themselves will not influence the integrated method proposed in this study. This paper focuses on the integration of three kinds of hydrodynamic models with spatio-temporal GIS. In the future, more hydrodynamic models will be integrated into the system and experiments will be conducted in other watersheds to improve the applicability of the method and the dynamic simulation system.

Third, from the view of model processes, the uncertainties include choice of model structures, model parameters, model inputs (e.g., channel geometry, floodplain), validation data, and land use change [5]. For this study, the main uncertainties were model structures, model inputs, and validation data. However, since the aim of the study was to propose a method for quickly establishing spatio-temporal GIS and put forward an integration framework, we did not quantify the uncertainties. The uncertainties will be identified, quantified, and represented in the future.

In addition, the platform's current infrastructure does not employ a big data architecture. Using the dynamic simulation of flood risk, population, economy, and society data will be added, so that the quantitative assessment of the affected population, structures, and economy under different flood scenarios can be analyzed [23,50]. Different emergency plans will be provided to promote the intelligent development of emergency management of flood risk. MongDB and Spark will be used to reconstruct the underlying data, and then the operation efficiency of the whole system will be further improved [51].

#### **6. Conclusions**

In this study, a spatio-temporal computing framework for flood risk simulation is proposed. A seamless integration of spatio-temporal GIS and hydrodynamic models is realized for the data flow process, computing process, business process, and system functions. This framework can perform dynamic and visual simulations of flood evolution and can identify flood risk in three dimensions, greatly improving the efficiency of flood risk simulations and decision making. The method described in this paper can quickly establish a flood simulation in a spatio-temporal GIS framework, which can provide up-to-date and high-resolution basic data to compute hydrodynamic models. Taking Xiashan Reservoir as an example, a visual decision support platform for the flood risk of the Xiashan Reservoir was established; this platform seamlessly integrated spatio-temporal GIS with one-dimensional and two-dimensional hydrodynamic models. Using the hydrodynamic models, the platform can simulate and display the dynamic flood evolution in three dimensions. The results show changes in the flood factors, such as the inundation range, the cross-section velocity and flow, the water depth, etc., over time. The platform can provide three-dimensional flood simulations of various scenarios for flood emergency administration. In addition, it provides more intuitive and scientific decision support information than the individual hydrodynamic model system provides.

A three-dimensional representation of the study area was quickly established using tilt photogrammetry. A spatio-temporal database was established and included spatial data (DEM, three-dimensional model of water conservancy facilities, submerged areas, towns, villages, etc.) and temporal data (temporal data of flood factors, meteorological temporal data, temporal data of sensing networks such as hydrological stations, etc.). Based on the digital earth platform, a spatio-temporal GIS of the study area was constructed to support model computation, spatial analysis and multi-dimensional visualization. A three-dimensional representation was conducted of the water conservancy engineering facilities in the research area, such as the reservoirs, river bodies, dams, pump rooms, etc.

Flood risk was simulated dynamically and in three dimensions. Through a seamless integration of spatio-temporal GIS and hydrodynamic models, a flood risk dynamic simulation platform was established. The platform integrated parameters input, computation processes, visualization of the computation results and emergency management decision making. In the three-dimensional representation, the temporal evolution of the key flood factors, such as the flood inundation range, the inundation depth, the flow velocity and the flow direction, was displayed. A dynamic and visual assessment of flood risk was conducted and the spatio-temporal characteristics of human factors, geographical objects, events, and flood processes were comprehensively displayed.

With the advent of the Internet of Things, big data, cloud computing, artificial intelligence, and 5G technology, the frequency, varieties, and volume of flood risk monitoring data will increase. Accurate flood risk simulations and intelligent emergency management will be the focus of future research. The integrated framework and methods proposed in this paper have laid a foundation for these types of future studies.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to express their gratitude to the project team.

**Author Contributions:** Conceptualization, Q.L.; methodology, Y.W. and Y.P.; software, Y.W. and F.P.; validation, Y.W.; formal analysis, Y.W. and F.P.; resources, Y.P. and Y.W.; data curation, H.L.; writing—original draft preparation, F.P.; writing—review and editing, Q.L., Y.W. and X.K.; visualization, Y.W. and F.P.; supervision, Y.P. and Q.L.; project administration, Y.P. and Y. W.

**Conflicts of Interest:** Yongxing Wu and Fei Peng contributed to the work equally and should be regarded as co-first authors. All authors declare no conflicts of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Integrated Hazard Modeling for Simulating Torrential Stream Response to Flash Flood Events**

**Bilal Ahmad Munir 1,\*, Sajid Rashid Ahmad <sup>1</sup> and Sidrah Hafeez <sup>2</sup>**


Received: 21 October 2019; Accepted: 10 December 2019; Published: 18 December 2019

**Abstract:** This study aims to monitor the flash flood response of Vidor/Wadore hill torrent in Pakistan by the integration of Personal Computer Storm Water Management Model PCSWMM (hydrologic) and HEC-RAS 5.x (hydraulic) models. The method leverages remote sensing and GIS derive estimates of measured and inferred parameters of Vidor rural catchment to quantify the flash flood events of the last four years: 2014–2017. The calibration of the PCSWMM is performed using the sensitivity-based radio tuning calibration (SRTC) tool. The Nash–Sutcliffe efficiency (NSE), coefficient of determination (R2), and relative error (RE) values were found between 0.75–0.97, 0.94–0.98, and −0.22–−0.09 respectively. The statistical indicators prove the accuracy of PCSWMM for rural catchments. The runoff response of Vidor torrent is also analyzed for 0.5/12.7, 1.5/38.1, and 2.0/50.8-inch/mm rainfall hyetographs. The generated hydrographs are used to simulate 2D-module in HEC-RAS 5.x for floodplain demarcation in the piedmont area. The accuracy of the flood extent is analyzed using spatial overlay analogy in the ArcGIS environment by comparing simulated and historically available flood extents. The simulated flood extent shows 76% accuracy with historic flood extent. The impact of flash flood events shows wheat, maize, and fruit orchards are the most effected agriculture in piedmont area. The results revealed that the integration of hydrological, hydraulic, and geospatial modeling approaches can be used to model a full picture of catchment response during flash flood events.

**Keywords:** flash flood; PCSWMM; curve number; rainfall-runoff model; HEC-RAS; Pakistan

#### **1. Introduction**

Flash floods are considered one of the most disastrous natural hazards because of their sudden and severe impact [1]. Intense rainfall on steep slopes of hill torrents engenders flash flooding with a short lag time that causes an unbearable economic loss. Human intrusion and variability in climate have modified the prevailing natural conditions and rainfall-runoff processes in mountainous catchments, ensuing an increase in the frequency of flash flood events. According to the statistical record of the National Civil Affairs (NCA), on average, a global loss of 3–6 % in the gross national product is caused by flash floods each year [2–8].

The rainfall-runoff behavior in steep mountainous catchment is a complex process solely depends upon the physical parameters (shape, size, type of stream, etc.) of catchment [9]. The accurate discharge measurement for such catchment remains a challenge due to the in-depth calculation of physical parameters. However, improvements and advancement in studies, on water balance and natural environment change, have propelled hydrological—as well as hydraulic—studies in past decades [10]. Different strategies have been developed by researchers [11–15] to integrate the effects of physical and human-induced factors on runoff for storm events. Among the different methods (Green Ampt, Hortons, Holton, etc.), Soil Conservation Service Curve Number (SCS-CN) has proven to be an enduring technique to quantify the flash flood discharge volumes in catchments [9,16].

In recent decades, hydrological models are proved as one of the effective measures used to predict and monitor flash floods. Globally, hydrometeorological approaches have been classified as the best management practices for flash floods. Hydrological models are classified into lumped, semi-distributed and distributed schemas. Lumped modeling incorporates uniform assumptive conditions; however, distributed models can incorporate distributed measured and inferential data (e.g., precipitation, solar radiation, temperature, soil moisture distributions, etc.). Distributed models perform calculations on both chorological as well as chronological scales. Semi distributed models perform in the same way as a distributed model with a difference of catchment scales. In semi-distributed model the smallest unit of a watershed is sub-catchment however, in the distributed model it further calculates the process on pixel (raster-based) scale. Distributed and semi-distributed models are equally implemented for flash flood monitoring, and have better performance than lumped models. Zoccatelli et al., Anquetin et al., and Jia, P et al., confirm that model results are significantly improved with distributed spatial forcing data (rainfall, temperature, soil moisture condition, etc.) [17–19]. However, hydrologic calculation remains a challenge due to the non-linear behavior of flash floods in spatiotemporal domain. Generally, different uncertainties in the source input data affect hydraulic models, particularly the available supporting data which determines the output of the model. Briefly, data limitations undermine the implementation of detailed physical-based models for the prediction of flash floods [20–25].

The monitoring of hydrological variables is complex in nature. Researchers make use of various computer-aided programs and models to interpolate and extrapolate diverse variables in spatiotemporal contexts [26]. Such platforms and models are valuable to explore and comprehend a system behavior as the models are used to identify errors and inconsistencies in the system. The advantage of different models is to improve user-defined scenarios. The developed setups portray the reality of a system with sufficient accuracy [26,27]. Researchers have simulated the flooding situation and perform flood risk zoning which provides technical support for flash flood control and disaster reduction in urban areas [28–31]. Personal Computer Storm Water Management Model (PCSWMM) is one of the models that offer hydrological (lumped, semi-distributed, and distributed) and hydraulic capabilities. However, it is mostly calibrated for urban catchments using gauge and flow data collected from the site [32–35]. Other models such as MIKE series (Denmark) and InfoWorks ICM (UK), with strong hydrodynamic ability and preprocessing, also have a good number of users. In comparison, the SWMM model (United States) is broadly used because its code is open-source and it is free to use [36]. SWMM has been extensively used for urban flooding, however, limited studies have shown it to be equally successful on rural catchments suing gauge data [7,37].

The frequency of flash floods in different mountain ranges of Pakistan has increased due to change in rainfall patterns in recent years [38–41]. The frequency of small to large range of flash floods is very high in D.G. Khan district. Vidor is one of the major hill torrents of the Suleiman mountainous range of D.G. Khan district. Frequent flash flood events in the Vidor catchment damage the infrastructure, valuable crops, and small villages in the piedmont area. The mighty flood of 2012 and 2015 also damage the irrigation canal systems in the downstream of piedmont plain at different RDs-reduced distances [7].

Based on historical records, and frequent flash flood events in D.G. Khan from the Vidor catchment, this paper validates the applicability of PCSWMM using semi-distributed modeling approach, and analyzes the impact of flash floods in downstream piedmont plains using HEC-RAS 2D hydraulic modeling. The accuracy of PCSWMM is cross-examined for rural catchments as it had been extensively used for urban catchments. The characteristics of the Vidor catchment are scrutinized using online data source in remote sensing and GIS domain. Finally, flash flood extents are demarcated for medium, high, and intense classes for vulnerability analysis using spatial overlay analogy. The proposed study will redound to the social benefits considering the importance of CCA (climate change adaptation) and DRR (disaster risk) reduction strategies.

#### **2. Materials and Methods**

#### *2.1. Study Area*

For this appraisal research, Vidor watershed/catchment in the administrative boundary of D.G. Khan district was selected to scrutinize the rainfall-runoff behavior for different curve number (CN) arrangements. Barren mountainous Vidor catchment with steep slopes receiving non-perennial flows during a storm is prone to disastrous flash floods [7]. The area was divided into two parts, Vidor watershed area, the source area of runoff generation, and the piedmont area of Vidor. The area has diverse terrain with elevation >200 m to <2500 m. The catchment consists of mostly hard rock strata with the barren land formation and dominating sandy soil [42]. The area lies in low annual average rainfall zones with rainfall schemes mostly adopting sudden surge from cloudbursts or thunderstorms and therefore categorizing the streams in the watershed area as episodic streams (Figure 1).

**Figure 1.** Study area map.

Geologically, the area was part of the structural trough which was occupied by an arm of old Tethys Sea continuously receiving the calcareous and argillaceous sediments. Quaternary and Tertiary deposits are the geologic formations of the area. Generally, sandy soil and alluvial/piedmont plain are present in the study area. The area consists of loam, silt loam, sandy loam, and sandy soils [43]. These soil types are considered as good for agriculture. However, sandy soils are an exception in the study area. The surface hydrology of the area is mainly controlled by torrential streams in the west. The streams pattern shows behavior in the west to east direction with an inclination in the south and tributaries in the catchment area shows a mixture of dendritic and contorted patterns. However, tributaries debouch from the piedmont line follows a distributary channel system in the piedmont plain. Stream network from main ordered stream dissipates into smaller tributaries scattered in all directions in the piedmont area. Most of the agriculture and urban areas are far away from the catchments that generate torrential flood. Sometimes, weather conditions recorded by the downstream weather station

in the city and the catchment areas are far different from each other. Mostly, these streams are activated suddenly without any weather symptoms or any warning sign [38].

The catchment area consists of shrubland in the west and diverse barren land with a mixture of rocky plains. The density of shrub cover is high in the western part; however, sparse vegetation cover is spread over the upper catchment area. The piedmont area contains seasonal patchy crop patterns including wheat, cotton, maize, and rice. In Pakistan crops are classified for two seasons namely: Rabbi (November–May) and Kharif (June–October). The exposed crop classes by flash floods are reclassified as single, double, and triple classes. The area with two and three crops in a year are classified as double and triple crops, respectively. The study area also contains double (wheat and cotton, wheat and maize, wheat and rice) and triple (rice) crop patterns in a calendar year.

#### *2.2. Models Description*

PCSWMM is an updated version of SWMM developed by the United States-Environmental Protection Agency (US-EPA). Many researchers have used SWMM/PCSWMM in the urban watershed for rainfall-runoff modeling [7,44]. However, ref. [10,45] also practice SWMM for hydrological modeling of rural watersheds and found it a suitable package for implementation over other rural watersheds. The model works on different blocks among which runoff block generates surface runoff governed by the equation of continuity.

PCSWMM runoff block's entire process for the generation of runoff is followed by the equation of continuity which tracks the volume of generated runoff within a sub-catchment. The runoff block deals each of the sub-catchment as a nonlinear rectangular reservoir as shown in Figure 2.

**Figure 2.** Concept of the sub-catchment as a non-linear rectangular reservoir in PCSWMM [36].

The rectangular sub-catchments only generate runoff when the water depth (d) exceeds the depression storage (dp) of sub-catchment. The system calculates the water depth (d) for the rectangular sub-catchment, nonlinear reservoir, by coupling the equation of continuity (1) with manning's Equation (2).

$$\frac{dv}{dt} = \frac{d(A.d)}{dt} = (A \times i\varepsilon) - Q \tag{1}$$

where *dv* = *d*(*A*.*d*) = change in volume store per unit time, *A*×*ie* = net inflow to the watershed, Q = runoff (outflow from the watershed)

$$Q = w \frac{1.49}{n} \left(d - d\_p\right)^{\frac{5}{3}} \times S^{\frac{1}{2}} \tag{2}$$

where *w* = sub-catchment width, *n* = manning's coefficient, *d* = water depth, *dp* = depression storage, *S* = average slope of the sub-catchment.

The coupling of Manning's and continuity equation into one nonlinear differential equation is used for the calculation of water depth (d). The equation was developed for calculating flows in open channels. However, the equation has been adapted for overland flow applications [45,46]. The nonlinear differential equation for sub-catchment behaving as a nonlinear reservoir is given below

$$\frac{\delta d}{\delta t} = i - w \frac{1.49}{n} (d - d\_p)^{\frac{5}{3}} \times S^{\frac{1}{2}} \tag{3}$$

The model solves the above nonlinear equation at each designed time step to quantify the response of each sub-catchment for a given hyetograph. Where <sup>δ</sup>*<sup>d</sup>* <sup>δ</sup>*<sup>t</sup>* = change in water depth, *i* = rainfall intensity.

The generated hydrograph from PCSWMM is used as a basic input for 2D hydraulic modeling through HEC-RAS version 5.0. The 2D system of the model discretizes the floodplain or river areas into a number of grid cells. The grid cells are attributed to friction (Manning's n) and elevation (DEM) values along the surface. The water surface elevation is developed using the topography of the individual grid cell, more commonly referred to as stage-storage curve analysis. The momentum equation (diffusion-wave form) is coupled with the continuity equation for the calculation of water surface elevation for individual grid cells at interval of time. The differential equation of mass conservation in the unsteady form given as

$$\frac{\delta D}{\delta t} + \frac{\delta (ud)}{\delta x} + \frac{\delta (vd)}{\delta y} + q = 0 \tag{4}$$

The model solves Equation (4) at each time interval for the individual grid cell. Where *t* = time; *D* = water surface elevation; *d* = water depth; *u* and *v* are the velocities in *x* and *y*-direction respectively; and *q* = source or sink term.

#### *2.3. Data and Processing*

The data used for hydrologic and 2D hydraulic modeling were collected from different online sources as well as from government organizations (Table 1).

For a particular storm event, the observed discharge data for Vidor torrent from the Irrigation Department is limited and was preprocessed to fill the data gaps. The cubical spline process is applied to the observed hydrographs to have continuous values at the 30 min interval. The time interval is set to match with simulated results of runoff. The process is a third-order polynomial equation. The higher-order polynomial equation is preferred over linear equations because runoff processes follow smooth trends. The general equation is given as

$$f\mathbf{j}\ \left(\mathbf{x}\right) = b\mathbf{0} + b\mathbf{1}\left(\mathbf{x} - \mathbf{x}\mathbf{0}\right) + b\mathbf{2}\left(\mathbf{x} - \mathbf{x}\mathbf{0}\right) + b\mathbf{y}\left(\mathbf{x} - \mathbf{x}\mathbf{0}\right)\left(\mathbf{x} - \mathbf{x}\mathbf{1}\right)\left(\mathbf{x} - \mathbf{x}\mathbf{2}\right)\tag{5}$$

where, x, x0, x1, x2 are the known values of variables; and b0, b1, b2, b3 are the unknown coefficient.


#### **Table 1.** Dataset used for modeling.

#### *2.4. PCSWMM Input Parameter Estimation*

There is only one station in the Vidor watershed for discharge measurement, and it has data gaps and limited records. Therefore, only peak events for the year 2014–2017 were used for analyzing the behavior of Vidor torrent. The maximum peak observed for each year is selected for simulating the PCSWMM.

The torrential floods are very fast and short events generated only from 3 h to 8 h rainfall. Hourly rainfall data was not available for the study area, therefore, satellite rainfall data from global precipitation measurement (GPM) was used for rainfall records due to its high temporal resolution (30 min).

Arc Hydro toolset of Arc-GIS platform was used for the delineation of Vidor watershed. The toolset uses DEM and flow direction as a major input for morphometric characteristic calculations of the watershed. For hydrodynamic modeling, the input parameters are classified as measured (surface elevation, catchment and channel geometry, node, etc.) and inferred parameters (infiltration parameters, depression storage, percentage impervious, roughness coefficient, etc.) [47]. The inferred parameters are calculated by the model and are used for the calibration of the model. Ebrahimian et al. use RS and GIS technique to estimate the sub-catchment area parameters for SWMM [48].

Land use land cover (LULC) of the Vidor catchment was classified through Landsat-8 satellite image. The image was processed under a maximum likelihood classifier (MLC) algorithm. The classified LULC and soil datasets were coupled to generate CN of the Vidor watershed. James et al. reports that in the watershed, an area with a slope greater than five percent should be adjusted for CN values [49]. Adjusted CN is calculated by multiplying the CN values obtained from the SCS general guide with a *K* factor.

$$\text{CN}\_{2a} = \text{CN}\_2 \times \text{K} \tag{6a}$$

where

$$K = \frac{322.79 + 15.63(a)}{a + 323.52} \tag{6b}$$

CN values, DEM (digital elevation model), and associated datasets were all in a raster format and were average out to small catchment level, therefore the slope factor was converted at such level by means of weighted slope factor.

$$Weight\ Sloop = \frac{\sum\_{i=1}^{n} \alpha\_i \times s\_i}{A} \tag{7}$$

where, α is the area (ha) of the slope, *si* is a slope in percentage, and *A* is the area of sub-catchment.

The calculation of input parameters is performed at the pixel level of 30 m and 12.5 m. The generated datasets are downscaled at sub-catchment levels for PCSWMM. The model has a unique database system that linked each sub-catchment with a reference identifier system. The complete system file is in .INP format.

#### *2.5. Calibration and Validation*

Vidor streams are non-perennial in nature and remains active for the monsoon season. Sudden response to short rainfall spells generates peaks greater than 0.1 million cusecs. Therefore, the PCSWMM was calibrated and validated for peak events only. The hydrologic model developed undergoes calibration process to adjust the input parameters and validation is performed using the real-time gauge data. A sensitivity-based radio tuning calibration (SRTC) tool was used for PCSWMM calibration. The generated model was calibrated for a peak event of 2014 and prediction was analyzed for 2015, 2016, and 2017 peak events. SRTC tool allows tuning the model for selected inferred parameters with a percentage of influence. However, the parameter selection is solely depending upon the performance of the model and uncertainties in the data. For the calibration width, Manning's (n) for the pervious area, drying time, and CN were used. The uncertainties for the calibration process based on [50–52] are shown in Table 2.

**Table 2.** Calibration parameter.


However, boundary conditions (BC) for 2D hydraulic simulations include: (1) BC for upstream is defined at outlet location of Vidor catchment (locally named as Darra Site) and (2) BC for downstream is defined at reduced distance RD-241 and RD-242 of DG Khan canal.

The calibration-validation results were evaluated using statistical criteria such as Nash–Sutcliffe efficiency (NSE), relative error (RE) (%), and coefficient of determination (R2). Whereas, the demarcated flood extents using a 2D hydraulic scheme was validated through overlay analogy of the ArcGIS platform using historic flood extent collected from Irrigation Department, Pakistan. Different researchers have discussed the Soil & Water Assessment Tool (SWAT) performance on the basis of NSE as: NSE > 0.65, 0.54 < NSE < 0.65, and 0.5 < NSE < 0.54 as very good, good, and satisfactory respectively [53–55]. However, ref. [56] have used RE for accuracy assessment of PCSWMM and found it a reasonable indicator for model performance assessment. The plus and minus sign of RE represents the underestimation and overestimation of the model respectively.

The methodological flow chart is shown in Figure 3.

#### *2.6. Flash Flood Exposure Assessment*

The flash flood exposure analysis is performed using spatial overlay analogy in the ArcGIS environment. It was preferred to compare the results with the field data from the Irrigation Department instead of the MODIS product due to the cloud cover in the imagery. The piedmont area is vulnerable to flash floods; however, local community manages the patchy crop production in the area. Small diversion structures are constructed in the area to overcome the impacts of flash flood events. LULC data acquired from the Punjab Irrigation Department was used to analyze the crop vulnerability in the piedmont plains. The 2D simulations result in different flash flood inundation scenarios.

**Figure 3.** Methodological flow chart.

#### **3. Results**

#### *3.1. Measured and Inferred Parameters*

Arc Hydro tools were used to derive the Vidor watershed morphometric characteristics. The generated hydrodynamic model consists of measured and inferred parameters. EI Alfy and Santhi et al. use Arc Hydro technique for delineation of catchment boundary and found it a suitable technique for deriving morphometric characteristics of watersheds [53,54]. Vidor watershed consists of a total of 88 sub-catchments with different characteristics. The ranges of measured and inferred parameters are shown in Table 3.


#### **Table 3.** Measured and inferred parameters.

#### *3.2. PCSWMM Simulations*

The integration of measured and inferred parameters at the sub-catchment level allows each catchment to behave independently for rainfall hyetograph. The system of sub-catchments generates a semi-distributed scheme as shown in Figure 4. Each sub-catchment behaves as the smallest unit of the system. Every catchment has a dual behavior in the designed system for a particular event including (1) the runoff generated by the individual sub-catchment, and (2) the contributed discharge from the connected upstream catchment known as run-on. The individual sub-catchment generates runoff for a particular rainfall event; however, due to upstream and downstream catchment connections the total runoff volume and peak is a combination of the runoff and the run-on. The discharge variations for each sub-catchment was observed for the peak events of 2014–2017 with a constant time interval Δt = 30 min as shown in Figure 4.

**Figure 4.** Sub-catchment behavior observed at arbitrary time t0 (**a**), t1 + Δt (**b**), t2 + Δt (**c**), and t3 + Δt (**d**).

The calibration and validation results show high correlation with the observed discharge measurements. Figure 5 shows the calibrated and validated scatterplots on event basis.

It can be seen from Table 4 that model performed efficiently for Vidor torrent when tested on peak events on the statistical criteria of NSE, RE, and R2. The model dropped NSE and R2 value for the year 2016 and remains in an excellent class for calibrated and validated periods. It has been observed that the model for whole events overestimates the results with fractions of difference (−0.22 < RE < −0.09).

**Figure 5.** Scatterplots for observed and simulated results.


**Table 4.** Hydrologic modeling results.

Flash flood event of the year 2012 was recorded as the maximum peak of Vidor catchment [7]. However, due to the lack of high accuracy data of rainfall and discharge historical events were not simulated. GPM data was available since 2014 at 30 min resolution and has been used exorbitantly all over the globe for flood studies [57–59]. GPM satellite data found suitable for storm burst event simulations for Vidor torrent. It was observed that GPM generates high accuracy (NSE, RE, and R2) for high-intensity rainfall events and accuracy reduces for low-intensity peak events (Table 4).

The dry and abrupt nature in the topography of Vidor torrent makes it a unique catchment. The runoff is entirely depending upon rainfall behavior in the catchment with no base flow. The abrupt nature of rainfall on Vidor torrent results in sudden response due to steep slopes. The result substantiates no influence of evapotranspiration on runoff due to very short time of flash flood events. The data interpolation of observed discharge and the rainfall-runoff direct relation results to achieve such high values of statistical indicators.

#### *3.3. Data Limitations and Interpolation*

The observed data record (2014–2017) has a limited number of observations and was interpolated using cubical spline statistics. The data was completed for the peak events and was used for accuracy assessments. The limited number of observed records are shown in Figure 6. Different researchers have used the cubical spline methodology for the calculation of missing datasets [60–66]. For uneven distribution of gap, it is preferred to apply piece-wise cubical spline interpolation. In this technique, each gap is interpolated individually. However, the size of the gap between two observations affects the accuracy of interpolation.

**Figure 6.** Observed records and simulated hydrograph.

High variation in the slope of Vidor torrent substantiates random and sudden response depending upon the shape of hyetograph. The intensity and duration of rainfall define the lag-time (time between the peaks of hyetograph and hydrograph) of torrent response to flash flood events. The Table 4 summarizes the storm duration, total rainfall, and resulting lag-time. The PCSWMM results show an excellent agreement between simulated and observed flash flood peaks for four years' data.

It was found that the slope adjusted CN scheme performs efficiently for event-based simulation [7,16,49,67]. However, soil physical infiltration models (Horton, Green Ampt, etc.) behave more efficiently for monthly and annual quantification of runoff [68]. The study demonstrates the supremacy of modern GIS and RS coupling techniques in the integration of different parameters at a pixel level. The slope shows a huge impact on the modification of CN and hence the surface runoff. [69–73] utilizes GIS and RS techniques for accurate measurement of CN values.

#### *3.4. Flood Extent Analysis (2D Simulations)*

The flood plains are demarcated with recursive 2D hydraulic simulations. The designed storms of 0.5, 1.0, and 1.5-inch intensity of SCS type at 6 h duration were generated in PCSWMM. The simulated hydrographs were further used for flood plain zoning in the down piedmont plains. Furthermore, changes in the depth and flow area were scrutinized for designed storms.

It was observed that the behavior of flash flood water in the piedmont plain is solely dependent on the shape of hyetograph, time of peak, lag time, roughness coefficient (soil × LULC), and DEM (elevation and slope). However, the main flow path remains the same with the increase of inundation depth. The flow in piedmont plain lacked in channel flow due to flash flood events and its non-perennial nature.

The behavior of flash flood water for 0.5, 1.0, and 1.5-inch SCS 6-h storm hyetograph in comparison with flash flood event 2015 is shown in Figure 7. The results for different storm intensities show the change in expansion of flood zone, branched from natural path of the nullah. With the increase in the intensity of storm, width of the flood plain is also increased with widespread, generating a fan-shaped zone in the piedmont area. Contrarily, low-intensity flood events also follow the main path of the natural nullah but the width of the flood plain decreases sufficiently, generating a small fan-shaped zone in the piedmont area. The depth of the flash flood water remains high in the existing nullah, but the depth diminished at the outskirts of nullah due to the steep topography of the piedmont plain (upstream to downstream) and the vague footprints of natural nullah. It was observed that for all developed scenarios the peak flood depth appeared near the outlet of the catchment where the depth of the channel is high. However, in the piedmont plain, the flood water depth decreases sufficiently due to the dispersion of flood in different directions.

The 2D hydraulic modeling procedures in HEC-RAS 5.x not only provide better and accurate (DEM resolution dependent) floodplain demarcation as compared to 1D hydraulic modeling but the methodology is also time and effort efficient in data-limited areas. However, the 1D hydraulic modeling model needs extensive data inputs including a storage area, cross-sections, weir structures that require a lot of time and effort [34].

In the current study, the multiple 2D hydraulic simulations result in the classification of flood extents for medium, high, and intense flash flood scenarios. The classification is based on the peak event discharges observed in the past years and the areal destruction recorded (local-knowledge) during the event. The discharge magnitudes of 0.05–0.1 million cusec, 0.1–0.15 million cusec, and 0.15–0.2 million cusec are classified as medium, high, and intense flood respectively. The medium, high, and intense flood covers 19557.58, 26875.28, and 29087.97 ha of area respectively as shown in Figure 8. Whereas historically available flood extent, extracted from MODIS data, for the year 2012 during a high flood event is 20519 ha.

#### *3.5. Flash Flood Exposure and LULC Vulnerability*

Simulated flood extent boundary, considering intense flood scenario, matched 76% with historically available flood extent boundary. The main flow area remains the same when compared with the available flash flood extent data collected from the Irrigation department. Ahmed et al. use multi-temporal Landsat-8 and MODIS satellite products to analyze the impact of the flash flood on rice crops [74]. However, ref. [75] uses Landsat 7 to delineate maximum flood extent in North Carolina. Tanguy et al. use high-resolution Radarset-2 dataset for flood monitoring [76].

**Figure 7.** Inundation behavior for 0.5/12.7, 1.5/ 38.1, and 2.0/50.8 inch/mm and flood 2015.

**Figure 8.** Classified floodplain map for Vidor torrent.

It was observed that barren land, rural settlements, and wheat classes cover most of the floodplain substantiated by intense flash flood events as shown in Table 5 and Figure 9.


**Table 5.** Vulnerable LULC for intense flood simulated class.

**Figure 9.** Floodplain LULC and exposed percentage.

#### **4. Discussion**

Pakistan is a country with varied topography—from alpine forests in the north to the barren mountainous catchment in the south and west regions. Mostly, low-lying areas including piedmont plains on the downstream of mountainous catchment face low to high-intensity flash flood events and during last few years, Pakistan has experienced an increased number of flash flood events owing to the changes in weather pattern [38,39]. Flash floods in Pakistan are mostly caused by monsoonal torrential rains, and a shift of monsoonal patterns has increased the frequency of flash flooding [30]. Hill torrents of southwest Punjab (Pakistan) are extremely vulnerable to flash floods. Historically, flash floods of small to large intensity frequently damage agriculture in the west of Dera Ghazi (D.G.) Khan District. The main source of flash flooding in these areas are the hill torrents situated nearby the area including Vidor, Chachar, Mithawan, Sakhi Sarwar, Kaha, Sanghar, and Sori Lund. Hill torrents of southern Punjab, after receiving torrential rains during monsoon season, frequently generate flash floods in the area. In September 2012, heavy spells of Monsoon rainfall on Suleiman hills (Vidor catchment) of D.G. Khan cause disastrous conditions in the city area and adjoining rural areas [7].

In Pakistan, flash flood management measures are done structurally by embankments, studs, spurs, and various flood protection techniques. For the safety from small as well as large-scale destruction flood walls with diversion structures are constructed. Structural measures delay action dams and channelization of flood water through small streams are also some techniques that are used for flood protection at the local level. Mirani Dam constructed in 2006 is an example of a delay action

dam in Baluchistan for the storage of hill torrent water for the sake of irrigation purposes and for the protection of the low-lying area from any flash event [38].

The study tries to combine multi-steps methodology for alternative cost and time effective parsimonious methodologies used for flash flood water assessments in rural areas. The main goal of the research was to use a semi-distributed modeling approach to calibrate the PCSWMM for rural Vidor' catchment prone to flash flood. Satellite rainfall, GPM, data with a temporal resolution of 30 min accurately predicts the peak discharges for flash flood events of 2014–2017. The reliable results of statistical tests prove the applicability of GPM data for torrential catchment in a semi-distributed modeling approach. Hydrologic modeling requires extensive data records for processing when performed for flash flood events. The input parameters were generated using remote sensing and GIS techniques. Different researches have shown the accuracy of semi-distributed and distributed hydrologic modeling for flash flood prediction and monitoring. The models are equally applicable for catchments with steep slopes and high topographic relief [17–19]. The generated slope adjusted CN scheme and designed methodology for the semi-distributed model generate consistent results with less data requirement and is replicable in other torrential catchments with ease. The CN scheme with different antecedent moisture conditions results in high accuracy of peak discharges. Rozalis et al., also used CN schemes on Merhavia watershed in northeastern Israel and found it a reliable approach for flash flood prediction [77]. However, ref. [78] used GIS-based model to evaluate the flood discharge of Xirolaki torrent in Northern Greece. Contrarily, the approaches for flash flood inundation modeling requires extensive data measurements. The results substantiate that for foothill areas, piedmont plains with vague cross-sections, 2D hydraulic modeling without incorporating cross-section data in piedmont plains can be replaced with 1D hydraulic modeling. However, 2D hydraulic simulation solely depends upon the accuracy and refinement of elevation data, therefore preprocessing including Fill sinks, depression evaluation, and DEM reconditioning must be performed prior to simulate 2D hydraulic model. The processed elevation data incorporates in 2D hydraulic model substantiates sufficient accuracy for flash flood extents. The integration of hydrological modeling (PCSWMM) with open source data and 2D hydraulic modeling (HEC-RAS 5.x) links the research for flash flood hazard modeling in data-poor areas.

This study also highlights the value of open-source information for flash flood hazard assessment. The freely available ALOS PALSAR elevation data, 12.5 m spatial resolution, generates a good prediction of flood extents in the piedmont plains. The multi-step scheme for flash flood assessment is simple and effective in nature for areas with limited data and hence could be implemented on other mountainous catchments with similar topography.

#### **5. Conclusions**

The study utilizes an integrated procedure for the quantification of flash flood events and demarcation of floodplain at foothill area of Vidor torrent. Hydrological investigations were performed using PCSWMM suit to scrutinize the rainfall-runoff relation. The generated hydrographs were tested for piedmont plain to analyze the simulated flood extents (medium, high, and intense) and depth distributions using 2D hydraulic modeling domain of HEC-RAS 5.x. The foregoing research tries to attempt the disaster risk reduction strategy in D.G. Khan for Vidor torrent. The integrated approach of hydrological, 2D-hydraulic modeling, and GIS-RS approach is used to identify the potential disaster of Vidor torrent. The study quantifies the flash flood extent for designed and peak events observed in the past years 2014–2017. The target area covers the Vidor catchment (source of flash flood) and the piedmont area in the down-plains of the source catchment. The rainfall-runoff simulation and flood extent demarcation with sufficient accuracy unveil the flash flood problem in the area. The demarcated flood extents are classified in medium to intense flood scenarios substantiates the basics for an early warning system in the area. The study concludes the following testimonials:


#### **6. Recommendations**

The changing behavior of climate is expected to affect flooding through rainfall disturbing patterns. The abnormalities in rainfall exacerbate the existing effects of flash flooding on community services, infrastructures, etc. in Pakistan. The priorities, therefore, of flood risk management should be changed to adapt the changing behavior of climate. Vidor is an active non-perennial torrent and frequently affects the settlements and crops in the piedmont area. The management of present-day and future risk from Vidor flash floods should be a combination of possible risk mitigation through structural or regulatory measures. In a study by [7], they proposed two small dams with a storage capacity of 13.31 and 14.29 million cubic meter capacity. The proposed structural measures by [7] are a part of disaster risk management and will lock the communities in the future of no risk from flash flooding.

**Author Contributions:** Conceptualization: Bilal Ahmad Munir, Sajid Rashid Ahmad, Sidrah Hafeez; Analysis: Bilal Ahmad Munir, Sidrah Hafeez; Writing—Original Draft Preparation: Bilal Ahmad Munir, Sidrah Hafeez; Writing—Review and Editing: Bilal Ahmad Munir, Sajid Rashid Ahmad, Sidrah Hafeez. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to acknowledge the Punjab Irrigation Department and Pakistan Meteorological Department for providing rainfall, discharge, historical flood extents, and LULC of piedmont plains datasets.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Van Westen, C.J.; Alkema, D.; Damen, M.C.J.; Kerle, N.; Kingma, N.C. *Multi-Hazard Risk Assessment: Distance Education Course Guide Book*; United Nations University—ITC School on Disaster Geoinformation Management: Tokyo, Japan, 2009.


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Consideration of Level of Confidence within Multi-Approach Satellite-Derived Bathymetry**

#### **René Chénier \*, Ryan Ahola, Mesha Sagram , Marc-André Faucher and Yask Shelat**

Canadian Hydrographic Service, Fisheries and Oceans Canada, 200 Kent Street, Ottawa, ON K1A 0E6, Canada; Ryan.Ahola@dfo-mpo.gc.ca (R.A.); Mesha.Sagram@dfo-mpo.gc.ca (M.S.); Marc-Andre.Faucher@dfo-mpo.gc.ca (M.-A.F.); Yask.Shelat@dfo-mpo.gc.ca (Y.S.)

**\*** Correspondence: Rene.Chenier@dfo-mpo.gc.ca; Tel.: +1-613-220-5026

Received: 4 December 2018; Accepted: 16 January 2019; Published: 19 January 2019

**Abstract:** The Canadian Hydrographic Service (CHS) publishes nautical charts covering all Canadian waters. Through projects with the Canadian Space Agency, CHS has been investigating remote sensing techniques to support hydrographic applications. One challenge CHS has encountered relates to quantifying its confidence in remote sensing products. This is particularly challenging with Satellite-Derived Bathymetry (SDB) where minimal in situ data may be present for validation. This paper proposes a level of confidence approach where a minimum number of SDB techniques are required to agree within a defined level to allow SDB estimates to be retained. The approach was applied to a Canadian Arctic site, incorporating four techniques: empirical, classification and photogrammetric (automatic and manual). Based on International Hydrographic Organization (IHO) guidelines, each individual approach provided results meeting the CATegory of Zones Of Confidence (CATZOC) level C requirement. By applying the level of confidence approach, where technique combinations agreed within 1 m (e.g., all agree, three agree, two agree) large portions of the extracted bathymetry could now meet the CATZOC A2/B requirement. Areas where at least three approaches agreed have an accuracy of 1.2 m and represent 81% of the total surface. The proposed technique not only increases overall accuracy but also removes some of the uncertainty associated with SDB, particularly for locations where in situ validation data is not available. This approach could provide an option for hydrographic offices to increase their confidence in SDB, potentially allowing for increased SDB use within hydrographic products.

**Keywords:** Canadian Hydrographic Service; Satellite-Derived Bathymetry; empirical; classification; photogrammetry; level of confidence

#### **1. Introduction**

The Canadian Hydrographic Service (CHS) is responsible for providing hydrographic products and services to ensure safe, sustainable and navigable use of Canada's waterways. As Canada contains the longest coastline in the world, CHS, with support from the Canadian Space Agency through a Government Related Initiatives Program (GRIP) project, has been exploring remote sensing technologies to help improve its nautical products. CHS's GRIP focuses on specific applications of Earth Observation (EO) data: shoreline extraction, change detection and Satellite-Derived Bathymetry (SDB).

To date, SDB research has focused on the development of novel empirical [1–3] and physics-based approaches [4,5], as well as the application of these techniques [6–8]. CHS investigations of empirical [9] and photogrammetric [10] SDB techniques have illustrated the potential of these approaches for deriving accurate bathymetry estimates in Canadian waters. While representing excellent progress for estimating water depth from satellite imagery, significant challenges remain with providing hydrographic offices (HOs) with sufficient confidence in SDB estimates to allow the information to be incorporated into official nautical products. Internationally, all sources of water depth information must meet accuracy specifications for specific depth ranges, defined by International Hydrographic Organization (IHO) CATegory of Zones Of Confidence (CATZOC) levels (refer to Section 4.3 for a detailed overview of CATZOC levels) [11]. As CHS demonstrated in [9] and [10], SDB obtained through various techniques can regularly achieve CATZOC level C. While this represents a reasonable accuracy for satellite-derived information, there is a desire by HOs to improve their confidence in SDB results. As well, CATZOC determinations can only be obtained when in situ bathymetric information is available, limiting HO confidence in SDB results when only limited or no in situ depth information is available.

This paper proposes a new technique for quantifying confidence in SDB estimates: a level of confidence assessment. This method operates by determining the level of agreement between multiple SDB techniques applied to a single site. Initially, agreement is assessed between all of the applied techniques (four for this study: empirical, classification and photogrammetric (automatic and manual)). Subsequent agreement is then evaluated for different combinations of techniques (e.g., three agree, two agree), allowing SDB coverage to be gradually increased while maintaining the highest possible degree of confidence in the results. Final SDB estimates are obtained via averages of the SDB results for the techniques incorporated into each combination.

The results of this work demonstrate the benefit of completing a level of confidence assessment and suggest a new approach for HOs to evaluate SDB results:


CHS believes that through the completion of level of confidence assessments, international HOs can increase their overall confidence in SDB, potentially allowing for increased use of SDB within international hydrographic products.

#### **2. Materials and Methods**

#### *2.1. Study Site*

CHS's primary use of SDB will be in the Canadian Arctic where the greatest concentration of gaps is present in CHS's hydrographic surveys. The site for this study is situated in the waters near Cambridge Bay (69◦07'N, 105◦02'W), a hamlet located on Victoria Island, Nunavut (Figure 1). Water in Cambridge Bay is generally clear with a depth visibility of around 15 m. The bottom is mostly composed of sand and rock but the benthic environment is more heterogeneous with numerous patches of vegetation, making the site somewhat complex for SDB.

**Figure 1.** Location of the Cambridge Bay study site on Victoria Island, Nunavut.

#### *2.2. Satellite Imagery*

A WorldView-2 stereo pair acquired on 20 September 2015 over Cambridge Bay was used for this study (Figure 2). The forward image was used for the empirical and classification SDB approaches. Table 1 provides details of the viewing geometry for the image pair.

**Figure 2.** WorldView-2 stereo pair used for each Satellite Derived Bathymetry (SDB) technique. (**a**) Forward image, used for empirical and classification approaches. (**b**) Backward image, only used for photogrammetric techniques. Imagery © 2015, DigitalGlobe, Inc.



#### *2.3. Hydrographic Surveys*

For this research, in situ hydrographic information was required for two aspects of the study: training the models used for the empirical and classification approaches and for assessing the accuracy of all techniques. The in situ dataset is composed of five CHS hydrographic surveys: three multibeam sonar surveys acquired in 2014, 2015 and 2017, as well as two Light Detection and Ranging (LiDAR) surveys acquired in 1985 and 1992 using the Larsen 500 system (Table 2). The age of the LiDAR datasets is less than ideal and results in some uncertainties relative to the multibeam datasets (refer to a detailed discussion in [10]). Nevertheless, the LiDAR surveys provide an important source of validation information for shallow depths, particularly from 0–2 m. The geographic coverage of the LiDAR measurements is also more widespread than the multibeam surveys (Figure 3), making it critical for understanding spatial patterns of uncertainty within the SDB estimates.


**Figure 3.** Spatial distribution of survey data up to 20 m in depth within Cambridge Bay. Note the wide coverage of the LiDAR survey relative to the multibeam datasets. Imagery © 2015, DigitalGlobe, Inc.

#### *2.4. SDB Methods*

SDB estimates were derived independently from the WorldView-2 imagery through four techniques as described below. The results from these techniques were then used to develop final SDB estimates through the level of confidence approach.

#### 2.4.1. Empirical

Empirical SDB methods develop a relationship between water depth (most commonly through in situ depth measurements) and optical band values to allow depth to be estimated across an image. For this study, the multi-band approach developed by Lyzenga [2] was selected for use as it has consistently achieved good accuracy within previous CHS investigations [9,12]. This method first identifies water-leaving radiance (*Lw*) values for optically deep waters (*L*∞). The assumption is that any pixels containing values above *L*∞ have contributions from the reflectance of the underwater surface, which is exponentially attenuated by the water column. Therefore, constants can be empirically determined to link the log-corrected difference between *Lw* and *L*<sup>∞</sup> of multiple bands to the depth, as shown in Equation (1) [2]:

$$z = a\_0 + \sum\_{i=1}^{N} a\_i \ln[L\_{\text{wi}} - L\_{\text{coi}}] \tag{1}$$

where *z* is depth, *i* indicates band specific parameters and *a* represents the empirically determined constants. Due to its characteristics, this method is strongly dependent on the input water depths used for training, as well as the characteristics of the imagery (e.g., degree of visibility to the bottom, sun glint presence, and atmospheric characteristics). The survey data used to supply training water depths represented the best available for the site, while the input WorldView imagery contained good bottom visibility, limited atmospheric influence and minimal sun glint.

The multi-band model assumes that the grey values of each band contain a linear relationship with water depth. To determine which bands were best to incorporate into the model for each study site, scatter plots of band values relative to survey depths were generated, allowing the degree of linearity of the relationship to be evaluated. For this study, WorldView-2's blue, green, red and yellow multispectral bands were used.

#### 2.4.2. Classification

A random forest decision tree classification [13] technique was also applied. As a first step, training areas were delineated within the WorldView-2 image based on similar spectral characteristics. These polygons were then associated with specific depth classes (0.5 m intervals) using available hydrographic survey depths. The random forest classification was applied using these training areas with the WorldView-2 multispectral bands (red, green, blue and yellow) using the R statistical language [14]. Similar to the empirical approach, the classification technique is sensitive to the water depths used for model development, as well as the characteristics of the imagery. However, the technique does not require each band to have a specific relationship with depth.

#### 2.4.3. Automatic Photogrammetry

Automatic photogrammetric extraction combines the techniques of image matching with photogrammetric triangulation principles. The correlation coefficient method for stereo matching was developed in the 1970s [15]. Since then, many improvements have been made to update image matching performance [16–18]. Hirschmuller [19] introduced the concept of a Semi-Global Matching (SGM) algorithm that is based on pixel by pixel matching between images forming a stereo pair. This algorithm is particularly useful when using high-resolution stereo pairs as identical features can be more easily identified between pairs within high-resolution data. This study used an SGM approach available in PCI Geomatica 2017 to automatically derive a digital surface model for the WorldView-2 image. A light refraction correction, as described in [10], was implemented to account for the influence of refraction and the air-water interface.

Unlike the empirical and classification approaches, the automatic photogrammetric approach does not rely on in situ water depths for SDB estimation, preventing potential biases resulting from training data. However, the method is highly sensitive to site and imagery characteristics as contrasting features are required to support image matching. Through CHS's initial evaluation of automatic photogrammetry for bathymetry estimation [10], it was noted that this technique experienced challenges when attempting to match stereo pairs over homogeneous seabed areas. As well, areas containing strong temporal variability (e.g., areas containing waves) cannot be easily matched through this approach.

#### 2.4.4. Manual Photogrammetry

Photogrammetry extracts three-dimensional information based upon two different scenes acquired in close succession with varying viewing geometry. These images are used to create a three-dimensional view, similar to the effect accomplished by human eyes. Photogrammetrists make use of binocular vision and depth perception to extract the XYZ coordinates of features in the stereo pair. Positions are computed through triangulation by considering the sensor's viewing geometry at the time of acquisition. Here, a conventional photogrammetric analysis was used to manually extract depth from the stereo images using SOCET SET software. The photogrammetrist used a small number of the in situ survey data at different locations within the image to correct for the influence of light refraction and tidal effects on the estimated water depths when extracting the XYZ coordinates. Compared to the other approaches, the manual photogrammetry technique eliminates biases from sources of training information and limitations for automatic pair matching. However, this technique is strongly influenced by the individual completing stereo matching, leading to difficulties with result replication. The capability of manual photogrammetry for estimating water depth has previously been shown in [10,20].

#### 2.4.5. Empirical and Classification Approach Training

For the empirical and classification approaches, a random sample of the available survey data was selected in order to develop the empirical model and train the random forest classification. For the empirical approach, survey data was restricted to areas of the bottom which were clearly visible within the imagery. Survey data over dark areas (e.g., underwater vegetation and deep water) were excluded. For the classification approach, survey depths from 0–20 m were retained for model training, as visibility through the water column was limited at greater depths. For both approaches, 10% of the survey data, which met these restrictions, was selected for training.

#### *2.5. Level of Confidence Approach*

The proposed level of confidence method operates through understanding the level of agreement between each of the applied SDB techniques. It aims to develop a final SDB estimate where the highest number of techniques agrees within a set difference range for the largest geographical area. This is achieved through the following steps:


combined techniques. For example, averaged SDB for locations where three techniques agree would not be calculated if for the same locations, four techniques agree, as the averaged result from the four techniques would be maintained.

#### 2.5.1. Multi-Approach Combination

For step 5, noted above, there are multiple combinations of three and two techniques which can be selected. To determine which combination of techniques to incorporate into the multi-approach SDB model, as well as when to use them, the authors iteratively assessed the overall accuracy of each combination against available survey data (refer to Section 2.5.2 below). Each combination was assigned a rank based on these accuracies. This rank determined the order in which each combination contributed to the final SDB estimate. Lower ranked combinations were only used for pixels which did not meet the defined 1 m agreement level within higher ranked combinations. Combinations incorporating a greater number of techniques were automatically assigned a greater rank, even if they achieved a higher root mean square error (RMSE). Pixels which did not demonstrate agreement within 1 m for any combination of techniques were excluded from the final SDB result, as were pixels which did not contain results from all techniques. The ranking technique provides the first iteration for automatically grouping the SDB of the different techniques into one model but as described in Section 3.1, other factors can impact the results of the different approaches. In order to leverage the strengths of each approach as much as possible, it is important to complete a visual QC of each combination. The areas that match with three or more techniques are given a high level of confidence and therefore do not require manual intervention. Areas which only correlate with two or fewer approaches are given a lower level of confidence and should be reviewed visually to ensure that the best results are used in the final products.

#### 2.5.2. Accuracy Assessment

Accuracy assessment was used to determine the rank of the various technique combinations (as described in Section 2.5.1). In situ CHS survey data described in Section 2.3 was used to calculate RMSE statistics for differences between SDB estimated depths and the in situ survey depths. Overall RMSEs were calculated for depths up to 10 m in order to understand the general accuracy of each combination. Greater depths were omitted from this assessment as some approaches experienced limitations when estimating deeper water.

Accuracy assessments were also completed for the final SDB estimate resulting from the level of confidence approach and for the individual SDB techniques. This allowed for an understanding of the degree of improvement in accuracy (if any) when implementing the level of confidence approach. For these assessments, the linear error at 90% (LE90) level of confidence was calculated for each approach, allowing for an understanding of the level at which 90% of SDB errors relative to survey data would be expected to occur overall and for specific depth ranges (e.g., 0–2 m). Calculation of LE90 allowed for a direct comparison with CATZOC requirements (see Section 4.3). Overall bias was also calculated to allow for an assessment of the error trend for each approach.

#### **3. Results**

#### *3.1. Individual Approaches*

In terms of the total surface of bathymetry extracted after application of the four SDB methods, the manual photogrammetric extraction was the approach that provided the most SDB coverage, followed by classification, empirical and the automatic photogrammetric extraction approaches (Figure 4). The manual photogrammetric technique extracted 100% of the area for depths from 0–20 m, which was followed by 81% from classification, 59% from the empirical method and 39% from the automatic photogrammetric approach.

**Figure 4.** SDB results from the (**a**) manual photogrammetry, (**b**) classification, (**c**) empirical and (**d**) automatic photogrammetry techniques. Imagery © 2015, DigitalGlobe, Inc.

Table 3 presents accuracy assessment results for each of the individual SDB techniques relative to the in situ bathymetric validation data. The individual SDB techniques generate good results in general, with most achieving LE90 near 1 m for individual depth ranges. The empirical approach provided the best overall result, with a LE90 of 0.95 m for depths from 0–10 m. Unfortunately, the limited stability of the empirical approach for depths of 0–4 m (possibly due to a lower number of survey points for these depths and/or an overreliance on LiDAR measurements) prevents application of the CATZOC A2/B classification. The automatic photogrammetry approach provided the second best overall accuracy followed by the manual photogrammetry approach. The automatic photogrammetry approach achieved good results for depths from 0–4 m, while the manual photogrammetry technique generated consistent results for all depth ranges, although at a worse level of confidence relative to most of the other techniques. The random forest results provided the worst level of confidence from 0–10 m but generated good results for shallow depths (up to 4 m). Based on these results, it is important to note that each of the approaches has strengths and weaknesses related to the depth ranges and the physical properties of the study site. While each technique is capable of generating good results for some depth ranges, there are clear inconsistencies in the stability of the results between techniques.


**Table 3.** Accuracy assessment results for individual SDB techniques.

Other physical factors can play a major role in SDB accuracy and should be considered when selecting an approach. One of these factors is seabed type. For each of the investigated techniques, characteristics of the seabed can have beneficial and detrimental impacts. Dark features, commonly caused by underwater vegetation, are of particular concern for the empirical approach as they confuse dark features with deep water (Figure 5). In contrast, the automatic and manual photogrammetric techniques are not affected by dark underwater features. The photogrammetric approaches actually benefit from such dark areas, as the edges generated by these features assist with image matching [10,20]. The benthic environment also creates challenges for the random forest classification but shows improved results over the empirical approach. This is mainly because the decision tree classification approach defines a particular rule set based on training area definition. The random forest decision tree model works better for capturing non-linearity in the data by dividing overall data into smaller sub-spaces (decision trees) based on the defined training classes. This impacts the overall result and is more apparent for distinction between dark features and deep water.

**Figure 5.** Results over dark features (e.g., underwater vegetation). (**a**) WorldView-2 image with underwater vegetation visible. (**b**) Empirical approach, (**c**) manual photogrammetric technique, (**d**) automatic photogrammetric technique and (**e**) classification approach. Note the greater depth variability in the empirical and classification results, likely due to dark feature presence. Imagery © 2015, DigitalGlobe, Inc.

While the photogrammetric approaches benefit from a variable bottom, homogeneous bottom types (e.g., sand) can have a negative impact, particularly for the automatic method [10,20]. As the automatic photogrammetric technique is based on image matching, heterogeneity within the image facilitates pixel matching. Within homogeneous areas, the algorithm encounters difficulties with matching pixels, preventing a correlation from being achieved.

#### *3.2. Level of Confidence and Multi-Approach SDB Model*

Table 4 presents the overall RMSEs for each combination of techniques assessed for this study, along with the rank assigned to each combination. The percentage of the area covered by four, three and two combinations is also listed. This determined the order in which each combination was applied to generate the combined multi-approach SDB model.


**Table 4.** Overall root mean square error (RMSE) and rank of each technique combination for the level of confidence approach. The percentage of the overlap area captured by four, three and two agreeing

<sup>1</sup> AP = Automatic Photogrammetry; EM = Empirical; MP = Manual Photogrammetry; RF = Random Forest Classification. <sup>2</sup> Coverage where multiple techniques overlap with an agreement >1 m is ~0.3%.

EM, RF 0.83 10 MP, RF 0.90 11

Figure 6 presents a geographical representation of the locations where four, three and two techniques agreed. 81% of the values were given a high level of confidence (i.e., at least three approaches agree) and therefore don't need or require less QC. The values that agreed with at least two approaches represented 19% of the total overlap area; for these areas, a lower level of confidence is assigned and should be reviewed in the future by a remote sensing expert. The strengths and weakness of each approach as described in Section 3.1 should be used to make appropriate decisions regarding which combination of techniques to use. Figure 7 presents an example where the low level of confidence mask can be used as a QC tool. In this example, as the benthic environment is creating issues for the empirical and classification techniques, priority should be given to the photogrammetric approaches.

Table 5 presents accuracy assessment results for the final multi-approach SDB model. In order to compare the results against IHO requirements, LE90 was calculated for each depth range. Compared with results for the individual techniques (Table 3), the averaging of the approaches reduced the bias to negligible levels, at less than 0.2 m for each quantity of combined combinations. By combining the agreement where four and three techniques agree we retain accuracy for 0–10 m depths. Similar results are also achieved for water deeper than 10 m; where four and three techniques agree, accuracies of 1.00 m and 1.24 m are achieved respectively. Figure 8 presents a visual overview of the multi-approach SDB results.

**Figure 6.** Visualization of locations where four, three and two SDB techniques agree within 1 m. Imagery © 2015, DigitalGlobe, Inc.

**Figure 7.** Illustration of potential use of the level of confidence technique for quality control (QC). (**a**) Section of WorldView-2 imagery and (**b**) the associated number of matching SDB approaches. Note that the number of matching techniques is lower for dark bottom areas, likely due to vegetation presence. Specific techniques can be targeted for areas where two or fewer approaches agree if they can be expected to perform better for those locations (e.g., photogrammetric techniques for heterogeneous bottom areas). Imagery © 2015, DigitalGlobe, Inc.



**Figure 8.** Overview of the multi-approach SDB results. Imagery © 2015, DigitalGlobe, Inc.

#### **4. Discussion**

#### *4.1. Level of Confidence Accuracy Improvements*

From the accuracy assessment results presented in Table 5, using the level of confidence technique to build the multi-approach SDB model resulted in notable improvement to LE90 values for most depth ranges compared to results for the individual SDB techniques (Table 3). While some of the individual techniques achieved better results for specific depth ranges, the multi-approach technique resulted in less variability across all depth ranges. This suggests that by combining separate SDB results through an agreement approach, the strengths of each technique are leveraged to generate an overall improved result. By averaging depth estimates from multiple SDB techniques, outliers from individual techniques are reduced, improving the overall quality of the SDB estimates. The proposed multi-approach technique not only improved the results and confidence in the data but it also reduced required QC effort, as the area that required more detailed QC was reduced to 19% of the site. In order to simplify the QC and make the process more automatic, a seabed type classification could be implemented in the ranking system. Rule sets could be established based on water depths, accuracy and seabed classification.

#### *4.2. Geographical Coverage*

The geographical coverage of individual SDB results (Figure 4) is greater, in some cases significantly so, than the coverage offered by the level of confidence approach (Figure 8). For some locations, this is due to only a single technique producing a result. Thus, an argument can be made that some results are unfairly excluded from the level of confidence technique due to the minimum requirement for the presence of results for at least two techniques at every pixel. However, for other areas, coverage is omitted from the level of confidence approach as none of the SDB techniques produced results within 1 m of each other. In such instances, the appropriateness of applying any SDB technique to these locations can be questioned. This is an important consideration for HOs where confidence in the information is a significant concern. While use of the level of confidence approach will result in reduced geographical coverage, HOs should be able to have greater confidence in the locations which retain coverage after the approach is applied.

#### *4.3. International Standard Compatibility*

The suitability of hydrographic surveys for charting applications is determined using the IHO S-57 standard [11]. This standard defines CATZOC levels which contain depth accuracy specifications for specific depth ranges. Table 6 presents a summary of the accuracy requirements for each CATZOC level for depths up to 30 m. SDB estimates will need to be assigned a CATZOC level to allow for incorporation into nautical products.

**Table 6.** Required depth accuracies for International Hydrographic Organization (IHO) CATegory of Zones Of Confidence (CATZOC) levels [11].


Comparison of the CATZOC definitions with SDB results derived for each technique (Table 3) shows that all SDB estimates meet the CATZOC C level. For certain techniques and depth ranges, CATZOC A2/B accuracies are achieved. However, inconsistencies in the results between techniques across various depths highlight the instabilities present within individual approaches, limiting the degree of confidence which can be placed in any single approach. While each approach has the potential to achieve CATZOC A2/B level accuracy, there remains the potential for the presence of outliers throughout the area of coverage, particularly where survey data for validation is not available. This represents a significant hurdle for the incorporation of SDB within hydrographic products, as HOs require a high degree of confidence in the data they are using.

The application of the level of confidence approach resulted in improvements to overall accuracy as well as accuracies for specific depth ranges (Table 5). While the improvement relative to individual SDB techniques is generally not substantial, result consistency and thus overall confidence is increased as more than one approach is known to be producing a similar result for the same location. This is particularly important for areas outside of in situ survey data coverage as it provides a partial form of validation. While not as robust as comparisons with in situ information, knowing that at least two or more SDB techniques agreed within a defined level is significantly better than having no understanding of SDB appropriateness outside of survey data locations. This represents an important consideration for HOs, especially as they work to assign CATZOC classifications to SDB estimates. By applying this technique, the results demonstrate that 81% of the total common SDB surface could now meet the CATZOC A2/B level requirements, the rest of the SDB surface could be classified as CATZOC C.

#### *4.4. Hydrographic Office Use*

CHS envisions several potential uses of the level of confidence approach by HOs to increase the incorporation of SDB within nautical products:


#### *4.5. Future Research*

While the proposed approach suggests an interesting method to increase overall confidence in SDB results, its application would benefit from additional investigation:


#### **5. Conclusions**

For many years, SDB has offered the potential of providing an additional accurate, inexpensive source of bathymetric information. While representing a potentially important data source for many HOs, uptake of SDB information for use in official nautical products has been limited. Much of the cause for this limited uptake is due to the level of confidence which HOs are comfortable placing in data sources. Identified as a common theme during the first IHO supported International Hydrographic Remote Sensing workshop (the workshop was hosted by CHS in collaboration with Service Hydrographique et Océanographique de la Marine (SHOM) and the National Oceanic and Atmospheric Administration (NOAA), 18 to 20 September 2018, Ottawa, Ontario, Canada.), HOs will not utilize SDB information unless it can be validated.

Within this paper, CHS has proposed a level of confidence approach to allow for combinations of SDB estimates obtained from multiple techniques where they agree within a defined level (e.g., 1 m). Through a ranking scheme, various combinations of techniques are built up to generate a final SDB estimate where agreement is highest amongst the most techniques for the largest possible area. Results presented in this work illustrate that the level of confidence approach improved LE90 statistics overall and for individual depth ranges in most cases. CHS's overall confidence in the level of confidence results was also increased, as for all locations containing SDB estimates, at least two applied techniques demonstrated agreement within 1 m.

CHS believes the proposed technique will allow HOs to obtain greater confidence in SDB results, allowing for wider implementation within nautical products. HOs will be able to define a required agreement level and also determine an appropriate number of techniques which should agree to allow SDB estimates to be produced for an area. The approach also allows for a better understanding of

the appropriateness of SDB techniques within areas where no in situ survey data is present, allowing for a form of validation for the entire geographic coverage of SDB estimates. While representing an interesting first step, future research will be required to understand the repeatability of the approach, the potential for adding weighting approaches and how the technique could be used to combine SDB results obtained from multiple images.

**Author Contributions:** Conceptualization, Supervision, Project Administration and Funding Acquisition, R.C.; Methodology, Formal Analysis and Investigation, R.C., R.A., M.S., M.-A.F. and Y.S.; Validation and Visualization, M.S.; Writing—Original Draft Preparation and Writing—Review and Editing, R.C., R.A. and M.S.

**Funding:** CHS received support through the Canadian Space Agency's Government Related Initiatives Program to complete this work.

**Acknowledgments:** The authors would like to thank Matúš Hodúl, a M.Sc. graduate from the University of Ottawa, who worked at CHS during his studies and helped in the development of the automatic photogrammetric extraction process used in this study. The authors would also like to acknowledge the efforts of Effigis Geo Solutions, who completed manual photogrammetric image processing through a contract with CHS. Additionally, the authors wish to sincerely thank the Canadian Space Agency for their continued support, without which this work would not have been possible.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Method of Watershed Delineation for Flat Terrain Using Sentinel-2A Imagery and DEM: A Case Study of the Taihu Basin**

**Leilei Li 1,2,\*, Jintao Yang 2,3 and Jin Wu 2,3**


Received: 4 September 2019; Accepted: 24 November 2019; Published: 26 November 2019

**Abstract:** Accurate watershed delineation is a precondition for runoff and water quality simulation. Traditional digital elevation model (DEM) may not generate realistic drainage networks due to large depressions and subtle elevation differences in local-scale plains. In this study, we propose a new method for solving the problem of watershed delineation, using the Taihu Basin as a case study. Rivers, lakes, and reservoirs were obtained from Sentinel-2A images with the Canny algorithm on Google Earth Engine (GEE), rather than from DEM, to compose the drainage network. Catchments were delineated by modifying the flow direction of rivers, lakes, reservoirs, and overland flow, instead of using DEM values. A watershed was divided into the following three types: Lake, reservoir, and overland catchment. A total of 2291 river segments, seven lakes, eight reservoirs, and 2306 subwatersheds were retained in this study. Compared with results from HydroSHEDS and Arc Hydro, the proposed method retains crisscross structures in the topology and prevented erroneous streamlines in large lakes. High-resolution Sentinel-2A images available on the GEE have relatively greater merits than DEMs for precisely representing drainage networks and catchments, especially in the plains area. Because of the higher accuracy, this method can be used as a new solution for watershed division in the plains area.

**Keywords:** watershed division; Sentinel-2A; Google Earth Engine (GEE); Taihu Basin; hydrology; plains area

#### **1. Introduction**

Drainage networks are essential to geospatial assessment, basin analysis, and applications for catchment delineation, flow statistics, flood risk assessment, and climate modeling [1,2]. They also play an important role in estimating the transmission of pollution and nutrients [3]. Drainage networks are a fundamental condition for watershed delineation and an indispensable component in hydrological modeling. Therefore, reasonable and accurate watershed delineation is a precondition for runoff, sediment transport, water quality simulation, and basin resource management [4].

With the arrival of DEM, there has been a growth spurt in studying DEM-based extraction algorithms for drainage networks [5–7]. DEM, which represents the continuous variation of relative elevation values in each pixel, can be used to identify a river direction or flow path on the ground surface. Many algorithms, including the drainage networks algorithm [8] and the delineating watersheds algorithm [9], are developed to automatically extract essential hydrological features with DEM. Those algorithms require filling DEM to avoid depressions on the slope and generate a reasonable stream network in the preliminary stage. The stream network process usually starts by allocating flow directions for each DEM cell, then analyzes the flow accumulation, and finally selects those cells that have a total higher threshold of flow accumulation than a defined value [10]. Many factors, including the DEM spatial resolution, the calculated algorithm, and the physical characteristics of the basin, can directly affect the accuracy of drainage networks derived from DEM data [11,12]. Using DEM with an algorithm has obvious advantages, because the watershed processing takes less time and is not influenced by human subjectivity.

Although DEM is widely used, many problems are inherent in its practical application. DEMs usually include shallow depressions, i.e., areas surrounded by higher elevation that could occur naturally or artificially in a river plain landscape. For the determination of stream networks, an effective method of eliminating pits or depressions is the stream burning algorithm [13]. The stream burning algorithm often identifies river channels or lakes that are invisible in the DEM, preventing serious errors in the streaming. However, many features, such as hardened roads, low dams, artificial rivers, and large lakes could result in significant depressions and limit the accuracy of the results, especially in low-relief areas [14]. In addition, the most freely available open-source DEM datasets have a medium resolution of 30 m, which limits the spatial accuracy. Therefore, traditional DEM processing can generate unrealistic drainage network models due to large depressions and subtle elevation differences in local-scale areas, especially on river plains [15]. Automated watershed delineation with medium resolution DEM is not as accurate as high-resolution raster images of drainage networks for a plains area.

Researchers have made great efforts to improve the accuracy. The manual division of watershed boundaries from topographical maps is a very intense, time-consuming process, especially for large-scale studies, but is considered to be a very effective and accurate method [16]. Similar extraction and comparison of drainage networks could be obtained automatically from LiDAR-derived DEM or manually by drawing rivers from aerial photographs [17,18]. Automated delineation from LiDAR-derived DEM in plains areas have many problems. Stream networks extracted from satellite images or aerial photography are held to be truthful channels [19], primarily because the more detailed hydrological features in aerial photographs guarantee better drainage networks as compared with results obtained from DEM. Although high-resolution aerial photographs are ideal for extracting precise drainage networks and catchments, they are very expensive and not available on a large scale.

Due to their high spatial resolution and multispectral bands, optical remote sensing images, such as TM, SPOT, and Sentinel-2 images, can be used for water extraction. There are two common methods of extracting waterbodies. The first method which includes single threshold segmentation and an Otsu threshold segmentation algorithm, is to process the waterbody index threshold on the image such as normalized difference water index (NDWI), the modified NDWI (MNDWI), or the enhanced water index (EWI) index. Another method uses supervised classification, which employs training samples to select different trains, identify the characteristic parameters, and establish a discriminant function to distinguish different classifications. The second method usually involves extensive calculations, manual intervention, and a long time to obtain training samples.

For this situation, high-resolution satellite images, freely available on GEE, are a better choice for extracting drainage networks and catchments, especially for plains areas [20,21]. We develop a method to correct drainage networks by extracting rivers, lakes, and reservoirs. In addition to DEM, high-resolution Sentinel-2A images are used as input data. First, we extract drainage networks from Sentinel-2A images of the Taihu Basin on the GEE Platform. Second, we obtain results that show improved quality of watershed delineation as compared with catchments derived from the DEM. Compared to medium-resolution DEM of drainage networks and catchments, the objective of this study is to demonstrate the relative advantages of the proposed method for making precise and accurate watershed delineation in plains areas using Sentinel-2A images available on the GEE, taking the Taihu Basin as a case study.

#### **2. Materials and Methods**

#### *2.1. Study Area*

As shown in Figure 1, the Taihu Basin, which occupies an area of about 36,500 km2, is situated in the downstream area of the Yangtze River basin, surrounded by water on three sides comprised of the Yangtze River in the north, the Hangzhou Bay in the south, and the East China Sea in the east. It extends across four provinces of China which are Jiangsu, Zhejiang, Anhui, and Shanghai. The Taihu Basin is divided into four geomorphologic types, high mountains, gentle hills, alluvial fans, and aggraded plains which are mainly dominated by gently sloping plains covering more than 70 percent of the total area.

**Figure 1.** The geographic location of the Taihu Basin.

Its overall terrain is distinctly variable, with elevations varying from 0 to 1600 m. The Taihu Basin has dense artificial rivers and numerous lakes, the largest of which has an area of 2338 km2. The river network in the plains area is dense, crisscrossed, and surprisingly similar to an urban road network in its overall pattern. The average annual rainfall is 1177 mm for the total basin, the spatial distribution of which gradually decreases from 1408 mm in the southwest Tianmu mountains to less than 1100 mm in the eastern coastal area and northern plains [22]. Since it has one of the highest gross domestic products (GDP) and the most urbanized basin area in China, many projects have been allocated for watershed-based research since the 2000s and earmarked until 2025. Due to serious historical water pollution and flooding in the plains area, precise watershed delineation is of great significance for water conservation in the region [23,24].

#### *2.2. Data Source*

#### 2.2.1. Sentinel-2 Images

GEE provides an online coding environment that can guarantee relatively rapid, distributed-cloud geospatial processing and analysis of large satellite image datasets that are convenient for everyone. It provides access to publicly available Sentinel-2A images and large-scale remote sensing analysis algorithms [25]. Therefore, GEE is used to code, analyze, and generalize the principal rivers, lakes, and reservoirs that could be derived from Sentinel-2A imagery of the Taihu Basin.

Corresponding to the latest temporal resolution and highest spatial resolution data available freely in GEE, Sentinel-2A OLI images that completely covered the whole Taihu Basin were collected from 1 May 2018 to 31 December 2018. Sentinel-2A images, obtained from the GEE platform throughout 2018, were considered to be Level 2 products and archived as top-of-atmosphere reflectance collections. The preprocessing work, including stitching images and masking clouds, was coded and executed through the developed application programming interface (API) in the GEE platform [26]. To ensure complete coverage with little cloud cover, less than 10% of scenarios were selected from the image sequences, using each specific path and row shown in Figure 2. To reduce the influence of seasonal variability on the identification of waterbodies, satellite images in summer were preferred, because plentiful water contributed to the extraction of rivers and lakes during this period. A total of three rows and four paths in the Sentinel-2A images, and approximately 15 scenarios in 2018, were obtained and used to satisfy the basic coverage needs of the Taihu Basin.

**Figure 2.** The (**a**) spatial distribution and (**b**) cloud coverage in the Sentinel-2A images for 2018.

#### 2.2.2. DEM Data

With approximately 30 m resolution and high vertical accuracy (less than 8 m), the SRTM 1 arc-second DEM data, acquired from the United States Geological Survey, is produced for the year 2000 and provided the GCS\_WGS\_84 geographic coordinates for the Taihu Basin. The SRTM data is easily downloadable from https://earthexplorer.usgs.gov/. The entire DEM is composed of eight scenario small DEM images of the Taihu Basin, whose respective names are n30\_e119\_1arc\_v3, n30\_e120\_1arc\_v3, n31-e119\_1arc\_v3, n31\_e120\_1arc\_v3, n31\_e121\_1arc\_v3, n32\_e119\_1arc\_v3, and n32\_e120\_1arc\_v3.

#### *2.3. Processing Methods*

The overall approach used in this study for the watershed delineation of the Taihu Basin is shown in Figure 3. Significant characteristics of the proposed method include the following: First, GEE, a distributed-cloud satellite image processing platform, was applied to process a number of Sentinel-2A images covering the Taihu Basin within a short timeframe. It is a relatively automatic, precise means of extracting waterbodies and determining the drainage network mainly from Sentinel-2A images, rather than from the DEM data. Second, it delineated catchments for each hydrological element by considering the modifying flow direction of each river, lake, and reservoir, in the plains area. These steps are discussed in detail in the following sections.

**Figure 3.** The overall method of watershed delineation in this study.

First, we obtained the actual rivers, lakes, and reservoirs from the Sentinel-2A images to make up a drainage network, as opposed to results that were extracted from the DEM. Secondly, for the key point of watershed delineation we acquired a drainage network and modified the flow direction of the hydrological features to obtain accurate catchments in the Taihu Basin. Therefore, the major process, shown in Figure 3, was as follows: (a) Calculate the MNDIW index with GEE; (b) use the Canny edge detector to extract water bodies; (c) determine the flow direction for each river, lake, and reservoir; (d) determine the flow direction of the overland flow; and (e) conduct subwatershed delineation for each river, lake, and reservoir.

#### 2.3.1. MNDWI Index Introduction

Most methods that detect waterbodies from multispectral satellite images rely on the fact that the absorbing radiation effect of waterbodies is evident in the near-infrared and even longer bands. This principle enables them to highlight waterbodies clearly with raster band calculation, using the NDWI index to represent differences in the spectral curve and highlight the waterbody information in the image. Due to its use of the near-infrared band to replace the shortwave infrared band, the MNDWI index appears to be more sensitive than the NDWI index, resulting in a greater ability

to highlight waterbody information [27,28]. We processed Sentinel-2A images to determine the surface reflectance and followed formula (1) to compute the MNDWI for surface water detection on GEE. In fact, the spectral reflection spectrum curve of water usually varies, owing to different subtle objects in the water. This makes it hard to adopt a spectral index with a fixed threshold to detect waterbodies on the ground. Numerous regional errors could be generated if a zero threshold is applied to detect waterbodies from water index images. There are regional gaps in the spectral properties of different waterbodies.

$$\text{MNDDWI} = \frac{\rho\_{\text{Grenen}} - \rho\_{\text{SWIR}}}{\rho\_{\text{Grenen}} + \rho\_{\text{SWIR}}} \tag{1}$$

Green indicates a green band (band 3, 10 m) and SWIR indicates a short infrared band (band 11, 20 m) in Table 1. Due to automatically and dynamically projecting the 20 m (band 11) to the 10 m (band 3) resolution in GEE, we were able to directly calculate the MNDWI between the green and SWIR bands without considering the resolution difference of the Sentinel-2A images. According to previous studies, varying thresholds in the MNDWI index have generally required manual threshold adjustment and limited their applicability. The challenge was to establish a varying optimal threshold method, allowing for more fully automated waterbody detection. On the basis of a statistical histogram of all the MNDWI values in a satellite image, histogram-based threshold detection methods help to overcome this problem [29].

**Table 1.** An introduction to the bands used in the Sentinel-2A images.


2.3.2. Canny Edge Detector for Extracting Water Polygons

The Canny edge detector is a widely used method applied in visualization computer systems, which allows for accurate edge detection in images [30]. The scope of this study employed histogram-based threshold methods, using a Canny edge filter to detect the boundaries between water and nonwater pixels. The Canny edge algorithm, which is applied to a MNDWI image with a threshold, revealed water edges that are in rapidly changing areas of an index value [31,32]. In terms of the MNDWI index, this situation generally involves obvious value changes in image pixels when shifting from one neighborhood to another. Waterbody and overland pixels are distinguished by applying buffer analysis to count the MNDWI values, generate a histogram, and determine a segmentation threshold. In a real situation, the MNDWI histogram presents bimodal statistical distribution, and therefore an obvious distinction could be made between waterbody and overland pixels.

In this study, we used binary images where the boundaries between water and nonwater pixels were clearly defined. The algorithm operated as follows: (1) For the image smoothing, generally, edge detectors are prone to noise, therefore, the image was smoothed with a square-sized Gaussian structural element, usually of 5 × 5 size; (2) for the gradient intensity calculation, the gradient magnitudes and directions were computed, and the gradient magnitude defined the edge, with a high gradient magnitude generating a rapid change in color (an edge), in contrast to a low gradient magnitude. The gradient direction, as the name implies, defined the orientation of the edge; and (3) for the hysteresis thresholding, this stage removed small pixel noises based on two threshold values of the intensity gradient, and discarded pixels below a certain threshold.

GEE supplies a built-in interpretation of the Canny edge detection algorithm, which is utilized to compute the edges [33]. To estimate the boundaries of a river under cloudy circumstances, multiple historical Sentinel-2A satellite images were used to extract composite waterbodies, representing the boundaries of rivers and lakes. Multitemporal composite results were used to identify segments of a river that were initially unknown due to the presence of clouds. As shown in Figure 4, the process of extracting waterbodies comprised the following five steps: (a) calculate the MNDWI index from the Sentinel-2A images, (b) detect sharp edges with the Canny edge detector, (c) count the threshold values, (d) extract the waterbodies using the statistical threshold in the MNDWI image, and (e) compose waterbodies from multitemporal Sentinel-2A images. The detailed code can be found in Appendix A.

**Figure 4.** The process of extracting waterbodies based on the modified normalized difference water index (MNDWI) index and Canny algorithm: (**a**) False color composition of the Sentinel-2A image, (**b**) MNDWI index image, (**c**) the histogram threshold based on the MNDWI values, and (**d**) extracted river polygons and their partial enlargement results

After extracting the waterbodies from the GEE platform, the rivers, lakes, and reservoirs were generalized in the ArcMap software. First, each river segment was simplified into a single centerline and the width calculated using the "Collapse Dual Lines to Centerline" tool in ArcMap. Second, the lakes were manually distinguished from reservoirs based on the fact that lakes in the plains area were usually less than 10 m deep or even below zero, while reservoirs in the mountain area were higher than 20 m. Finally, waterbody polygons were preliminarily processed so that they could be analyzed and organized as linear rivers, polygonal lakes, and polygonal reservoirs based on their hydrological features in ArcMap.

2.3.3. Flow Direction Determination for Each River, Lake, and Reservoir

The well-known D8 algorithm was used to determine the flow direction from each cell to its downslope neighbor or neighbors on a topographic surface [34]. This method can track and represent the flow from each pixel to neighboring pixels by using eight discrete flow angle values. Due to the large numbers of ring-shaped and well-shaped river structures in the plains area, a traditional D8 algorithm, used to express dendritic nodes, could not properly digitize the flow direction of river nodes connecting more than three rivers [35]. In order to digitize multiple directions uniformly into a raster value, river nodes connected to more than three river segments had to be recorded. Therefore, the multi-D8 algorithm, which is an extension of the traditional D8 algorithm, was used to express multiple flows from a crossed node. The multi-D8 algorithm could express multiple flow directions from a node cell in a binary exponential manner, so there was no confusion when multistreams were encoded and decoded in a digital representation [36]. As shown in Figure 5, there was a detailed explanation for 30 (2<sup>1</sup> + 22 + 2<sup>3</sup> + 24) equaling four flow directions in a node cell.

**Figure 5.** The multi-D8 algorithm extended from the D8 algorithm in a multiple flow node cell.

Because the spatial resolution of Sentinel-2A images was 10 m, the principal rivers represented by image pixels in the GEE were usually wider than 100 m in the results. In a sense, this was a generalization of rivers from satellite imagery, neglecting narrow rivers of the same scale. In this study, a variety of means were used to comprehensively assess and analyze the flow direction of each river, such as hydrological observation data, DEM height differences, and node directions for each river segment. First, hydrological observation data from observation stations directly determined the flow direction of each generalized river, but it could not cover all the rivers in the Taihu Basin. Second, node junctions consisting of rivers, lakes, or reservoirs in a hydrological network were useful for determining the flow direction of undefined river segments. If a node junction was composed of three main rivers, an undefined river flow could be estimated from the flow directions of two other rivers [37]. Third, the flow direction of each main river was estimated according to the average elevation difference between the start point and end point values, which could be extracted from the DEM data of a river segment [38]. Following the above three steps, we determined the flow direction of each main river in the Taihu Basin. The hydrological observation data was provided by the Taihu Lake Basin Hydrological Information Service System on the internet, which was maintained by the local water conservancy bureau and can be visited at http://218.1.102.107:8100/indexWater.html. The Taihu Basin Hydrological Information Service System covered the 98 hydrological observation sites that collected daily hydrological data for the study area, including flow rate and water flow direction.

Since the catchment of a lake or reservoir usually had a specific shape, it was not necessary to determine flows for a cell overlapped by itself. Instead, a unique identifier was assigned to each cell that was overlapped by an integer value of 500 for a lake and 600 for a reservoir. As shown in Figure 6, the flow direction of a lake or reservoir was represented by its edge line with an arrow.

**Figure 6.** Flow direction for each (**a**) river, (**b**) lake, or reservoir, and node using the multi-D8 algorithm.

2.3.4. Flow Direction Determination for the Overland Flow

All the flow direction grids comprised river segments, lakes, reservoirs, and overland flows in the Taihu Basin. The flow direction of the overland flow derived from the overland DEM is shown in the Figure 7. First, the overland DEM need to be pretreated by masking waterbody grid cells that include the main rivers, large reservoirs, and lakes, then filling the pits or depressions of the overland DEM in ArcMap. Second, the flow directions of the overland flow can be digitized through the D8 algorithm and calculated from the overland DEM in ArcMap. Therefore, we cut the original DEM into pieces of overland polygons, which excluded the geometry of rivers and lakes; calculated the D8 flow direction for each overland flow; and merged the flow direction grids using the hydrological tools in ArcMap. Compared to calculating the flow direction directly from the whole DEM, this method prevented a large number of straight flow directions in depression areas and made streams more precise in the overland flows.

**Figure 7.** Illustration of the flow direction digitization process for the overland flow.

#### 2.3.5. Subwatershed Division

The Taihu Basin can be divided into three kinds of catchments: Reservoirs, lakes, and onshore catchments. The sub-basin partitioning method was distinct from traditional methods that made use of flow direction, flow accumulation from DEM, and outlets to the complete catchment division in the Taihu Basin. This method extracted waterbodies based on the high-resolution Sentinel-2A satellite images and composed a drainage network for the plains area. The catchment area was delineated for each actual main river in this section.

On the basis of the above modified flow directions of each river, node, lake, reservoir, and overland flow, an entire flow direction raster grid covering the Taihu Basin was generated by the overlapping process of ArcMap, as shown in Figure 8. First, the sub-basin division searched the inflow paths of each grid in the upward direction, from the end to the start point of a river segment. Second, an inflow path for each grid in a river segment linked all the same inflow adjacent grids of the overland flow into a chain. Third, all inflow paths flowing into a unique river were marked with the same integer identification (ID). The sub-basin of each river segment consisted of grids with the same ID. Finally, a few sub-basins with an area of less than 10 km2 threshold were combined into neighboring sub-basins to meet the requirements. In this way, we could accurately divide the two pavement catchments of each river, the catchment of each lake, and the catchment of each reservoir.

**Figure 8.** Illustration of the subwatershed extraction process for river segments.

It was clear that each river segment had only two catchments, on the left and right. Each river node was connected to at least three catchments, which had the same directions as their river segments. The catchments' topology, which was important for the calculation and analysis of the distributed hydrological modeling, was based on the flow directions in the river network and organized according to hydrographic nodes. Following the determined flow directions of the stream network, the relationships between upstream and downstream nodes became clear. Therefore, it was easily possible to derive the up–down relationship of the catchments from the downstream and upstream nodes in the river network. In this way, we established a watershed delineation with clear upstream and downstream relationships and provided accurate catchment data for watershed application in the Taihu Basin.

#### **3. Results**

#### *3.1. Accuracy Analysis for Waterbody Extraction*

We compare the effects of waterbody extraction among rivers, lakes, and reservoirs in different terrains. Although some very narrow rivers have broken, the overall structure of rivers can be shown in the river plain. The effects make it better in reservoirs and lakes.The effect is shown in Figure 9.

**Figure 9.** A comparative analysis of waterbody extraction between rivers and lakes, mainly including true color images, MNDWIs, and extracted polygons. They are listed as: (**a**) artificial river extraction for the plains river network, (**b**) natural river extraction for the plains river network, (**c**) reservoir extraction at high altitude, and (**d**) small lake extraction in the plains.

The Global Land Cover 30 (GLC30) is produced using multi-seasonal Landsat data based on the random forest-based mapping framework of Tsinghua University, for which the overall accuracy of the global land cover mapping is 73% at 30 m spatial resolution. The classification types include cropland, forest, grassland, shrub land, wetland, water, impervious surface, bare land, and snow. To validate the effect of waterbody extraction, this study selected validation points from GLC30, which can be downloaded from http://data.ess.tsinghua.edu.cn/ [39]. A random sampling method was used to select 500 water surface points and 500 nonwater surface points from GLC30, in which the water and nonwater verification points were evenly distributed across the Taihu Basin. The sample points that were extracted from the GLC30 data were checked and modified manually using the Google Earth software.

$$OA = \frac{1}{N} \sum\_{i=1}^{M} n\_{ii} \tag{2}$$

$$P\_{\mathcal{E}} = \frac{1}{N^2} \sum\_{i=1}^{M} n\_{i+} \ast n\_{i+} \tag{3}$$

$$Kappa = \frac{P\_o - P\_c}{1 - P\_c} \tag{4}$$

The overall accuracy (OA) was calculated from the confusion matrix using Equation (2), which essentially represented the proportion that was classified correctly from all the reference sites. Cohen's kappa coefficient was generated from a statistical test to evaluate the accuracy of a classification, which was calculated from Equations (3) and (4). In Equation (4), Po represented the ratio of the total sample, i.e., the ratio of the correct classification of each category to the total sample, which meant that Po was the overall accuracy. The user's accuracy was calculated by taking the total number of correct classifications for a particular class and dividing it by the row total. In this study, the user's accuracy was calculated from the above error matrix. As shown in Table 2, the user's accuracy of water was 90.3%, the overall accuracy was 89.3%, and the Kappa coefficient was 0.793.

**Table 2.** The confusion matrix and accuracy evaluation of waterbodies in the Taihu Basin.


#### *3.2. Drainage Networks and Catchment of the Taihu Basin*

This section introduces the results of watershed division analysis through the proposed method in the study area, which occupies a large plain in southeastern China. The northwestern and southern parts of the Taihu Basin are dominated by undulating hills and low mountains, accounting for approximately 20% of the total area, while the central and eastern parts are flatlands. The wide plain area is filled with dense river channels, large lakes, and a considerable number of reservoirs, which makes the hydrological elements relatively complex in the Taihu Basin. The hydrological elements were obtained and identified from Sentinel-2 images using the Canny algorithm on the GEE platform. Considering the resolution limitation of the DEM and Sentinel-2 images, some small lakes and reservoirs with an area of less than 10 km<sup>2</sup> were replaced by the centerlines of each waterbody. It was important to maintain the actual topological connections between each linear river, polygonal lake, and polygonal reservoir, which were vital for the pretreatment processing of the drainage network.

As shown in Figure 10, almost 2291 river segments with a width greater than 100 m were extracted, filtered, and generalized from the Sentinel-2 images. Seven lakes and eight reservoirs, which had clearly topological entrances and outlets to river segments, were extracted. The final catchments in the Taihu Basin were divided into the following three types using the proposed approach: Lake catchments, reservoir catchments, and overland catchments. The flow directions in a stream were expressed by arrows along the centerline of a river, or the edge line of a lake or reservoir. Because the streams used in the catchment division were the real drainage networks from the Sentinel-2A images, they retained the original crossing river structures of the Taihu Basin. Each river segment had two overland catchments on the left and right. The divided subwatersheds were detailed, comprising a total of 2306 catchments across the entire basin. These actual catchments with topology connectivity provided more accurate fundamental spatial data for distributed hydrological analysis.

**Figure 10.** Drainage networks and catchment results extracted from Sentinel-2A and DEM.

#### *3.3. Comparison of the Stream Burning Methods in Arc Hydro and HydroSHEDS*

Simultaneously, this paper compared the sub-basin division results of the stream burning methods in the Arc Hydro and HydroSHEDS tools. In order to contrast this way with the stream burning method frequently used in the plains area, we used rivers and DEM to divide the Taihu Basin through Arc Hydro tools with an integrated stream burning method. The stream burning method in Arc Hydro mainly includes the following steps: (1) generating filled and flawless DEM, (2) using the DEM reconditioning tool to reduce a certain elevation artificially for vector rivers and burn the rivers into 30 m DEM, (3) generating flow direction grids and extracting streamlines based on a flow accumulation threshold, and (4) using streamlines and outlets to extract catchments using the Watershed tool. According to several processing experiments, a good result was obtained using the burning algorithm when the flow accumulation was set to 4000 ha in Arc Hydro. This result brought the source of the drainage network closer to the actual situation. The results of the stream burning method are shown in Figure 11a.

HydroSHEDS, a hydrological product that offer local- to global-scale catchment data, is mainly generated from SRTM DEM data at three arc-second spatial resolution [40]. Several procedures, such as eliminating small pits from elevation data (the sink filling method), combining digital river maps (the stream burning method), and calculating flow directions, are used in the data generation process and have been well documented in the literature. HydroSHEDS generates high-precision drainage network and subwatershed models, because it uses the above methods to produce hydrologically conditioned elevation and reduce errors [41]. As shown in Figure 11b, the streamlines derived by the proposed method were contrasted with HydroSHEDS datasets to evaluate the accuracy of the streamlines, which were downloaded from the official website https://www.hydrosheds.org.

The comparative results in Figure 11 showed that the streamlines obtained by stream burning in Arc Hydro and HydroSHEDS were generally similar in their geometric shape and spatial distribution. However, compared to actual river streamlines drawn from the Sentinel-2A images, HydroSHEDS did not solve the problem of ring-shaped and well-shaped river networks in the plains area. Meanwhile, it can clearly be seen from Figure 11 that streamlines obtained by stream burning in Arc Hydro and HydroSHEDS could not eliminate the large numbers of incorrect streamlines in some large lakes in the Taihu Basin. These fallacious streamlines, which were inconsistent with the actual situation, could introduce serious errors into hydrological applications used in the Taihu Basin.

**Figure 11.** Comparative effects of (**a**) stream burning and (**b**) HydroSHEDS using the proposed method.

#### **4. Discussion**

#### *4.1. Di*ff*erences from the Stream Burning Method*

From the perspective of data precision, we find that there are some obvious differences between the catchments extracted using a traditional method and those obtained by using the proposed method in the Taihu Basin. The traditional method depends on 30 m DEM to define the flow accumulation and extract the stream network using a critical threshold value in Arc Hydro [42,43]. However, a digital hydrological network, which is composed of rivers, lakes, and reservoirs, is precisely extracted from 10 m high-resolution Sentinel-2A images covering the entire Taihu Basin. The spatial resolution of the satellite images and actual local hydrological elements have a significant impact [44]. The accurate and realistic rivers and lakes, which are important components of hydrological networks, are extracted from the Sentinel-2A images on the GEE and prevent many cumbersome manual vectorizations.

Compared with previous studies, the stream network extracted by the DEM in Arc Hydro is a simple dendritic structure, with flow direction represented by the standard D8 algorithm for each cell [45,46]. In fact, the usually circular and well-shaped river structures of the southern plains in the Taihu Basin are ignored. Meanwhile, multi-flow directions in the crossed nodes of rivers are redefined and digitized by the multi-D8 algorithm in the flow direction determination process of this study. The proposed method directly adjusts the flow direction of each river, lake, and reservoir in the real stream network, hence, it could avoid many unrealistic river channels in the lake area and retain the complete structure of the plains area. With the help of high-resolution Sentinel-2A satellite images, this automated method proves conducive to the timely and accurate division of the catchment area.

#### *4.2. Catchments for Hydrological Simulation*

The detailed catchment results are derived from rivers and lakes in the Sentinel-2A images and retain their crossed and annular structure. From the perspective of hydrological simulation, a simple dendritic structure is usually needed. In fact, the rivers could simply be changed to a dendritic structure by interruption and used for the hydrological slope and sub-basin processes in the SWAT model for a smaller area [47–49]. Simultaneously, the detailed extracted catchments that corresponded to these rivers could be merged and used in the hydrological runoff simulation. This finer division of catchments in the local plains may provide improved accuracy of rainfall runoff simulations.

The drainage network results could also be used for the main river channel calculation in the HEC-RAS model, which has the ability to calculate crossed and annular rivers in the plains area [50,51]. This offers the possibility of performing two joint models to improve flood simulation accuracy based on the extracted drainage network and catchments [52].

#### *4.3. Further Improvement*

Some aspects of the proposed method require further improvement in order to more readily obtain accurate catchments in the future. The finer the degree of the stream network, the greater the workload required to preprocess the flow direction cells. The flow direction of rivers that are assessed using hydrological data from observatory stations may be the most accurate method, and therefore more extensive data should be collected and utilized to correct the flow direction of rivers in the future. In addition, the consideration of artificial highways and railways, which form relatively high-terrain pseudo-topographic valleys in the plains area, could make the division of overland flow catchments more reliable and reasonable [53,54]. Water bodies derived from satellite images at different resolutions is corresponded to catchment areas of different scales. In addition, the multiresolution satellite images used in the GEE could be further considered and tested for their usefulness in extracting multiscale drainage network models and watershed delineation in the future.

#### **5. Conclusions**

Due to large depressions and subtle differences of DEM in the plains area, catchments division are difficult. The traditional method, using stream burning to delineate the watershed, is not sufficiently accurate for actual local hydrological purposes. This study proposes a method to delineate the watershed using high-resolution Sentinel-2A images as the data source for the plains area. On the basis of the Canny edge detection algorithm, rivers, lakes, and reservoirs are automatically extracted from Sentinel-2A satellite images, then combined to form the drainage network. Catchments are extracted for each hydrological element by directly modifying the flow direction of each river, lake, reservoir, and overland flow. A watershed is delineated into three types of catchment, i.e., lake, reservoir, and overland flow catchment in the Taihu Basin. The proposed method is found to be suitable for the Taihu Basin plains study area. Compared to the traditional method, the new method is more consistent with the actual precise watershed delineation. Therefore, this study provides a new method for extracting the subwatershed from Sentinel-2A images and providing more accurate basic data for subsequent hydrological applications, especially in the plains area.

**Author Contributions:** L.L. and J.Y. developed and planned all the experiments; L.L. and J.Y. performed the code programming work and method validation, provided the results discussion, and wrote the majority of sections in this manuscript; J.Y. and J.W. supplied recommendations and comments and edited the manuscript; all the authors revised and adjusted the manuscript collaboratively.

**Funding:** This research was funded jointly by the National Natural Science Foundation Project of China with funding reference 41590845. We also thank the funding sponsor from the State Key Laboratory, the GEE platform for the Sentinel-2 satellite imagery, and all those people who have made valuable comments and suggestions for this article.

**Acknowledgments:** The authors are grateful to ChenHu Zhou for his valuable suggestions on improving our work. We thank our colleagues in the State Key Laboratory of Resources and Environment Information System at the Institute of Geographic Sciences and Natural Resources Research for their assistance with this study. We thank the journal's editors and anonymous reviewers for their kind comments and valuable suggestions to improve the quality of this paper.

**Conflicts of Interest:** The authors declare that they have no conflicts of interest.

#### **Appendix A**

The scripts that were created to estimate the geometry of a river in a more or less automated manner are presented below. The scripts are hosted on the GEE servers and are openly available for anyone to read. Comments are limited to ensure the clarity of the code https://code.earthengine.google. com/f4e3e8ff6aab575a4a83a54f6cd0d93f.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Monitoring the Water Quality of Small Water Bodies Using High-Resolution Remote Sensing Data**

#### **Zehra Yigit Avdan 1,\*, Gordana Kaplan <sup>2</sup> , Serdar Goncu <sup>1</sup> and Ugur Avdan <sup>2</sup>**


Received: 4 November 2019; Accepted: 1 December 2019; Published: 2 December 2019

**Abstract:** Remotely sensed data can reinforce the abilities of water resources researchers and decision-makers to monitor water quality more effectively. In the past few decades, remote sensing techniques have been widely used to measure qualitative water quality parameters. However, the use of moderate resolution sensors may not meet the requirements for monitoring small water bodies. Water quality in a small dam was assessed using high-resolution satellite data from RapidEye and in situ measurements collected a few days apart. The satellite carries a five-band multispectral optical imager with a ground sampling distance of 5 m at its nadir and a swath width of 80 km. Several different algorithms were evaluated using Pearson correlation coefficients for electrical conductivity (EC), total dissolved soils (TDS), water transparency, water turbidity, depth, suspended particular matter (SPM), and chlorophyll-a. The results indicate strong correlation between the investigated parameters and RapidEye reflectance, especially in the red and red-edge portion with highest correlation between red-edge band and water turbidity (r<sup>2</sup> = 0.92). Two of the investigated indices showed good correlation in almost all of the water quality parameters with correlation higher than 0.80. The findings of this study emphasize the use of both high-resolution remote sensing imagery and red-edge portion of the electromagnetic spectrum for monitoring several water quality parameters in small water areas.

**Keywords:** RapidEye; water quality; red edge; remote sensing

#### **1. Introduction**

Water quality monitoring in surface water is important, in order to obtain quantitative information on the waters characteristics. A challenge to measuring and monitoring water quality in situ is that it can be prohibitively expensive and time consuming. An alternative method to monitor water quality is by using satellite imagery. In the past few decades, remote sensing techniques and capabilities have been studied for monitoring several water quality parameters. Remote sensing is an effective tool for synoptic soil moisture assessment [1], water and water level monitoring [2], water demand modeling [3], groundwater management, flood mapping [4], and water quality monitoring [5]. Depending on the study area, different satellite sensors have been used for water bodies monitoring. Thus, moderate resolution data from Sentinel-1 has been used for water resource management applications [6], Landsat-8 [7,8], and Sentinel-2 [9–12] have been used for water bodies extraction and water quality monitoring, MODIS (Moderate Resolution Imaging Spectroradiometer) has been used for water quality assessment [13]. Also, unmanned aerial vehicle (UAV) data have been used for water quality measurements [14]. Landsat data have also been used for assessing detailed water quality parameters, such as suspended sediments, water transparency, chlorophyll-a, and turbidity [15]. Water quality parameters help in decision-making regarding the further use of water based on its quality. However, the use of middle or low spatial resolution data (10–30 m) is not useful in monitoring small lakes, dams or rivers.

The main indicators of water quality are chlorophyll-a, water transparency, turbidity, and suspended particulate matter (SPM) [16]. Chlorophyll-a is a commonly used measure of water quality as eutrophication level. Higher concentrations indicate poor water quality, usually when high algal production is maintained due to high chlorophyll-a and nutrient concentrations. Water transparency can also be used as an indicator of water quality. It is measured with the simple Secchi disk depth. Turbidity is the cloudiness or haziness of the water caused by large numbers of individual particles that are generally invisible to the naked eye. The measurement of turbidity is a key test of water quality. The measurement of total dissolved soils (TDS) is also important as it shows the concentration of dissolved solid particles in the water.

The use of satellite data to monitor water quality parameters has been mainly focused on developing algorithms with band combinations. The remote sensing algorithms for water quality assessment generally involves visible and infrared portion of the electromagnetic spectrum. Thus, turbidity, as one of the most important water quality parameters, has been monitored through remote sensing and several remote sensing studies have quantified the relationship between actual and estimated turbidity from satellite remote sensing data [17,18]. Similar investigations have been made for chlorophyll-a [19], SPM [20], water transparency [21], etc. In order to establish the reliability of the remote sensing data, the results are compared with in situ measurements.

Building on this previous work, this study elaborates on monitoring water quality of small water bodies using high-resolution (5 m spatial resolution) data from RapidEye and implements indices widely used in the literature. However, since there are number of indices, in this study we only used indices that gave significant results or indices that needed clarification about their use. That is the case for one of the indices, that was used in the literature under different names, and was used for both turbidity and chlorophyll-a estimation.

Thus, in this study, we use eight different indices retrieved from high spatial resolution satellite sensors for estimating the correlation between the in situ and estimated water quality parameters: pH, electrical conductivity (EC), TDS, water transparency, water turbidity, depth, SPM, and chlorophyll-a. Specific objectives were: a) to retrieve water quality parameters from RapidEye reflectance, and b) to correlate satellite retrieved parameters with in situ measurements. The paper has been structured as follows: The materials and methods sections provide information about the study area, in situ and laboratory measurements, satellite data, and methodology. In the third part of the paper, the results are presented, followed by detailed discussion and conclusion.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Borabey Dam, built in the period of 1991–1992, is located in Emirce Village on Bozdag slopes north of the Eskisehir city center (Figure 1). The lake covers an area of approximately 166,559 m<sup>2</sup> at an altitude of approximately 924 m. Borabey Dam was built in order to support the local agriculture. Firstly, the lake was used as the Water Sports Center of Anadolu University in 1999. Later, it was planned to contribute to the drinking and utility water network of Eskisehir. The lake is mainly fed through rain and snow-melting, and one main stream. However, in 2011, this purpose was canceled. Details about Borabey Lake are given in Table 1 [22].

Eskisehir, which is located in the Central Anatolian region, generally has low rainfall. The annual rainfall varies between 300–500 mm. When it comes to the distribution of precipitation according to the seasons, the summer season is generally dry, autumn is rainy, while in spring and winter, there are relatively heavy rains. In winter, a significant proportion of precipitation falls in the form of snow. The climate is continental with very high temperatures during the summer days, but with cool nights, and cold and hard winters.

**Figure 1.** Study Area, Borabey Lake (left; RapidEye image, RGB - 3,5,1), Eskisehir, Turkey.


**Table 1.** Borabey Dam Technical Information.

#### *2.2. In situ and Laboratory Measurements*

Water sample collections were initiated during the summer season, on August 12, 2014. Following a random sampling technique, water samples were collected and transferred to the laboratory to perform further analyses. In total, twenty-one water samples were collected from the Borabey Lake, and for each sample, the coordinates were recorded with a Global Positioning System (GPS) (Universal Transverse Mercator - UTM coordinates (World Geodetic System – WGS 84) Zone 36S). The samples were collected from the water surface. All of the samples were collected in glass bottles, while the samples used for the Chl-a analyses were transported in opaque bottles. For each sample, the following parameters were estimated: pH; electrical conductivity (EC); total dissolved solids (TDS); Water transparency (Secchi Disk Depth—SDD); depth; turbidity; suspended particulate matter (SPM); and chlorophyll-a (Chl-a). For all of the parameters, two measurements were made and then averaged. Since there was no significant difference between the two readings of the Chl-a, only the average values were given.

pH, EC and TDS values were measured with portable multimeter (Hach HQ40d). Turbidity was measured with a portable field type turbidity meter—(WTW—Turb®355 IR) as NTU unit. A black/white stripe and 25-cm diameter Secchi disk was used for measuring SDD values. Especially, midday hours are preferred for the measuring SDD values. Chl-a values are measured with using Standard Methods (10200 H) analyzed with a spectrophotometer; SPM values were measured with a gravimetric method at the laboratory. In order to determine the SPM amount, the water samples were filtered through a membrane filter (0.45 μm pore size) and Sartorious filtration apparatus was used.

#### *2.3. Satellite Data*

One RapidEye product acquired three days after the in situ measurements, August 17, 2014, was used in this study. The RapidEye instrument acquires data in five different spectral bands (Table 2) with 5 m spatial and approximately six days temporal resolution. The products included a standard geometric correction resampled to a 5 × 5 m pixel size. The data were collected from the Planet Explorer webpage [23]. The RapidEye image with the in situ measurements are given in Figure 2.

**Figure 2.** RapidEye-3 image of Borabey Lake, and locations of in situ measurements.


**Table 2.** Spectral bands for the RapidEye Satellite Constellation.

#### *2.4. Methods*

In situ measurements were used to evaluate eight different remote sensing indices that were successfully used and proposed in the literature for different remote sensing sensors. For simplicity, the indices have been renamed as I1-I8 (Table 3). The structure of the first four indices is similar to the widely used Normalized Difference Vegetation Index (NDVI). The first index I1, Normalized

Difference Water Index (NDWI) [24], has been widely used for water bodies extraction. Even though the index has been developed to delineate open water features, it is also assumed that it can provide turbidity estimations of water bodies. The index is calculated from the green and near-infrared bands, which are highly sensitive to water contents. I2 is similar to the widely used Normalized Difference Chlorophyll Index (NDCI) [25], but instead of red, uses the red-edge band. I3, or Normalized Difference Turbidity Index NDTI [26,27], uses red and green reflectances for estimating the turbidity in water bodies. However, the same index has been used for chlorophyll-a estimation [25,28,29]. I4 was formed from the red and red-edge portion of the electromagnetic spectrum [29]. I5-I7 are two band ratios, successfully used in the literature for water quality assessment. Namely, the green–red ratio, I5, and the red–red edge ratio have been successfully used for water quality monitoring, especially in chlorophyll-a estimation [30,31]. Similar to I5, I7 is a red–green ratio, used for total suspended matter (TSM) estimation [31]. I8 is a three-band index used for turbidity estimation [32]. The indices were then adjusted according to the RapidEye satellite sensor bands. Afterward, the results from the retrieved indices were statistically analyzed against the water quality parameters from the in situ and laboratory measurements. In addition, correlation analyses were made between every spectral band from RapidEye and every water quality parameters used in this study, and between the in situ and laboratory measurements of the water quality parameters.


**Table 3.** List of indices tested in this study.

#### **3. Results**

#### *3.1. In Situ Data*

The results of the in situ and laboratory measurements are summarized in Table 4. The minimum, maximum and mean values of the eight studied water quality parameters of all 21 measurements collected on August 12, 2014 are presented. In addition, a polynomial relationship was obtained between in situ measurements of SDD and turbidity (r2 = 0.79) and in situ measurements of SDD and laboratory measurement of chlorophyll-a (r<sup>2</sup> = 0.51). The results are shown in Figure 3.

**Table 4.** In situ data of the investigated water quality parameters collected on 12 August 2014.


**Figure 3.** Relationship between: (**a**) SDD and turbidity; (**b**) SDD and chlorophyll-a.

#### *3.2. RapidEye Data vs. In Situ and Laboratory Data results*

In situ and laboratory water quality parameters were collected five days before the satellite image. Band by band correlation was calculated between the RapidEye data and the in situ/laboratory measurements. Also, significance levels for every correlation were calculated. The results are given in Table 5. The results of pH were not presented, as there was no correlation with the satellite data.


**Table 5.** Pearson's correlation coefficient between RapidEye reflectance and in situ and laboratory measurements.

Although high correlations were noticed in the near infra-red (NIR) part, all of the water quality parameters gave highest correlation with the red and red-edge parts of the spectrum. The highest correlation was noticed between the in situ turbidity and the red-edge band with R = 0.9, and p = 0.00. The chlorophyll-a had lowest with R = 0.39 in the red-edge portion. However, after analyzing the results, it was concluded that the first measurement was divergent from the others, thus, after the first measurement was excluded from the statistical analyses, the correlation was R = 0.53, and p = 0.02. The results are shown in Figure 4. For better visualization, normalization of the data was performed (0–1).

**Figure 4.** Relationship between red edge wavelength and in situ turbidity (blue); and chlorophyll-a (red).

The in situ and laboratory measurement values were also regressed against the different band combinations (Table 3) successfully used in the literature for different sensors. The correlations between the indices retrieved from RapidEye data and the in situ and laboratory measurements were calculated. Also, significance levels for every correlation were calculated. The results are given in Table 6 and some significant results are shown in Figure 5.


**Table 6.** Correlation (R) RapidEye indices and in situ measurements.

**Figure 5.** Correlation between estimated and actual water quality parameters (I1-I8 are shown in Table 3).

The EC has the highest correlation with the I8 index, index that is generally used for estimating water turbidity, and I7, index used for TSM estimation. While the first index is a combination of RGB bands, the second index is a red–green ratio, indicating that EC is sensitive in the visible spectra. It should also be noted that EC showed high second-degree polynomial correlation with I8 of r2 = 0.81. Similar is the case with TDS, where the highest correlation is noticed in the I8 index of R = –0.49. All indices other than I4 gave significant results in the water transparency estimation (SDD) with the highest negative correlation of R = –0.86 for the I3 and I7 indices. I3 and I5 were successful in estimation of the depth of the water with correlations higher than R = 0.8. Also, the same indices had good second-degree polynomial correlation higher than r2 = 0.82. Turbidity was best estimated using I3 and I7 with correlation higher than R = 0.9. However, the other indices other than I4 and I6 gave significant correlations higher than 0.84. The highest coefficient of determination was achieved with second degree polynomial correlation of r2 = 0.88 between actual turbidity and I8. The green–red and red–green ratios were significant for SPM estimation. As expected, chlorophyll-a was best estimated using I2 (R = −0.61) and I5 (R = –0.54). I2 is a combination of green and red-edge portion, while I5 is a

green–red ratio, widely used in the literature for chlorophyll-a estimation. I3 also gave satisfactory results, with R = 0.53. Although in the literature I3 is defined as turbidity index, it has also been used as a chlorophyll-a index. Also, I2 gave highest linear correlation of determination of r2 = 0.37.

#### **4. Discussion**

The main objective of the represented study was to investigate the ability of high spatial resolution satellite data for estimating water quality in small water areas. For that purpose, several successful algorithms that have been used by many authors to examine characteristic relationships between remotely sensed data and water quality parameters were selected and applied to high-resolution imagery. In this study, we investigated the possibility of using RapidEye products to monitor water quality parameters in Borabey Lake, Eskisehir, Turkey. Although it would be more suitable if the in situ and satellite data were acquired on the same day, since there was no significant hydrological event (e.g. rain) between 12 and 17 August 2014, and the main stream was dry during the summer season (July–August), it was assumed that the water quality did not changed during the short period of time. The study includes the elaboration of regression analysis to examine the relations between the reflectance of RapidEye bands and the water quality parameters. The parameters used in this study are: EC, TDS, water transparency, depth, turbidity, SPM, and chlorophyll-a. In situ measurements of the mentioned parameters were collected few days before the acquired satellite image. Since there are number of algorithms/indices in the literature for water quality monitoring, only the most successful were selected in this study. Also, it should be noticed that some of the indices have been used for monitoring different water quality parameters, such as the I3, that has been used as both chlorophyll-a and turbidity monitoring. For simplicity, the indices used in this study have been renamed as I1-I8 (Table 3).

Several studies have explored the possibility of determining water quality with remote sensing images. Many of these studies consider Landsat imagery (30 m spatial resolution) to be high-resolution [16], which may not be satisfactory for monitoring small water bodies. The comparison between two different remote sensing sensors indicates that higher spatial resolution offers higher accuracy in water quality parameters [34,35]. Taking into consideration that results of water quality are limited due to mixed pixels [36], in this study we evaluate high-resolution (5 m spatial resolution) for monitoring water quality in small water bodies. To our knowledge, only few studies can be found in the literature exploring high remote sensing imagery (1–5 m) for water quality parameters, namely exploring only single parameter [37].

The main results of the study showed a significant relation between RapidEye bands and in situ and laboratory measurements. Similar to different studies, it was concluded that the most useful spectral information, needed to retrieve different water quality parameters are the red and NIR part of the spectrum [38,39]. However, our study also showed that red-edge part of the spectrum is more useful than the NIR for some of the parameters. For example, while EC and TDS showed higher correlation with the NIR band, all of the other parameters showed significant higher correlation with the red-edge band. The 690–730 wavelength range has been widely used in water quality monitoring [31], supporting our findings in Table 5.

As expected, I1 or NDWI, gave significant results in the water quality parameters. Although it was strongly suggested that it can be used for both water bodies extraction and turbidity monitoring [24], the results in this study indicate that turbidity can be better monitored through a combination of red and green bands. Taking into consideration the higher correlation of red-edge bands with chlorophyll-a rather than NIR, similar to NDCI, I2 was constructed from green and red-edge bands. I5, a green–red ratio, which is also widely used for chlorophyll-a monitoring, gave lower correlation than I2. Also, I2 performed best with a linear correlation of R2 = 0.37. In the literature, I3 (a combination of red and green bands) has been used for both chlorophyll-a and turbidity estimation. The results in this study showed that I3 gave second-best result in the turbidity estimation among all of the investigated indices (R = 0.92), while in the chlorophyll-a estimation was ranked as third (R = 0.54).

The results in this study are slightly better than similar investigations; for example, Masocha et al. [8] found that blue–red ratio provided strong positive relation between measured and retrieved turbidity in two different lakes (R = 0.90; R = 0.80). Different turbidity retrieval approaches showed a good correlation for two different study areas (R = 0.81) [40].

The findings emphasize the use of both high-resolution remote sensing imagery and the red-edge portion of the electromagnetic spectrum for monitoring several water quality parameters in small water areas. Also, this study represents another case study that confirms the use of satellite remote sensing in water quality mapping and monitoring.

#### **5. Conclusions**

The presented study assessed the ability of high-resolution satellite sensor to monitor several water quality parameters in small water bodies, such as Borabey Lake in Eskisehir, Turkey. The relation between the in situ measurements, as well as the relation between each RapidEye band, and indices retrieved from RapidEye data were investigated. The main findings of this study proved the ability of RapidEye sensor in water quality estimation with high correlation between the in situ data and the RapidEye reflectance. The results showed a high sensitivity of the water quality parameters in the red, red-edge, and NIR part of the spectrum. The highest correlation has been noticed between in situ turbidity and the red-edge band (r2 = 0.92). Also, among the investigated water quality indices, the highest correlation was noticed for turbidity estimation (r<sup>2</sup> = 0.88).

Investigating different indices, widely used in the literature, it was concluded that some are more successful than others. Thus, although some of the indices were used for both chlorophyll-a and turbidly estimation (I3), the finding of this study showed that chlorophyll-a can be best observed with combination of green and red-edge bands, while good turbidity estimation can be obtained with combination of red and green portions of the spectrum (I3 and I6). As expected, the same combination can be used for estimating water transparency. No significant correlation between red and red-edge bands was observed.

In order to improve the results of this study, a higher number of in situ measurements may be needed. Also, different water quality parameters with high-resolution satellite imagery should be observed. As the algorithms used in this study are easy to implement, we recommend further analyses in different study areas at different seasons of the year in order to get a wider range of values of water quality. Similar studies should be conducted in both turbid and clear water in order to support the results and encourage remote sensing data in water quality monitoring.

**Author Contributions:** Zehra Yigit Avdan and Serdar Goncu did the in situ measurements and water quality assessments. Gordana Kaplan and Ugur Avdan did the remote sensing data processing and analyses. All author contributed equally in the writing of the paper.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors want to thank Özge Bilget for her contribution during the data collection and Planet Labs, Inc. for providing the RapidEye imagery.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Estimation of Crop Water Deficit in Lower Bari Doab, Pakistan Using Reflection-Based Crop Coe**ffi**cient**

#### **Muhammad Asif Javed \*, Sajid Rashid Ahmad, Wakas Karim Awan and Bilal Ahmed Munir**

College of Earth and Environmental Sciences, University of the Punjab, Lahore 54000, Pakistan; principal.cees@pu.edu.pk (S.R.A.); waqaskareemawan@gmail.com (W.K.A.); bilalahmedmunir35@gmail.com (B.A.M.)

**\*** Correspondence: m.asifjaved13@gmail.com; Tel.: +923214700241

Received: 24 January 2020; Accepted: 1 March 2020; Published: 13 March 2020

**Abstract:** There is a global realization in all governmental setups of the need to provoke the efficient appraisal of crop water budgeting in order to manage water resources efficiently. This study aims to use the satellite remote sensing techniques to determine the water deficit in the crop rich Lower Bari Doab Canal (LBDC) command area. Crop classification was performed using multi-temporal NDVI profiles of Landsat-8 imagery by distinguishing the crop cycles based on reflectance curves. The reflectance-based crop coefficients (Kc) were derived by linear regression between normalized difference vegetation index (NDVI) cycles of the Moderate Resolution Imaging Spectroradiometer (MODIS) MOD13Q1 and MYD13Q1 products and Food and Agriculture Organization (FAO) defined crop coefficients. A MODIS 250 m NDVI product of the last 10 years (2004-2013) was used to identify the best performing crop cycle using Fourier filter method. The meteorological parameters including rainfall and temperature substantiated the reference evapotranspiration (ET0) calculated using the Hargreaves method. The difference of potential ET and actual ET, derived from the reflectance-based Kc calculated using reference NDVI and current NDVI, generates the water deficit. Results depict the strong correlation between ET, temperature and rainfall, as the regions having maximum temperature resulted in high ET and low rainfall and vice versa. The derived Kc values were observed to be accurate when compared with the crop calendar. Results revealed maximum water deficit at middle stage of the crops, which were observed to be particularly higher at the tail of the canal command. Moreover, results also depicted that kharif (summer) crops suffer higher deficit in comparison to rabi (winter) crops due to higher ET demand caused by higher temperature. Results of the research can be utilized for rational allocation of canal supplies and guiding farmers towards usage of alternate sources to avoid crop water stress.

**Keywords:** crop water requirement; reflectance-based crop coefficients; normalized difference vegetation index; evapotranspiration

#### **1. Introduction**

The irrigated agriculture sector is the prime user of freshwater resources around the world and consumes approximately 69% of the freshwater withdrawal [1]. Asia has the largest consumption of around 56% of global fresh water for irrigation purposes [2]. Freshwater is becoming increasingly scarce in general and in Asia in particular [3]. The existing canal irrigation system in Pakistan is a supply-based system working on the principle of equitable distribution independent of the crop cover [4]. Major constraints faced by irrigated agriculture in Pakistan are canal water scarcity and its uneven distribution. Irregular rainfall pattern and water logging also have severe impact on the crop productivity [5,6]. The water resources of Pakistan, both groundwater and surface water, have become inadequate to fulfill the growing demands of the irrigation-based agriculture sector [7]. Canal water

is not sufficient to solely satisfy the crop water needs as it fulfills only 37.5% of the crop demand in Punjab. Rainfall contributes 15.5% and groundwater contributes 47% of the remaining demand. Crops receive only 70% of their demand in Punjab despite consumptive use of all sources including rainfall [8]. Punjab province, the largest irrigation area of Pakistan, is an intensely cultivated region covering an area of about 8.4 million ha [9]. This area observes cultivation throughout the year in two crop seasons namely summer (kharif) and winter (rabi). The high flow irrigation period falls into the summer from June to August, while late winter season from February up to early June is observed as low flow season [10]. A rotational 8-day irrigation scheme is implemented to equally satisfy the water need in all irrigation districts [11]. There is a realization globally of the need to efficiently utilize water resources to avoid crop water stress and increase crop productivity [12].

Crop water consumption, generally known as evapotranspiration (ET), can be calculated accurately at pixel level using remote sensing technology [13]. ET is the key parameter used for improvement in agriculture water management. Several methods are being used for direct or indirect measurement of crop ET; these include weighing lysimeter-based measurements, eddy covariance method, Bowen ration surface energy balance (BREB), surface energy balance algorithm for land (SEBAL), and using reference ET by multiplying it with crop coefficients [14]. Estimation of reference evapotranspiration based on crop coefficients is being used commonly for farm level agriculture management [6]. Reference evapotranspiration (ETo) is mostly calculated using weather parameters by common equations known as the Hargreaves equation and the Penman–Monteith equation. The Hargreaves equation requires minimum and maximum temperature as input for ETo calculation, while the Penman–Monteith equation requires temperature, wind speed, humidity, and sun radiation as input weather parameters [15]. ETo calculated through Penman–Monteith or Hargreaves equations reflect the values of a reference crop in ideally suited climate conditions and needs to be converted to crop specific consumption, i.e., actual evapotranspiration (ETa) using crop coefficients [9]. Crop coefficients reflect variable water need as per stage of the crop from sowing up to harvesting and vary for each crop as per its leaf area density and climate conditions [16].

The crop coefficient approach is considered as reliable and efficient in determination of the spatiotemporal growth of crop water needs and crop specific variations show strong correlations with the satellite derived spectral index, normalized difference vegetation index (NDVI) [6,17–19]. Through regression calculation between crop coefficients and the NDVI, a set of crop specific parameters is created that allows the calculation of reflectance crop coefficients (Kcr) [17]. By applying Kcr to both, an ideal reference crop and the monitored crop, the difference between these two can be translated into crop water stress. This approach is quite significant as it provides crop coefficients (Kc) for each stage of the crop rather than providing a static estimate of Kc found in published reports or manuals. This method includes the data of land use for different dates throughout the crop growing season derived from the freely available satellite images of the whole year. This method is considered cost-effective in terms of data requirements, i.e., satellite images are available freely and climate data are available online in open access in comparison to the data collected through field surveys [20].

The spectral crop coefficient approach was used to monitor crop water consumption and to test the accuracy of crop coefficients for three fields in the Texas High Plains, all planted with cotton. Comparison was made of crop water use calculated from the reflectance-based crop coefficients approach and eddy covariance measurement approach (measured from a flux tower). In investigations, researchers used Landsat data for NDVI calculation, and ETo was calculated after the Penman-Monteith equation. They concluded that the use of reflectance crop coefficients is effective under a variety of different irrigation conditions and also stressed the adaptive character of the approach accounting for specific field conditions as opposed to standard crop coefficients [21]. NDVI-based Kc calculation approach is significant in crop water requirement estimation. The crop specific water requirements are expressed in terms of Kc and have been quantified for different crops. Water requirements can then be found through multiplying crop coefficients times the reference evapotranspiration [22–26]. Using modeled potential evapotranspiration (PET) and the relationship between NDVI and Kc,

NDVI differences between ideal and monitored crop can be translated into irrigation water needs. Successful applications of this technique or approach have been reported in many studies and projects [25–28].

The purpose of this study is to provide the spatial distribution pattern of crop water requirement by monitoring the crop health in the Lower Bari Doab Canal (LBDC) command of the Punjab province. The results from this study would provide useful information for allocating water in main and secondary canals in order to avoid crop water stress and provide equitable water across the canal command. In summary, the objectives of this study are: (i) identification and quantification of different cropping patterns in LBDC; (ii) developing a decision support model for crop water efficacy at monthly intervals by using hydro-meteorological, geographical, and solar parameters; and (iii) scrutinizing the potential water deficit along the canal by counting upon the contribution of irrigation supplies and rainfall for both rabi and kharif seasons.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area comprises four irrigation divisions, i.e., Balloki, Okara, Sahiwal, and Kahnewal, which are located in the longitudinal direction from the LBDC head. The climate of the study area is semi-arid. Mean daily maximum and minimum temperature varies between 32–36 ◦C and 11–15 ◦C, respectively. Temperature differs along the length of the canal command and gets warmer from east towards west. The mean yearly rainfall varies progressively from around 350 mm in the east to 200 mm in the west, the greater part of which falls in the storm, i.e., in the months of July and August. The winter precipitation generally occurs in the months of January and February in patches because of cyclonic storms from the south west. The winter season is gentle to chilly and stretches out from November to March [29]. Figure 1 describes the historical patterns of temperature and rainfall in study area for the years 2003–2012. The general slope of the territory is mild towards the south-west bearing, while the normal slope reaches from 1:4000 to 1:10,000. The LBDC canal command, overwhelmingly agricultural, lies at an altitude of 394–640 feet over the mean sea level. Texturally, the soil of the region is 70% medium (loam and sediment loam), 20% tolerably coarse (sandy soil and fine-sandy loam), 4% coarse (sandy and sandy loam), and 3% respectably fine (sandy-earth soil and silty-earth loam). The remaining 3% region possesses miscellaneous land types. The soils of the zone are intrinsically low in natural (organic) matter and accessible phosphorus [2,29].

LBDC originates from the Balloki headworks located southwest of Lahore, approximately 65 km apart, and ends in Khanewal division near Jahanian Tehsil as shown in Figure 2. A cultivable command area of about 1.7 million acres is being served by LBDC falling in the districts of Kasur, Okara, Sahiwal, and Khanewal. The design discharge of the canal is 9841 cusecs. The present carrying capacity of the canal is 8600 cusecs in the high flow season (June to August). The canal is usually not run on its designed capacity because of the shortage of water and also due to safety measures keeping in mind the old canal infrastructure [30]. Branch and distributary canals associated with the main LBDC are usually operated on a rotation basis to ensure equitable distribution of canal water. The rotation plan is prepared by each irrigation division at start of every cropping season by making groups of the canals based on their discharge, and priorities are fixed to release water on 8 daily bases. Generally, 3 groups are made of more or less equal discharge to rotate canals within a division, and each group gets its full share at least once every 24 days. The group falling in first priority gets full discharge, and remaining groups get a balanced share as per second or third priority fixed in the rotation plan [10]. There are approximately 275,000 farm families in this command area deriving their livelihoods directly from crop cultivation including wheat, maize, cotton, rice, flowers, vegetables, sugarcane, fodder, orchards, and citrus crops. The LBDC command area is considered as progressive and an important agriculture area in Punjab province and provides substantial potential for increased productivity. Wheat has the

major share amongst all the crops in study area, while cotton and maize have the second and third largest share, respectively [30].

**Figure 1.** Time series graph of temperature and precipitation.

**Figure 2.** Study area map.

#### *2.2. Data Collection*

Multiple datasets were collected from different sources to be used to demonstrate the crop water deficit in the LBDC area. Hydrometeorology, irrigation, and remote sensing satellite information collected from different sources including organizational and open-source platforms were pre-processed in a GIS environment (Table 1). However, methodological flowchart of the research is shown in Figure 3.


**Figure 3.** Methodological framework of the study.

#### *2.3. Data Processing*

The methodology used the concept of reference evapotranspiration (*ETo*) at a spatiotemporal scale. The empirical approach of the Hargreaves method has been used for the estimation of *ETo* in the LBDC command area. The Food and Agriculture Organization (FAO) of the United Nations proposed the Penman-Monteith equation, Equation (1), as the standard method for estimating *ETo*. The Penman-Monteith equation requires multiple climate parameters, which were not available with weather stations existing in the study area on the required frequency that hence could have yielded wide gaps and/or inaccuracies in *ETo* estimation. To solve this problem, the Hargreaves equation (HARG) is recognized by FAO and is often used in data poor areas. The Hargreaves equation is based on average, minimum, and maximum air temperature and extraterrestrial radiation. It tends to overestimate *ETo* in humid conditions and requires a local calibration. Moreover, this approach is simple and effective and can be used for daily seasonal and monthly calculations of *ETo* [31–33].

Different schemes have been employed for calculation of *ETo*, and the Hargreaves method was observed to be simple and highly correlated for seasonal and daily bases [34]. The mathematical equation of *ETo* calculation is given as:

$$ET\_o = 0.0023 \times \left(\frac{T\_{\text{max}} + T\_{\text{min}}}{2} + 17.8\right) \times \sqrt{T\_{\text{max}} - T\_{\text{min}}} \tag{1}$$

where, T-max(◦C) and T-min(◦C) are average daily maximum and minimum temperatures, respectively, and R is extra-terrestrial solar radiation.

Daily temperature data available for weather stations located inside the LBDC region were used to calculate the *ETo*. To get the *ETo* monthly measurements, monthly mean temperatures were supposed to resemble that of an average monthly day, so Equation (1) was applied, and then result was multiplied by the total number of days of the corresponding month. Afterwards, for spatial interpolation of *ETo*, inverse distance weight (IDW) was chosen over other complex methods (e.g., cokriging) that gave unsatisfactory results. This operation provides the spatial distribution of *ETo* at monthly intervals [35].

However, among the meteorological information, rainfall data were collected from the open-source platform of NASA (National Aeronautics and Space Administration: Earth Data). GPM (Global Precipitation Mission) and TRMM (Tropical Rainfall Measurement Mission) datasets were collected at daily temporal resolution and processed to monthly information. All information was standardized at the same scale and was restricted to the irrigation division level of LBDC command area. The datasets were pre-processed and resampled, if needed, for monthly calculation of crop water requirements for winter (rabi) and summer (kharif) seasons.

#### *2.4. Crop Classification*

For accurate quantification of water requirement and consumption, crops were classified into two seasons: rabi (winter) crops and kharif (summer) crops. This required a differentiation of major crop types grown in a calendar year. The crop pattern in the study area was identified using Landsat satellite imagery. The product was used for its free availability and reasonable spatial, spectral, and temporal resolution. The time series NDVI is widely considered as a reflection of crop phenology [36]. NDVI time series were calculated from near infrared and red bands of Landsat, and the phonological shapes were used for crop identification. To cover the LBDC area, two Landsat footprints/tiles (Path 149 Row 39 and Path 150 Row 39 as shown in Figure 4) were processed under the object-based classification algorithm (OBCA).

**Figure 4.** Landsat tiles coverage of the study area.

A small part of the Balloki division was lying in the 3rd tile of Landsat (Path 149 Row 38); however, it was stitched with the corresponding tile (Path 149 Row 39). For the composition of NDVI time series, each footprint was processed using 21 Landsat time series images of years 2013–2014.

The stacked NDVI of the study area was used for pattern generation representing the crop phenology throughout the year. Each crop behaves differently in NDVI representing cycles. The clusters of different vegetation covers are grouped in the remote sensing environment to generate the crop pattern in the area. The classification technique using NDVI series is used by different researchers and has been found efficient for crop type mapping and monitoring [6,18,19,37]. Specimens of NDVI-based phonological profiles of wheat and cotton crops are shown in Figures 5 and 6, respectively.

**Figure 5.** Phonological profile of wheat crop.

**Figure 6.** Phonological profile of cotton crop.

#### *2.5. MODIS, NDVI and Reference Crop Cycle*

NDVI is an index that evaluates biomass calculated from a near infrared band and a red band. The NDVI product of MODIS (AQUA/TERRA) is used for reference crop cycle generation. The reference crop cycle defines the maximum value of crop reflectance observed in past 10 years. These cycles of individual crops behave near ideally throughout the season from crop sowing to harvesting stage. To meet the optimum condition, it was of the utmost importance to minimize the data gaps. MODIS provides NDVI from two platforms: AQUA and TERA. These platforms deliver NDVI imagery with an offset of 16 days, however, harmonizing the two dataset results in 8 daily by-products of NDVI. These NDVI images are stacked over a period of 10 years for reference cycle identification using Fourier filter technique.

#### *2.6. Crop Coe*ffi*cients (Kc) and Crop Water Requirement*

Crop coefficients have a strong correlation with satellite derived *NDVI* values, and both are directly proportional to each other. Reflectance-based crop coefficients (*Kcr*) are calculated through linear regression between *NDVI* and FAO defined crop coefficients (Kc) values [6]. The shaded grey area in Figure 7a represents *NDVI* values included in the linear regression of sugarcane crop, and Figure 7b represents corresponding FAO defined Kc values. After applying the linear regression, the resultant coefficients named *C*1 and *C*2 were determined for calculation of reflected crop coefficients (*Kcr*) as shown in Figure 7c. The linear equation 2 used for *Kcr* calculation for both the monitoring period and reference crop period is given below.

$$\text{Kcr} = \text{Cl} \times \text{NDVI} + \text{Cl} \tag{2}$$

The *C*1 and *C*2 coefficients are crop specific and were determined individually. However, there exist generally applicable *C*1 and *C*2 variables also [6]. The reflectance-based *Kcr* were generated for reference (ideal) crop cycle and actual NDVI collected at 8 daily intervals. *Kcr* of both actual and reference NDVI were then translated to crop water deficit calculation using the following Equation (3):

$$\text{CWR} = (\text{Kcrr} - \text{Kcra}) \times \text{ETo} \tag{3}$$

**Figure 7.** (**a**) NDVI profile of sugarcane crop, (**b**) Corresponding FAO defined crop coefficients (Kc values), and (**c**) Linear regression calculation between NDVI and Kc values.

#### **3. Results**

#### *3.1. Cropping Pattern in LBDC*

The vegetation index (NDVI) time series approach generates a total number of 12 land covers in the study area including six crop types. The classification scheme follows the crop calendar and defines summer (rabi) and winter (kharif) season crops separately.

The time series NDVI profiles of each class behave differently due to the intrinsic properties of crops. For example, wheat and rice are the perennial crops of rabi and kharif seasons substantiated in the study area, respectively. The NDVI profiles of both crops are closely related to each other, but the seasonal crop calendar (sowing and harvesting dates) is completely different. However, different weather conditions may lead to slight variations in the NDVI profile of same crops among different years. The results show that wheat, maize, and potato crops cover a major portion for rabi season in the LBDC region illustrated in Figure 8.

**Figure 8.** Crop classification-rabi season.

On the other hand, maize and cotton are the major crops of the kharif season as shown in Figure 9. The area comprises a combination of perennial and non-perennial crops as shown in Table 2.

**Figure 9.** Crop classification-kharif season.


**Table 2.** Crop types statistics in study area.

Spatial distribution shows that the northeast (NE) region of the LBDC command comprises maize; however, a huge concentration of cotton was found while moving towards the southern side of the area. The spatial distribution describing rice, fodder, and sugarcane classes covering small portions in the study area as shown in Figures 8 and 9.

#### *3.2. Crop Calendar and Coe*ffi*cients of LBDC*

The sowing and harvesting periods and stages are different for different spatial extents. However, the LBDC area comprises a large spatial extent and has variable stages of crops. The annual crop cycles for each crop in the study area were collected from the field survey and by considering the local expert knowledge. The crop calendar describes the temporal variation of individual crop growth over a year in the study area. Different stages of crops according to their sowing and harvesting time periods are presented in Table 3. The reflection-based crop coefficients (Kcr) were generated for individual crops using Equation (2). The average Kcr values of individual crops for years 2014–2016 are shown below in Figure 10.

**Table 3.** Crop calendar of different crops in the Lower Bari Doab Canal (LBDC) region.


**Figure 10.** Crop coefficients of various crops for a crop calendar.

The smaller NDVI values resulted in smaller Kc values and vice versa. The maximum value of Kc for each crop was observed between the middle to final stage whereas the initial stage reflected a low value. For instance, following the rice–wheat crop combination, the Kc value curve attains maximum value for wheat crop in the months of February (middle) and March (final stage). However, the same line goes to minimum in the months of May and June when wheat is harvested, and it again starts to climb when rice is cultivated. It again attains peak and goes to bottom after the harvesting of rice. The same is the case with other crop combinations corresponding to their sowing and harvesting periods. The case of sugarcane is different as it is considered a whole year crop, and its Kc curve attains peak value depending upon its sowing period.

The reflection-based Kc trend corresponded with the stages of crops (initial, middle, and final stages) tabulated in the crop calendar as shown in Table 3, which was prepared using crop phonological profile. Crop calendar was also verified during field survey by interviewing farmers.

#### *3.3. Spatial Distribution of Meteorological Parameters (ETo, T, and P)*

Temperature, evapotranspiration, and precipitation are very crucial parameters in calculation of crop water requirement. The parameters for the quantification of water deficit are used on a monthly basis for the years 2014, 2015, and 2016. Figure 11 shows the variation of evapotranspiration, temperature, and rainfall for the years 2014, 2015, and 2016.

**Figure 11.** Distribution of reference evapotranspiration (ETo), temperature, and precipitation for 2014, 2015, and 2016.

The results show that high values of evapotranspiration, temperature, and rainfall lie between the months April-August, April–September, and July–September, respectively, for the years 2014, 2015 and 2016. Maximum values of evapotranspiration and rainfall were observed in the months of June and July, respectively. These two parameters are in inverse relation to each other by indicating a decreasing trend of ETo compared to the rise in rainfall in the month of July. However, with the increase in temperature, the values of evapotranspiration also increased as shown in Figure 11, which indicates a positive relation between the two parameters.

Chronological behavior of evapotranspiration, temperature, and rainfall shows variability in the spatial distribution pattern in LBDC as shown in Figure 12. It has been observed that temperature and evapotranspiration have an increasing trend from east towards west in the study area. On the contrary, precipitation has a decreasing trend moving from west towards east. In general, similar pattern and

correlation between evapotranspiration and temperature parameters have been observed for all the months of three years in the study area.

**Figure 12.** Spatial distribution of meteorological parameters.

#### *3.4. Crop Water Deficit and Requirement*

The incorporated meteorological parameters and crop cycles resulted in chronological crop water deficit assessment in the LBDC region. Seasonal mean water deficit for the years 2014, 2015, and 2016 has been mapped in the study area as shown in Figure 13. It depicts the spatial distribution pattern of water deficit for rabi and kharif seasons in millimeters. Generally, the same spatial distribution pattern of water deficit for rabi and kharif seasons was observed in the study area. However, a deficit value in the kharif season is much higher than rabi season for all three years. This is because the kharif season lies in the hottest months of the year (i.e., June and July), which causes maximum evapotranspiration resulting in high water deficit.

**Figure 13.** Season-wise spatial distribution pattern of water deficit.

Figure 13 illustrates extremely high deficit observed at the tail (southwest) of the study area due to less irrigation supply because of maximum conveyance losses occurring and being located far away from the canal head. Similarly, southeast parts show high deficit for both seasons of these three years. These areas also lie along the border of the canal command and at the tail of the secondary canals, whereas northeast and northwest parts face considerably less deficit due to the advantage of lying close to the canal head and near to the Ravi river.

Irrigation supply, groundwater, and rainfall are the available sources to fulfill the crop water needs. However, groundwater pumping is not feasible or available for all farmers on demand, and small land holders usually rely on canal supply and rainfall. Influence and correlation of irrigation supply and rainfall parameters with crop water deficit measured on a monthly basis in LBDC are shown in Figure 14.

**Figure 14.** Crop water balance heat map for years 2014, 2015, and 2016 (developed in R-Studio).

The results show that less than 10 mm mean deficit is observed in the winter months (October-February) for the years 2014, 2015, and 2016. The low values of water deficits are mainly due to the lower evapotranspiration rate because of temperature falling down in these months. Specifically, December and January show relatively the least deficit due to the teleconnection of water deficit with evapotranspiration and temperature. Irrigation supply is not available in the months of January due to annual canal closure, but rainfall has supplemented to fill the gap. The crop calendar and trend graph of Kc also complement the deficit result because rabi crops are at an initial stage of their growth period during these months and require a substantially lower amount of water as shown in Figure 15.

**Figure 15.** Crop water balance heat map for year 2014 (developed in R-Studio).

The annual cycle of water budget for 2014 (Figure 15) shows that more than 40 mm of irrigation water was supplied to the LBDC throughout the year except for the months of January and February. This constant volume of canal supply is not sufficient for the crops in the peak summer months of July and August. Although aid of precipitation, around 45 mm, was available in the month of July, it was not able to meet the crop demand. However, it is obvious that when the peak summer months were over and temperature started to fall in the month of September, rainfall was also available in a good proportion of 50 mm, the deficit fell considerably. The deficit values dropped to a range of 5–10 mm from the month of October onwards because of the sufficient decrease in Kc values, hence low water requirement of crops.

In the water budget graph of the year 2015 (Figure 16), water deficit did not attain peaks in peak summer months of June to August, because the LBDC command area received high precipitation in these months. It was observed that a very high value of rainfall (>100 mm) for the month of July restricted the deficit to as low as 5 mm. However, the irrigation supply remained in the range of 40–45 mm for April-November. It is obvious from the water budget graph of the year 2015 that rainfall volume in the study area directly influences the crop water deficit. The water deficit pattern in the winter months is much more similar to the year 2014.

**Figure 16.** Water balance heat map for year 2015 (developed in R-Studio).

Similar rainfall and irrigation supply patterns were observed in the year 2016 as shown in Figure 17. However, rabi crops are observed with a shortage of water in the peak demand month of March due to a lower amount of rainfall as compared to years 2014 and 2015. High volume of rainfall in June and July for the kharif season kept the deficit under control, similar to the year 2015.

**Figure 17.** Crop water balance heat map for year 2016 (developed in R-Studio).

From results of all three years, the fact is established that the study area observes relatively high water deficit in the kharif season due to higher temperature, which causes maximum evapotranspiration in the months of June, July, and August. Although irrigation supply remains relatively high during these kharif months, crop demand is not satisfied unless sufficient volume of precipitation is available or groundwater is pumped in ample proportions. For the rabi season, the months of March and April are critical to meet crop water demand because of wheat being the dominant crop at its peak growth. Moreover, increase in ET values and decline in precipitation during peak demand months of rabi and kharif seasons favor water deficit in crops. On the contrary, the months of December, January, and February, being the coolest months of the year, face minimal ET deficit in the study area because crop water demand is fulfilled by the precipitation and irrigation supply and groundwater abstraction.

#### **4. Discussion**

Canal water scarcity and inefficient management of irrigation water channels are resulting in maximization of the use of groundwater to fulfill the crop water needs. Groundwater pumping is at maximum during the driest months of the year because of the higher rate of evapotranspiration. This fact is also confirmed in different published research that shows that in south Asia, Pakistan is the country with the highest groundwater usage of around 53% [38]. It is imperative to improve crop water allocation by identifying the particular cropland under stress in irrigation districts and disseminating it to canal managers and farmers in a timely manner. Visualizing the impact of spatial precipitation variations is necessary to rationalize canal water allocation, which is possible by establishing an approach that is capable of providing estimation of the spatiotemporal dissemination of crop water requirements [39]. Availability of spatial and temporal crop distribution data is one of the major components in devising policies for the agriculture sector regarding sustainable irrigation water management. A decision support tool for water resources management should be developed by acquiring the understanding of agricultural practices and the availability of reliable assessments of crop water requirements [40]. Particular challenges in the development of a crop water management system are the size of the Punjab irrigation area, spatial and temporal crop variability, the location specific stress situation of crops that cannot be compared to standard crops, and the developmental status of technically viable monitoring techniques/models. Spatial analyses at such scale and at such short repetitive steps are only feasible with remotely sensed images.

This pilot research study attempts to provide a cost- and time-effective approach for efficient monitoring of crop health over one of the largest canal command areas in Punjab, Pakistan. The major goal of this study was to develop a spatial decision support tool to integrate and standardize the various parameters for estimation of crop water needs. Multiple datasets including hydrometeorology, irrigation, and remote sensing satellite information collected from different sources including organizational and open-source platforms were employed in GIS and remote sensing environments. Reference evapotranspiration (ETo) and crop coefficients were of major concern for quantification of crop water requirements in the study area. The Hargreaves method, considered as simple and reliable, especially in data poor areas, has been employed for calculation of reference evapotranspiration (ETo) [35,41], while reflected crop coefficients (Kcr), which proved highly correlated with the NDVIs, were derived using MODIS NDVI product and FAO defined crop coefficients (Kc) for calculation of reference and actual crop coefficients [42]. Crop type and pattern were identified using Landsat imagery for accurate crop water quantification.

The research results verified the worth of satellite data for proficient appraisal of crop health monitoring. They also highlighted the worth of open-source information for crop water requirement estimation. The multi-step scheme for near real-time irrigation water stress estimation is quite effective in nature for large scale studies with limited data availability and therefore could be employed on other regions of the same agro-climatic conditions.

#### **5. Conclusions**

This study presents a technique for crop water deficit modeling in data poor areas and concluded the following:

(i) NDVI-based crop classification revealed that wheat and cotton are the major cultivated crops for rabi and kharif seasons with an area of approximately 60% and 30%, respectively. However, maize crop is non-perennial in nature, showing a variable cultivated area in rabi and kharif seasons. Sugarcane and rice, being the perennial and non-perennial crops, respectively, are also cultivated in a small portion of the study area. It is clear from the crop statistics that the kharif season has more variation in crop cover than the rabi season, hence requiring more dynamic water allocation.

(ii) The reflection-based Kc values remained consistent with crop calendar of LBDC generated by the local experts and published reports. The Kc values remained high at middle and final stages for both rabi and kharif season crops due to the maximum NDVI value, hence maximum crop water requirement at these stages.

(iii) The LBDC command area receives a continuous canal water supply throughout the year, except the month of January, which is observed as the annual canal closure. A relatively low volume of rainfall, i.e., 50 mm, was observed in July 2014, which caused the deficit to rise up to 30 mm in July, one of the hottest months, while the year 2015 received higher rainfall, especially in July, of up to 120 mm, minimizing the deficit to as low as 5 mm. However, the study area received 80 mm rainfall in July for the year 2016, and the deficit climbed up to 40 mm because of relatively low rainfall accompanied with a reduction in irrigation supplies due to maintenance work going on in LBDC [8].

(iv) The consumptive use of irrigation supply and groundwater is vital to meet crop water needs, though consumptive use is still not able to meet crop demand as observed in the results. However, the rainfall pattern directly influences the quantitative behavior of the deficit. The research results endorsed the water budget published by the On Farm Water Management (OFWM) Agriculture Department that stated that groundwater and irrigation supply are not sufficient to meet crop demand in the LBDC, and rainfall plays an important role especially in the summer to reduce deficit. Moreover, the deficit values remained high in March and April for rabi crops and in August and September for kharif crops due to the high crop water demand. Water deficit was observed to be consistent towards the tail end of the main canal and secondary canals, which needs to be addressed to ensure equitable distribution of irrigation water resources.

There is a dire need to use modern technology of remote sensing and GIS to efficiently utilize irrigation water resources and avoid crop water stress by monitoring crop health and climate conditions simultaneously. A spatial decision support system (SPSS) needs to be developed for better management of canal water and alternate sources available for crop water use. The current study has laid the foundations for using independent and open-source technology to monitor crop water stress on a near real-time basis and guide decision makers to take informed decisions.

**Author Contributions:** Conceptualization: M.A.J., S.R.A., W.K.A., B.A.M.; Analysis: M.A.J., W.K.A., B.A.M.; Writing—Original Draft Preparation: M.A.J., S.R.A.; Writing—Review & Editing: B.A.M., M.A.J., S.R.A., W.K.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to acknowledge the Punjab Irrigation Department for providing irrigation discharge and crop samples for LULC validation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*ISPRS International Journal of Geo-Information* Editorial Office E-mail: ijgi@mdpi.com www.mdpi.com/journal/ijgi

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18