A Novel Algorithm to Delineate Surface Water Paths on Digital Elevation Model Image with Boundary Element Method

It is well-known that the delineated surface water paths (SWP) from a DEM image are essential in hydrology. Both grid-based and contour-based algorithms were proposed in the literature to improve the delineated SWP position’s accuracy. Compared to the grid-based algorithm, the contour-based algorithms usually give more accurate results but require higher computation resources, especially in case of the wide catchment. This paper introduces a new contour-based algorithm which is formulated with the semi-analytical solution of Laplace’s partial differential equation with boundary element method. This approach allows the algorithm to determine the SWP in any direction in order to keep the delineated SWP smooth. The proposed algorithm was validated with the standard synthetic surfaces, where theoretical SWPs are known for accuracy evaluation. The obtained SWPs are more accurate than that of the popular grid-based algorithms. Moreover, the proposed algorithm requires less computation resources when considering very long contour lines. In experimentation with a real DEM image, all delineated SWPs are absolute (no broken part) when the contour interval is equal to or less than 20 meters, and the spacing between adjacent discrete elevation data is 20 meters or narrower. This algorithm helps the hydrologists estimate the catchment area, which is useful for water management in flood and drought prevention.


I. INTRODUCTION
D IGITAL elevation model (DEM) is a digital representation of a topographic map from remote sensing technology. In some parts of the world, DEM is synonymous with DTM (Digital Terrain Model). In the United States, DTM is defined as a model representing the bare earth surface, where all-natural and built features were removed [1]. Since current data capture technologies such as field survey, photogrammetry, remote sensing or Light Detection and Ranging (LiDAR) enable the high quality DTM, the DTM elevation data is very engaging for land-use planning, geographical modeling, hydrological modeling, and estima-tion of natural three-dimensional catchments [1], [2]. The size of the catchment area is an essential parameter in the hydrological model for estimating the water flow rate and evaluating the flood risk level. With the concept of a physicalbased model (mechanistic or white box model), the size of the catchment area can be physically observed from a catchment model that is drawn by the surface water path (SWP). The catchment model's effectiveness depends on the algorithm that is used to delineate the SWP from DTM. For the last few decades, there have been two main approaches to develop the SWP delineation algorithm, namely i) processing the DEM elevation data in the raster form [3], [4] and ii) processing the data in contour form [5]. To demonstrate that these algorithms are able to identify the water flow path, the SWP position error indicator is often used to measure the accuracy of these SWP delineation algorithms. It calculates the error of the SWP position relative to the theoretical SWPs or the real position of the SWPs. The relation between the resolution of DTM and the accuracy of SWP is not straightforward. For example, increasing the resolution in simple terrain case does not improve accuracy [4]. Therefore, a new algorithm is needed to delineate the SWP from the very large DEM image with adequate accuracy and affordable complexity. To improve the accuracy, an algorithm should exploit the water flow path model in addition to relying on raw data [6].
Numerous researchers give precedence to the macroscopic SWP characteristic, which is the effect of gravitation; therefore, the algorithms are usually formulated with only the streamline flow concept, not including the turbulent flow concept. A streamline is a curve whose tangent at any point gives the direction of the fluid velocity at that point and is time-invariant. This concept is analogous to a line of force in a static electric or magnetic field [7]. The papers [8]- [10] initiated this concept and demonstrated that their algorithms could generate the SWPs using the spatial distribution of slopes in the DTM. These algorithms became the fundamental tools in the hydrological models, and their extensions were grouped into three main approaches, namely the gridbased, TIN-based and contour-based algorithms [11].
In the approach to delineate the SWP with the grid-based algorithm, the D8 algorithm proposed by [12] was the pioneering work, followed by several papers that improve the algorithm's accuracy when the DEM image contains missing pixels. The extensions of its algorithm are divided into two main groups [13]: i) single flow direction (SFD) algorithms, e.g. Rho4 [14], Rho8 and D8-LTD [15] and ii) multiple flow direction (MFD) algorithms, e.g. FDFM [16], DEMON [17] and D∞ [18]. The enhanced D8 algorithms are employed widely in the hydrological model because the implementations are simple and require less computational resources. However, the errors of the grid-based algorithms is sensitive to the accuracy of the DTM [19]. The accuracy of a DTM depends on the pixel size, slope and quantization level [20]- [24]. They also found that the accuracy of the grid-based algorithms is limited with the oversimplified assumptions on the behavior of the water flow [10], [25], [26]. Later, [4] evaluated the accuracy of the six popular algorithms, namely D8, Rho8, D8-LTD, FDFM, MFD, and D∞. They compared the results simulated by these algorithms with the theoretical solutions of SWPs on natural watersheds. Unfortunately, they found that the error of the SWP positions falls in the range of 46.0% to 161.4%. The maximum is gained by D8 or FDFM algorithms, and the minimum by D8-LTD or D∞ approach.
In another approach, the algorithms delineate the SWP from contour lines, which are extracted from a DTM. The [9] initiated the contour-based algorithm with the streamline flow concept of [7], [27] to eliminate the spurious SWPs delineated by the grid-based algorithms in the case of divergent areas. With the algorithm's efficacy, many hydrologists, who need to estimate a sub-catchment area based on the SWP, implemented this algorithm in the catchment model [23], [28]- [32]. To further enhance this algorithm, a procedure that increases the algorithm's robustness while determining the SWP near the ridge is proposed in [29]. Later, the [30] adopt that procedure for partitioning a catchment into contour-based elements as simple and automated as possible. Furthermore, [31], [33]- [35] extend this approach to solve their problems in the hydrological model. However, [5], [32] pointed out that the algorithms mentioned above were still semi-automated and had some problems when delineating the SWP on the complex terrains such as ridges, saddles, and peaks. Finding a fully automated algorithm to deal with complex terrains is one of the most significant research challenges. Then, [5], [32] proposed the modified algorithm by using skeleton construction techniques, whereas [36] created another algorithm that is formulated with Triangulated irregular networks (TINs) structure. TINs were proved to be useful in many hydrological applications [37]- [41]. Unfortunately, the existing algorithms require high computational resources, especially for complex terrains, where the regions of high heterogeneity must be fitted with more data points. Although the existing contour-based algorithms give accurate SWP results, they are not attractive for a hydrological model because the contour-based algorithms still cannot be implemented with affordable computational resources [42].
Even until recently, the accuracy-complexity tradeoff of SWP delineation is still an open problem [4], [42]. Each application of SWP delineation requires a different level of accuracy. Existing algorithms can delineate SWP with sufficient accuracy for applications in geomatics, geographic information systems (GIS), mapping, and surveying [43]. On the other hand, the rainfall-runoff model application, which simulates the process to learn the water cycle for decision making and future prediction, requires higher accuracy, which is not achieved by existing algorithms [44]. One of the important applications is the sub-catchment area estimation, which is the essential input parameter of all hydrological models for studying hydrologic phenomena, hydrologic cycle, the impact of climate change, and soil properties on hydrology and water resources. This research aims to develop an algorithm with reasonable complexity to delineate the SWP with sufficient accuracy for sub-catchment area estimation. Obviously, there are two possible approaches: i) improving the accuracy of the grid-based algorithm and ii) reducing the complexity of the contour-based algorithm. Among the existing algorithms, both grid-based and contourbased algorithms delineate the SWPs by using elevation and slope data over the whole region, which is discretized into hydraulically connected elements such as grids or triangular elements. The grid-based algorithm used to delineate the SWP is more popular than the contour-based algorithm because of its simplicity and fewer computational resource requirements. However, it selects the most appropriate pixel as the endpoint of the SWP segment by observing only the VOLUME 4, 2016 surrounding eight pixels, which causes a local error. The local inaccuracy propagates to other pixels and accumulates into a serious error. With the limited resolution of the DEM image, the accuracy at each pixel cannot be improved, and hence the error issue of the grid-based algorithm cannot be mitigated. On the other hand, the contour-based algorithm gives high accuracy but is not affordable, and hence it is less popular in the implementation of the distributed hydrological model [42]. The high complexity part of the contour-based algorithm is the determination of the SWP direction at all grid points between two adjacent contour lines because the number of grid points considerably increases with the size of the whole region.
Recently, it was pointed out that the existing algorithms for delineating the SWP over DTM image do not meet the required accuracy of hydrological model development which is the important tool for understanding, predicting, and managing water resources [4], [42]. Increasing the DEM resolution in grid-based algorithms does not improve the position accuracy (error floor) since these algorithms exploit only neighboring pixels without extracting topographic information from larger surrounding area [4]. This motivated us to propose a new flow direction algorithm that exploits topographic information to significantly improve the position accuracy and reaches the error floor at higher DTM resolution. Also, the new algorithm must not incur an excessive computational time. The boundary element method (BEM), which provides an outstanding accuracy but demands high computational performance and large memory, becomes more practical with modern computers because it solves a 3-D problem with 2-D domain. For example, BEM was successfully used to develop a theoretical flow model for any potential function in [45]- [49]. Therefore, the practical BEM-based algorithm is proposed to delineate SWP over DEM image. The proposed algorithm is derived from numerical solutions of the governing partial differential equations.
The paper's central research objectives are as follows.
• The elevation function in the DEM image was known as the contour lines obtained from image processing technique [50]. The elevation function between adjacent contour lines is unknown and can be described by 2-dimensional Laplace's partial differential equation, where the known contour lines serve as the boundary values. Our first objective is to formulate the semianalytical solution of Laplace's partial differential equation with BEM, which discretizes the problem domain along the boundary into piecewise linear curves, where the domain can be parameterized by one variable. The semi-analytical solution includes the topographic information, which improves the delineated SWP accuracy. Moreover, comparing to the TIN structure, the problem domain is reduced from an area between contour lines to contour lines. With the significant accuracy improvement and practical computational complexity, the semianalytical solution makes a major contribution to the urgent needs of hydrologists to cross the current limit of hydrological modeling. • The next step is to delineate every SWP, which is a path starting from a source and ending at a sink, one by one. All sources are the points on contour lines, while all sinks are the points between contour lines. For each SWP, we first determine the SWP direction at a source and then iteratively determine the SWP direction until reaching a sink. In doing that, we propose an algorithm that uses the previously obtained elevation function as well as its first and second derivatives in the polar coordinates to determine the SWP direction at every iteration from the source to the sink. Note that we can also delineate the SWP in the reverse direction from the sink to the source because our algorithm basically finds the steepest downslope or steepest upslope direction. Compared to other contour-based algorithms, which delineate an SWP by connecting two points on a pair of adjacent contour lines, our algorithm gives greater continuity along the SWP. Compared to gridbased algorithms, which limit the number of directions to eight at every pixel, our algorithm gives a more exact direction at every point along the SWP. The continuity and smoothness of the proposed algorithm mitigate the issue of grid-based algorithms that cease delineating SWP when they reach a flat plane surface. • In practice, the proposed algorithm is applied with real terrain and must be able to cope with two causes of complexity rising. First, the contour line encircling a foothill can be very long. Second, DTM with high resolution should be sampled at high sampling rate. Both of them creates a great number of sampling points, which incur a time-consuming computation. The third objective is to design the algorithm that can delineate SWPs with incomplete contour lines, which have missing contour segments and much less sampled points as shown in Fig. 1. Therefore, the proposed algorithm requires less computational resources than other contour-based algorithms.
The remainder of this article is structured as follows. Section II proposes the numerical solution of 2-dimensional Laplace's partial differential equation formulated with BEM. This solution is the elevation function that gives the elevation at any coordinate and is defined by contour integrals, where the elevations along the contours are known. Next, in Section III, we present how to use the elevation data in a DEM image to identify the gradient of the elevation along the contour lines which are the unknown boundary value functions. Both the elevation and the obtained gradient of the elevation along the contour lines are used to interpolate the elevation at any point between the adjacent contour lines in Section IV. Next, Section V presents how to use the contour integrals and the gradient terms for determining the SWP direction at any point along the delineated SWP, and describes the proposed algorithm for delineating the SWP from DEM elevation. In Section VI, we demonstrate and discuss the efficiency of the proposed algorithm on delineating SWPs over both synthetic surfaces and the real terrain. The synthetic surfaces are used to evaluate the position error of the delineated SWPs. The real terrain with very long contour lines is used to validate the proposed algorithm. Last, we conclude the paper in Section VII.

II. PROPOSED FORMULATION OF LAPLACE'S EQUATION WITH BEM
This section proposes the formulation of elevation function, which is the essential function for analyzing the directions of SWPs at any coordinates between the DEM contour lines. To improve the elevation interpolation in the contourbased algorithm of SWP delineation, numerous researchers pay attention to the elevation function formulation with numerical techniques. Reference [42]  solved. If the unknown variables are elevation gradients on a segment of DEM contour lines, and that segment is not close to SWPs, then they are not needed to be solved.

A. ANALYTICAL SOLUTION OF LAPLACE'S PDE
A flow line is a path starting at a point on a contour line and ending at a point on the adjacent contour line. With the streamline flow concept, the flow line directions at starting and ending points are perpendicular to their contour lines. Hence, the flow line directions at starting and ending points can be found by calculating the elevation gradients at these points. However, the flow line directions at other points on the path are unknown because the elevation data outside contour lines are not available. There are several algorithms that interpolate these elevation data [5], [9], [27], [29]- [32], [34]. Different types of interpolation give different elevation functions and different patterns of SWP. The most accurate algorithm finds the analytical solution of Laplace's partial differential equation in the context of the conservative vector fields. The solution is used to interpolate terrain surface from the elevation data along the contour lines [51]. Laplace's partial differential equation is considered, because the boundary value functions, including elevation data and their spatial derivatives, are time-invariant. The solution is solved by the boundary integral equation method. Because of its accuracy, we apply this analytical solution of Laplace's partial differential equation for interpolation. The interpolated elevation data and spatial derivatives are then used to determine the directions at all points along the flow line as shown in Section V.
Many kind of techniques can be used to solve the Laplace's equation. In [51], they selected the boundary integral equation method with Green's theorem [52] to solve the analytical solution of Laplace's equation. The obtained particular solution can be written as The boundary C contains any point where the elevation is known and that point is on any contour line. The boundary C may contain points on several contour lines. For example, the boundary C contains the points on both the lower contour line H L ( ⃗ Q) and the upper contour line H U ( ⃗ Q). The boundary C is not only the domain of the known elevation function H( ⃗ Q) but also the domain of the potential gradient ∇H( ⃗ Q). The direction of the boundary C at the point Q is specified with a unit normal vectorn( ⃗ Q). The value of ψ is determined by a point ⃗ P 's neighborhood, which is the part of a small disc around a ⃗ P point inside the boundary C. For example, considering a point ⃗ P on a contour line, the neighborhood of a point ⃗ P is a half-disc, and so ψ = π. If the point ⃗ P is between the contour lines, then the neighborhood of a point ⃗ P is a full disc, and so ψ = 2π.
The elevation function h( ⃗ P ) in (2) can be used to model the SWP accurately with an assumption that no contour line is missing, when the domain Ω is bounded by the contour lines C j , where j is the index of a contour line. For example, considering two adjacent contour lines with different elevations, the obtained SWP will be a path from the upper contour line to the lower contour line. As another example, if a contour line obstructs between two adjacent contour lines, the path of SWP will make a detour around the obstructing contour line. Moreover, this analytical approach can generate the continuous SWP over the domain Ω and the characteristics of the obtained SWP conform to the streamlined concept. The obtained SWP represents the water flow from the upper to the lower contour line. The flow is guaranteed not to be in a reverse direction because the value of the elevation function h( ⃗ P ) is bounded between the local minima (the elevation value of the contour line H L ( ⃗ Q)) and the local maxima (the elevation value of the contour line H U ( ⃗ Q)).

B. PROPOSED SOLUTION WITH BEM FOR DEM CONTOUR LINES
The analytical solution mentioned in (2) can reproduce the elevation function h( ⃗ P ) from the known elevation on the boundary H( ⃗ Q). However, it cannot promptly interpolate the elevation data from DEM contour lines. Here, we aim to take advantage of both the analytical solution and the DEM image.  There are two ways to do interpolation between the contour lines for SWP delineation.
First, we must represent the DEM elevation data along the contour lines into the elevation function in analytic form or closed-form. Then, we calculate its gradient in the analytic form. The available elevation data may be employed to approximate the function H( ⃗ Q), but it is very complicated to approximate the gradient function along the fixed contour lines. In practice, the contour-based DEM elevation data cannot be automatically described with an analytic function on the computer because the elevation is random and vast. Hence, this paper proposes the second way, where we approximate the DEM elevation data along a contour line with many piecewise Lipshitz functions [53]. Then, we apply these piecewise functions in the BEM process to obtain the numerical solution by reformulating the analytical solution in (2). In the mentioned analytical solution, a whole boundary integral is replaced with the definite boundary integrals on the piecewise intervals. The number of intervals equals the number of piecewise elevation functions. The numerical solution is formulated as follows.
Consider a boundary C. The boundary integral of (2) is split into N definite integrals, each of which is indexed by j = 1 . . . N . Each interval is denoted by C j . Then, we obtain the semi-analytical solution: For example, we implement (3) for two adjacent contour lines. In Fig. 3, the two contour lines from a DEM image are the boundary which is surrounding the interior region Ω. When we select N sample points along two contour lines, the lower and upper contour lines contain N C L and N C U intervals, respectively. Then, we obtain The boundary values defined in the semi-analytical solution comprise the elevation function H( ⃗ Q) and the normal component of the potential gradient ∇H( ⃗ Q). The normal and denoted by ∇H n . Both functions on a point in the interval C j are approximated with the polynomial functions of a dummy variable s: where ζ j , ξ j , α j , β j , γ j and η j are constants, and s is the distance from the middle of segment C j , and L j is the length of segment C j and is written as The accuracy of this approximation depends on two conditions. First, the elevation function is smooth within each interval, so that the potential gradient is continuous. Second, the number of definite split integrals N is large enough that the difference between the original and the approximated contour line is small. Note that the larger number of split definite integrals incurs a longer computational time. In fact, increasing the polynomial orders of H j (s) and ∇H n,j (s) also improves the approximation accuracy, but the computational time sharply increases. Note that this approximation can be applied with both complete and incomplete contour lines as shown in Fig. 3 (a) and (b), respectively. This approximation gives the polynomial coefficients that will be used to determine the boundary value function in the numerical solutions in Section III-A.
With polynomial approximation, we substitute H j (s) and ∇H n,j (s), defined in (5) and (6), respectively, for H( ⃗ Q) and ∇H n ( ⃗ Q) in (3). The semi-analytical solution in (3) becomes where The F αj , F βj , F γj , F ηj , F ζj and F ξj are geometric terms, referred to as shape functions. They are the functions of the boundary (contour line's path) and the position of the point ⃗ P . Both boundary and the point ⃗ P will be found in Section III.

III. ANALYSIS OF BOUNDARY VALUE FUNCTION
This section presents how to use the DEM contour line for identifying the boundary value functions in (8). The DEM contour line is a discrete sequence of referenced points located along the contour line, and is available from both grid-based DEM and contour-based DEM [54]. The discrete sequence of points is used to calculate the elevation function H( ⃗ Q) along the boundary and the six shape functions F * j ( ⃗ P , ⃗ Q) in Section III-A according to our numerical solution formulation in (8). Then, the obtained H( ⃗ Q) is used to numerically compute the gradient of the elevation function along the contour line as shown in Section III-B. The gradient will be used to interpolate the elevation function between the contour lines and to delineate the streamlines in Section V.

A. CALCULATION OF ELEVATION FUNCTION ON THE CONTOUR LINES AND THE SHAPE FUNCTIONS
To implement the proposed numerical solution, we use the information of the DEM contour lines which comprise the coordinates of the 2D points along the DEM contour lines VOLUME 4, 2016 and the elevations at those points. We use the position information of the points along the boundary C j for calculating the integration of the shape functions as shown in (9) to (14), and use the elevation levels at those points for approximating the elevation function H j ( ⃗ Q) along subdomain C j as shown in (5).
In Fig. 3 (a), the DEM contour lines were systematically sampled into points, whose coordinates define the points ⃗ Q along the subdomain C j . The DEM was also sampled into points to define the points ⃗ P in the domain Ω. The coordinates of both ⃗ Q and ⃗ P are used to calculate the integration of the shape functions, namely F αj , F βj , F γj , F ηj , F ζj and F ξj . The subdomain C j comprises the points along the straight line from point ⃗ Q j to point ⃗ Q j+1 . Its direction is from ⃗ Q j to ⃗ Q j+1 , and is used to identify the unit normal vectorn( ⃗ Q). Before integrating the shape functions, we transform the coordinates of ⃗ P and ⃗ Q into the local coordinates which are defined with the variable s so that the line integration can be done in closed-form.
After defining the coordinates of the point ⃗ Q on each subdomain C j , we use the elevation data at those points for approximating the elevation function H j ( ⃗ Q). The elevation levels at ⃗ Q j and ⃗ Q j+1 are used to calculate the coefficients of the approximated functions in (5). Usually, DEM image contains both ordinary contour lines and incomplete contour lines. In case of ordinary contour lines, each contour has the constant elevation data denoted H 0 , so the coefficients η and ζ are zero and H 0 , respectively. In case of very long or incomplete contour lines, the elevation levels over the gap of the contour lines are not constant, and we approximate the varying elevation levels on subdomain C j by piecewise linear function as shown in Fig. 3 (b). For example, the elevation levels at the point ⃗ Q j and the point ⃗ Q j+1 are H j and H j+1 , respectively. The coefficients η j and ζ j are the slope and H j ( ⃗ Q)-intercept, respectively.

B. ANALYSIS OF ELEVATION GRADIENT ON THE CONTOUR LINES
To present the process to determine four unknown coefficients for any subdomain of the boundary integral, we divide the whole boundary C into the N small boundary C j , where j = 1, . . . , N . The process is structured as follows.

1) Observing along subdomain Ci
In Section II-A, the function reproduced by the analytical solution is the elevation function h( ⃗ P ), where the point ⃗ P is in either domain Ω or boundary C. In Section II-B, we put the observed point ⃗ P in the domain Ω, and we employ the reproduced function from the numerical solution for interpolation at the points between adjacent contour lines. Here, we put the observed point ⃗ P along the boundary C and apply the numerical solution for determining the gradient terms. Now, we have the 4N unknown coefficients of the function ∇H n ( ⃗ Q); therefore, we will put the observed point ⃗ P at 4N different positions along the boundary C.
To identify the gradient terms in the numerical solution, we resolve the four unknown coefficients, which comprise α j , β j , γ j , and η j , of the function ∇H n ( ⃗ Q) (defined in (6)) in each boundary C j by solving a set of linear equations. We formulate a set of four equations with all unknown coefficients by selecting four different observed points ⃗ P along a subdomain C i , where C i denotes the i-th subdomain of the boundary C and is the location of four selected points. Along the boundary C i , we put the observed points ⃗ P at For example, when putting ⃗ P at s = −1/2 on C i , we obtain: where h( ⃗ P − 1 2 ) can be determined with the elevation function H i (s = − 1 2 ) in (5): The shape functions of ∇H n ( ⃗ Q) can be determined with the integrations: and the shape functions of H( ⃗ Q) can be determined with the integrations: Similarly, when putting ⃗ P at s = − 1 4 , 0 and 1 4 on C i , h( ⃗ P − 1 4 ), h( ⃗ P (0)) and h( ⃗ P 1 4 ) can be determined with the elevation function H i s = − 1 4 , H i (s = 0) and H i s = 1 4 , respectively.

2) Rewriting with vector equation
Now, we obtain a set of four equations after setting ⃗ P at four different positions on the subdomain C i . Obviously, they are vector equations and can be rewritten in an alternative form: where i = 1, · · · , N is the index of each set of four ⃗ P positions on the same subdomain C i , and j = 1, · · · , N is the index of subdomain C j , which serves as the j-th subdomain of the boundary integration. [H H H i ] is 2 × 1 boundary vector, which contain two coefficients of a function H( ⃗ P ) along C i and is defined as [∇H ∇H ∇H n,j ] is a 4N × 1 boundary vector, which contains four coefficients of a function ∇H n ( ⃗ Q) along C j and is defined as (22) [ψ ψ ψ i ] is a 4 × 2N block diagonal matrix, which contains four neighborhood factors of each set of four ⃗ P positions along C i and is defined as is 4 × 2N geometry matrix, which contains two integrals of shape functions F ζj ( ⃗ P i ) and F ξj ( ⃗ P i ) at ⃗ P on C i and at ⃗ Q on C j and is defined as

3) Solving equation system
To create the equation system for determining the [∇H ∇H ∇H n,j ] vector, we repeat the process of equation formulation with (19) for all subdomains C i , with i = 1 . . . N . Then, we obtain N sets of four vector equations, which can be written as a linear equation system: Next, we obtain the gradient vector [∇H ∇H ∇H n,j ] by solving the equation system: . . .

IV. CALCULATION OF ELEVATION AND ITS GRADIENTS BETWEEN THE CONTOUR LINES
The existing contour-based algorithm delineates SWP by employing the interpolation techniques that must discretize the problem domain, which is the interior area between two adjacent contour lines. The unknown elevations over the whole domain must be solved to do the interpolation. This incurs high complexity in a process to delineate the SWP. To reduce the computational resource requirements, we propose the process to delineate the SWP by interpolating only some points between the contour lines. The proposed process focuses only the points which are used to delineate the SWP, and interpolate the elevation data at those points and their neighborhoods. We solve the elevation function, representing the change of elevation levels at those points and their neighborhoods by the proposed numerical solution in (8) in Section II-B. The obtained elevation function will be used to determine the direction of an SWP in Section V. The elevation function comprises the elevation function h( ⃗ P ) and its derivatives at the point ⃗ P . Therefore, this section presents how to interpolate the elevation at the points between the contour lines and how to calculate the first and the second derivatives of the elevation function.

A. INTERPOLATED ELEVATION FUNCTION OVER REGION BETWEEN DEM CONTOUR LINES
Delineating SWP between DEM contour lines is to recursively calculate the next ⃗ P , starting from a chosen point between contour lines. The recursion will lead ⃗ P to the contour line with lower elevation. The recursive algorithm involves with the elevation function and its derivative. The elevation function h( ⃗ P ) is formulated with the proposed numerical solution in (8) in Section II-B. Firstly, we rewrite (8) in the vector form: where G G G H,j ( ⃗ P ) is 1 × 2N geometry vector which contains two integrals of shape functions: F ζj ( ⃗ P ) and F ξj ( ⃗ P ). The  (13) and (14) in Section II-B. Therefore, G G G H,j ( ⃗ P ) is defined as [H H H j ] is 2N × 1 boundary vector which contains two coefficients of a function H( ⃗ Q) along C j . These coefficients are equal to the components of boundary vector [H H H j ] in (26) and are used in solving the boundary value function in Section  III-B. [G G G ∇H,j ( ⃗ P )] is 1 × 4N geometry vector which contains four integrals of shape functions: F αj ( ⃗ P ), F βj ( ⃗ P ), F γj ( ⃗ P ) and F ηj ( ⃗ P ). The integration results are in (9), (10), (11), and (12) in Section III-A. Therefore, [G G G ∇H,j ( ⃗ P )] is defined as [∇H ∇H ∇H n,j ] is 4N × 1 boundary vector which contains four coefficients of a function ∇H ∇H ∇H( ⃗ Q) along C j . These coefficients are solved with the equation system in (26) in Section III-B3.

B. GRADIENT OF INTERPOLATED ELEVATION FUNCTION OVER REGION BETWEEN DEM CONTOUR LINES
There are several ways to find the derivative of elevation function. The slope at a point is the essential information to delineate SWP. Most methods calculate the slope at a designated point from the interpolated elevations of two points. The interpolation in existing methods is done in several directions, and only the one that gives the steepest slope is chosen. In this section, we present a different approach, where the slope can be immediately calculated from a closed-form formula. Specifically, the closed-form formula is derived by taking the derivation of the numerical solution in (3) at the point ⃗ P . After obtaining the boundary value functions in Section III, namely the elevation function and its gradient on the contour lines, we use them to calculate the interpolated elevation function and the slope at same point ⃗ P , mentioned in Section IV-A. The slope ∇ P h( ⃗ P ) at point ⃗ P comprises the derivative of the elevation function h( ⃗ P ) with respect to x (h x ) and the derivative of the elevation function h( ⃗ P ) with respect to y (h y ). To obtain two derivatives at point ⃗ P , we determine the gradient of a function h( ⃗ P ) by applying the vector operator ∇ to the scalar function h( ⃗ P ). The scalar function h( ⃗ P ) is reproduced from the proposed semi-analytical solution in (3). The process to calculate the gradient of elevation function h( ⃗ P ) at the point between the DEM contour lines is structured as follows.

1) Formulating the gradient of proposed solution
Again, we apply (3) in Section II-B to formulate the gradient of semi-analytical solution at point ⃗ P . Taking the gradient operator to (3) gives where ∇ P Θ and ∇ P Ξ are gradient terms in brackets, and can be calculated analytically. Then, we obtain the resulting expression as

2) Approximating the boundary value function
To formulate the numerical solution of ∇ P (h( ⃗ P )), we substitute the polynomial approximations of H( ⃗ Q) and ∇H n ( ⃗ Q) functions (which are defined in (5) and (6), respectively) into (29). It is similar to obtaining the numerical solution of h( ⃗ P ) in Section II-B. Then, the semi-analytical solution in (29) becomes the numerical solution in (47) in Appendix A.

3) Rewriting with vector form
After obtaining the numerical solution of ∇ P (h( ⃗ P )) in (47), we can use it to determine the value of the gradient function ∇ P (h( ⃗ P )) with the the boundary value functions, namely H( ⃗ Q) and ∇H( ⃗ Q) ·n( ⃗ Q). We substitute the coefficients of the boundary value functions from the vector [H H H j ] and [∇H ∇H ∇H n,j ] in (26) into (47) and rewrite it in a vector form as where 2N geometry vector, which contains two integrals of shape functions: F ′ ζj ( ⃗ P ) and F ′ ξj ( ⃗ P ). Both shape functions are of the point ⃗ P locating between contour lines and of the point ⃗ Q locating on C j , and are defined as ] is 1 × 4N geometry vector, which contains four integrals of shape functions: F ′ αj ( ⃗ P ), F ′ βj ( ⃗ P ), F ′ γj ( ⃗ P ) and F ′ ηj ( ⃗ P ). These four shape functions are of the point ⃗ P locating between the contour lines and of the point ⃗ Q locating on C j , and are defined as

C. GRADIENT OF GRADIENT OF INTERPOLATED ELEVATION FUNCTION OVER REGION BETWEEN DEM CONTOUR LINES
In this paper, a second-order series expansion will be used to determine the direction of the SWP over any point in the catchment basin with the proposed algorithm in Section V.
The first derivative at any point ⃗ P was found in Section IV-B. This section finds the gradient of the gradient of interpolated elevation function, which is a tensor and can be written as where h xx = ∂h x /∂x, h xy = ∂h x /∂y, h yx = ∂h y /∂x and h yy = ∂h y /∂y. The vector function ∇ P h( ⃗ P ) is reproduced from the proposed semi-analytical solution as shown in (29) in Section IV-B1. The process to calculate the gradient of gradient of the elevation function h( ⃗ P ) at the point between the DEM contour lines is structured as follows.

1) Formulating gradient of gradient of proposed solution
From (29) in Section IV-B1, we formulate the gradient of ∇ P (h( ⃗ P )) as where ∇ P Λ and ∇ P Ψ are gradient terms, and can be calculated analytically as follow.
and I I I is an identity tensor, which can be defined as and ⃗ P − ⃗ Q 2 ,n ⃗ P − ⃗ Q and ⃗ P − ⃗ Q n are also tensors, where 2) Approximating the boundary value function To formulate the numerical solution of ∇ 2 P (h( ⃗ P )), we substitute the polynomial approximations of H( ⃗ Q) and ∇H n ( ⃗ Q) functions (which are defined in (5) and (6)) into (36). Similar to obtaining the numerical solution of h( ⃗ P ) in Section II-B, the semi-analytical solution in (36) becomes the numerical solution in (54) in Appendix B.

3) Rewriting with vector form
After obtaining the numerical solution of ∇ 2 P (h( ⃗ P )) in (54), we can use it to determine the value of the gradient of gradient function ∇ 2 P (h( ⃗ P )) with the the boundary value functions, namely H( ⃗ Q) and ∇H( ⃗ Q) ·n( ⃗ Q). We substitute the coefficients of the boundary value functions from the vector [H H H j ] and [∇H ∇H ∇H n,j ] in (26) into (54) and rewrite it in a vector form as where G ′′ G ′′ G ′′ H,j ( ⃗ P ) is 1 × 2N geometry vector, which contains two integrals of shape functions: F ′′ ζj ( ⃗ P ) and F ′′ ξj ( ⃗ P ). Both shape functions are of the point ⃗ P locating between contour lines and of the point ⃗ Q locating on C j , respectively. They are defined as and [G ′′ G ′′ G ′′ ∇H,j ( ⃗ P )] is 1 × 4N geometry vector, which contains four integrals of shape functions: F ′′ αj ( ⃗ P ), F ′′ βj ( ⃗ P ), F ′′ γj ( ⃗ P ) and F ′′ ηj ( ⃗ P ). These shape functions are of the point ⃗ P locating between contour lines and the point ⃗ Q locating on the C j . They are defined as

V. DELINEATION OF SURFACE WATER PATH
This section proposes an algorithm to delineate the SWP over any point in the catchment basin. Unlike other algorithms, the proposed algorithm can estimate the size of a sub-catchment area, which is only a part of the whole area, because it can generate the SWP without interpolation over all points VOLUME 4, 2016 between the contour lines. This feature helps us reduce computational time because generating the SWP requires only the information in the area that the SWP is passing. The proposed algorithm comprises three main steps. The first step is to prepare the information to predict an endpoint of the SWP segment after identifying the starting point on the subcatchment area. The starting point information comprises the elevation and its gradients at that point. It can be calculated with the numerical solution in Section IV. In the second step, we predict the endpoint of the SWP segment based on the approximated surface in polar coordinates. We solve a second-order series expansion to approximate the elevation of a neighborhood of the starting point by using the numerical results obtained from the first step. The third step is to predict the next SWP segment. We use the endpoint of the SWP segment as the starting point of the next SWP segment. Therefore, this section presents the second and third steps of the algorithm.

A. DETERMINATION OF THE SWP SEGMENT
Before delineating the whole SWP from an upper contour line to a lower contour line, we must determine the location of each SWP segment. We choose a point between contour lines as the starting point of the first segment and choose the endpoint from a neighborhood of the starting point where the endpoint has the maximum or minimum elevation. In doing that, we approximate the elevation over the neighborhood of the starting point by using the surface in polar coordinates.

1) Creating surface over the neighborhood of the starting point
In Section II, we use the solution of Laplace's equation to interpolate the elevation at any point along the SWP. Here, we still solve its solution in a polar coordinate system, and apply that solution for creating surface over the neighborhood of the starting point of the SWP segment. Laplace's equation in a two-dimensional polar coordinate system is written as where r is the distance from the starting point of the SWP segment to the endpoint, and θ is between 0 and 360 degrees.
We can obtain its solution in general form by using the separation of variables. Then, the solution in a particular form is obtained by fitting this application and can be written as h(r, θ) = a 0 + ∞ m=1 r m (a m cos mθ + b m sin mθ). (45) Even though the series expansion in (45) is exact, it is an infinite series. For computation efficiency, we substitute a m and b m for the finite number of terms. In doing this, we keep the expansion area small (r < r 0 ) to limit the approximation error, where r 0 denotes the length of the SWP segment. Hence, SWP is a concatenation of several short SWP segments, each of which has its own series expansion. At the beginning, the surface over the neighborhood of the  starting point of SWP segment in (45) is approximated as (71) in Appendix C. The direction over this surface will be chosen to delineate the first segment of SWP in the next step.

2) Determining SWP direction
The SWP direction is emanating from the starting point of the SWP segment and toward the endpoint of that segment or the starting point of the next SWP segment. Given that the small disc around the starting point is contained in the region between the contour lines, we select the endpoint of the current SWP segment from the neighborhood of the starting point. Otherwise, no endpoint is selected, and SWP delineation finishes. When setting the length of the SWP segment as r 0 , we consider the elevation at any point on the boundary of the disc from the approximated surface in the polar coordinates h(r, θ), which can be determined in (71). Fig. 4 shows an example of approximated surface in the polar coordinates h(r, θ) for r < r 0 . We select the path, which starts from the point ⃗ P toward the endpoint that gives steepest descent. To do that, we find the endpoint, which minimizes h(r, θ) on the circle with r = r 0 . In fact, another approach is to select the path, which starts from the point ⃗ P toward the endpoint with the steepest ascent. To do that, we find the endpoint, which maximizes h(r, θ) on the circle with r = r 0 . In this paper, we use the steepest descent approach. Therefore, we identify the endpoint that minimizes h(r, θ) by finding θ min , where dh(r, θ)/dθ = 0 on the circle r = r 0 , as shown in Fig. 4 (b). An example of θ min finding with standard technique is given in Appendix D.

B. PROPOSED ALGORITHM
In this section, we propose a recursive algorithm for delineating an SWP, which is the concatenation of several SWP segments over a catchment area, in which the elevation is represented with contour lines in a DEM image. The algorithm comprises five phases as follows.

1) Allocating the memory on a computer
To delineate the SWP over a large catchment area, the algorithm should efficiently use computer memory to collect the contour line data and the position of the delineated SWPs. The memory management must allocate portions of memory and free them for reuse when no longer needed. The input data comprises the coordinates of the contour lines and a list of zoning. To apply the contour lines for being the boundary condition in the proposed algorithm, we regularly sample some points along the original contour lines obtained from a DTM and reconstruct the polylines (chain of line segments) as the input data. The number of segments depends on the sampling rate. A list of zoning is information describing the boundary of each sub-area. Fig. 5 presents the elevation over a catchment area, which is represented with N contour lines and K zones. The zone z 1 is bounded by the C 1 and C 2 , and the zone z K is bounded by the C N −1 and C N . The algorithm collects each coordinate, denoted as ⃗ Q Cn j , of the contour line C n . The output data comprises the boundary values along the contour lines C n , namely the elevation functions and the gradient of elevation function, denoted by H( ⃗ Q Cn j ) and ∇H( ⃗ Q Cn j ), respectively. These output data will be used to interpolate the elevation at the points, which the SWP is passing. This point is determined from the starting point of the m th segment as explained in Section V-A. The starting point and the endpoint of the m th segment are denoted by ⃗ P m and ⃗ P m+1 , respectively. The starting point of the 1 st segment ⃗ P 1 is arbitrarily designated from any point in the region lying within two adjacent contour lines (or any zone; e.g., z 1 , z 2 , . . . , or z K ). Notice that the 1 st segment starts from the upper contour line, and the last segment ends at the lower contour line.

3) Preparing information at ⃗ Pm
Before the algorithm computes the endpoint of the current SWP segment, it requires the elevation and its derivatives at point ⃗ P m . In this phase, hence, the algorithm computes h( ⃗ P m ) and its derivatives (∇h( ⃗ P m ), ∇ 2 h( ⃗ P m )) with the numerical solutions which are explained in Section IV. These computations require a part of the input data collected in the memory. For example, in Fig. 5, the point ⃗ P m is in the zone z 1 , so the algorithm requires only the input data involving the contour lines C 1 and C 2 . The required input data comprise the coordinates of points, i.e., the ⃗ Q C1 j and ⃗ Q C2 j , as well as the elevation and gradient functions along the contour lines, i.e., The calculations of both H( ⃗ Q Cn j ) and ∇H( ⃗ Q Cn j ) are explained in Section III.

4) Determining direction of the SWP segment
After the elevation and its derivatives at point ⃗ P m , namely h( ⃗ P m ), ∇h( ⃗ P m ) and ∇ 2 h( ⃗ P m ) were obtained, this phase employs those calculated results for determining θ min . The calculation of θ min is explained in Section V-A. Then, the algorithm uses the θ min to determine the endpoint of current SWP segment. Note that the length of the SWP segment (r 0 ) can be adjusted, and this will give a different endpoint.

5) Determining the starting point of the next SWP segment
In the last phase, the algorithm collects the coordinates of the endpoint of the current SWP segment into the memory and starts to delineate the next SWP segment by using the endpoint of the current SWP segment as the starting point ⃗ P m+1 . The algorithm repeats from the second phase until the endpoint of the SWP segment is not located in the zone z N .
All five phases can be written as pseudocodes as shown in Algorithm 1 and Algorithm 2, respectively.
To summarize, the computation is initiated by determining the endpoint of the first SWP segment from the assigned point with Algorithm 1. Then, Algorithm 1 is looped until the endpoint of the SWP segment locates outside the catchment area. Finally, each loop calls Function Next_SWP in Algorithm 2 to determine the endpoint of the current SWP segment, which is the starting point of the next SWP segment.

VI. RESULTS AND DISCUSSIONS
This section demonstrates how to delineate the SWP over the synthetic and natural surfaces with the proposed algorithm, and presents the validation of the extracted SWP. First, we study the position accuracy of the SWP extracted by using the elevation data derived from the formulas of standard synthetic surfaces, namely ellipsoid surface (hill), inverted ellipsoid 1 (pit), and inclined plane (flat land). Then, we use the elevation of natural terrain based on a DEM image which comprises the hill's features, pit, and flat land, to study the proposed algorithm's performance and efficiency. 1 Inverted ellipsoid is referred to as inverse ellipsoid in [4].
This section demonstrates how to delineate the SWP over the synthetic and natural surfaces with the proposed algorithm, and presents the validation of the extracted SWP. First, we present the performance of the proposed algorithm when there is no incomplete contour line, and analyze the position accuracy of the delineated SWP by using the elevation data derived from the formulas of standard synthetic surfaces, namely ellipsoid surface (hill), inverted ellipsoid (pit), and inclined plane (flat land). Then, the case that there are incomplete contour lines is considered. We use the elevation of natural terrain based on a DEM image, which comprises the hill's features, pit, and flat land. The SWPs delineated from the incomplete contour lines are compared with the SWPs delineated from complete contour lines to evaluate the effectiveness of the proposed algorithm.

A. SYNTHETIC SURFACE
Reference [4] presented the method for measuring the position error of extracted SWPs to systematically evaluate the accuracy of various flow direction algorithms, namely the D8 algorithm and its extensions. The position error is the difference between the extracted SWP position and the theoretical (true) SWP position derived from the formulas of the standard synthetic surface. Furthermore, the results in [4] reveal that the D8 algorithm and its extensions have nonnegligible error floors. Therefore, it is of interest to use this evaluation method to also measure the position error of our proposed algorithm, so that the comparisons can be made.

1) Theoretical SWPs
We use the theoretical SWP positions as the references for measuring the error of the extracted SWPs. The theoretical SWP is the closed-form solution which is derived from the formulas of the synthetic surface mentioned in [4]. Table 1 lists the formulas of three synthetic surfaces and the equations of theoretical SWPs. The 3-D graphics of the synthetic surfaces are illustrated in Fig. 6.
The positions of theoretical SWPs in three cases, namely ellipsoid, inverted ellipsoid, and inclined plane, are plotted with red lines in Fig. 7, Fig. 8 and Fig. 9, respectively. First, the ellipsoid that models a hill in Fig. 6 (a) gives the theoretical SWP, radiating out from the top of the ellipsoid as shown in Fig. 7. Second, the inverted ellipsoid that models a pit in Fig. 6 (b) gives the theoretical SWP, sinking from the border of the inverted ellipsoid to its center as shown in     Fig. 9.

2) Extracted SWPs over ellipsoid
This work uses the formulas of the ellipsoid listed in Table 1 to generate the contour lines which are in a circle shape. The ellipsoid comprises five contour lines, namely 2000,1998,1964,1854, and 1653 meters. Along a contour line, we select a point at every θ s degrees. For a fair comparison, θ s is set at 10 degrees, so that the distance between adjacent points is not shorter than that of D8 resolution (4.5 × 4.5 meters). ordinates of ⃗ Q points. We set the length of an SWP segment as r 0 = 1.51m., and locate the starting point ⃗ P 1 in the zone z 1 , which is bounded by the contour lines with h = 2000m and h = 1998m. To apply the method of [4], we put the starting points very close to the position, used to evaluate the error of the D8 algorithm. The distance between the starting points and the ellipsoid center is set at 30 meters, so that the coordinate of the starting points is close to the corner point of a DEM pixel. We select three starting points along a circle with a radius of 900 meters to delineate three extracted SWPs in the different directions, namely 135, 60, and 270 degrees, where the extracted SWPs are labeled by "Line-1", "Line-2", and "Line-3", respectively and are plotted with blue markers in Fig. 7.
To evaluate the accuracy of three extracted SWPs, we determine the position error of the extracted SWP by using the formula where θ Exact is the constant angle of the theoretical SWP, and θ m Cal is the angle of the extracted SWP which is converted from the coordinates of the ending point of the m th segment of extracted SWP. First, the positions of the "Line-1" extracted SWP in Fig. 7 are compared with θ Exact = 135 degrees, and the position error E(%) versus the distance of the SWP flowing from the ellipsoid center are plotted in Fig.  10 (a). Obviously, the value of E(%) is below 0.001 at all distances. Second, the positions of the "Line-2" extracted SWP are compared with θ Exact = 60 degrees, and the position error E(%) versus the distance of the SWP flowing from the ellipsoid center are plotted in Fig. 10 (b). The value of E(%) is below 0.1 at all distances. In case that θ s is increased to 20 degrees, the value of E(%) is below 0.2 at all distances. Last, the positions of the "Line-3" extracted SWP are compared with θ Exact = 270 degrees, and the position error E(%) versus the distance of the SWP flowing from the ellipsoid center are plotted in Fig. 10 (c). The value of E(%) is below 0.0002 at all distances.

3) Extracted SWPs over inverted ellipsoid
This work uses the formulas of the inverted ellipsoid listed in Table 1 to generate the contour lines, which are in a circle shape. The inverted ellipsoid comprises five contour lines, namely -2000, -1998, -1964, -1854 and -1653 meters, respectively. Notice that the coordinates of contour lines are the same as those in ellipsoid, but their elevations are negative. They are plotted in Fig. 8. Totally, a boundary C in the numerical solution in (3) consists of 5 × (360/θ s ) intervals (C j ).
Since we generate the coordinates of the contour lines with the same sampling rate, used in ellipsoid, an equal amount of memory is required. We set the length of an SWP segment at r 0 = 1.51m and locate the starting point ⃗ P 1 in zone z 4 , which is bounded by the contour lines h = −1653m and h = −1854m, so that the starting point is very close to the position, used to evaluate the error of D8 algorithm [4]. The distance between the starting points and the inverted ellipsoid center is set at 900 meters. We select three starting points along a circle with a radius of 900 meters to delineate three extracted SWPs in three directions, namely 135, 60, and 270 degrees, where the extracted SWPs are labeled by "Line-1", "Line-2", and "Line-3", respectively and are plotted with blue markers in Fig. 8. Fig. 11 presents the position error E(%) versus the distance of the SWP flowing from the border of inverted ellipsoid to the center. The position error E(%) is calculated with (46). First, the positions of the "Line-1" extracted SWP are compared with θ Exact = 135 degrees, and their error E(%) are below 0.1 and are plotted in Fig. 11 (a). Also, the positions of the "Line-2" and "Line-3" extracted SWPs are compared with θ Exact = 60 and θ Exact = 270 degrees, respectively, and their error E(%) are below 0.0002 and 0.00015, and are plotted in Fig. 11 (b) and (c), respectively.

4) Extracted SWPs over inclined plane
Unlike ellipsoid and inverted ellipsoid, the formula of the inclined plane listed in Table 1 is used to generate only one incomplete contour line which is in a square shape with the dimension 900m × 900m. The four corner points along the incomplete contour have different elevation levels, namely 0, 270, 630 and 360 meters, respectively. The sampled points are plotted in Fig. 9. Then, a boundary C in the numerical solution in (3) consists of 4 intervals (C j ).
After collecting the information of the boundary into the memory, we set the length of an SWP segment at r 0 = 20.0m. To apply the method of [4], we put the starting points very close to the position, used to evaluate the error of the D8 algorithm. Hence, we put the starting points at the coordinates (−360, 600), (240, −540) and (600, 150) to delineate the "Line-1", "Line-2" and "Line-3" extracted SWPs in the

5) Comparison to other algorithms
In previous section, the position errors of the delineated SWPs compared to the theoretical SWPs are calculated by using the evaluation method in [4], and found that each surface shape requires different resolution to achieve small position errors. For example, an ellipsoid and an inverted ellipsoid require higher resolution than a inclined plane. In this section, we present the accuracy comparison between the proposed algorithm and existing algorithms, including D8, Rho-8, D8-LTD, FDFM, FMD-md and D∞. For fair comparison, the proposed algorithm is tested with the same set of surface shapes available in [4], namely, ellipsoid, inverted ellipsoid and plane. Table 2 shows the average of position errors of the proposed algorithm in Fig. 10 to Fig. 12 and the average of position errors of the existing algorithms.
The proposed algorithm significantly outperforms all existing algorithms. In the ellipsoid case, the average position error of the proposed algorithm is 0.120% whereas the best existing algorithm (D8-LTD) obtains 47.5%. In the inverted ellipsoid case, the average position error from the proposed algorithm is 0.120% whereas the best existing algorithm (D∞) obtains 88.8%. In the plane, the average position error from the proposed algorithm is 0.001% whereas the best existing algorithm (D8-LTD) obtains 48.6%. Furthermore, while increasing the resolution of input data from 90 to 180 does not improve the existing algorithms [4], it helps the proposed algorithm achieve much smaller average position error in case of ellipsoid and inverted ellipsoid.

B. NATURAL TERRAIN
This work considers the real terrain, referred to as the Bang Wad area (Phuket, Thailand). This area often lacks the water supply for the tourism demand. Therefore the water management system must be improved. The center of the Bang Wad has latitude 7.88708 • N and longitude, 98.33433 • E. In this area, the geologic substratum is constituted of igneous rock. The topography can be described as fairly complex as a result of surface water erosion and active landslide. To make a quantitative analysis of the extracted SWPs on real terrains, a 30 × 30 DEM matrix is chosen and is downloaded from the Earth Science Data Systems (ESDS) (https://earthdata.nasa.gov/learn/articles/newaster-gdem) with a horizontal resolution of 30 meters and a vertical resolution of 1 meter. The elevation ranges from 515.0 to 45.0 meter amsl (above mean sea level). Fig. 13 presents a considered area of 5000 × 10000 m 2 in which elevations are illustrated with contour lines from 500.0 to 100.0 meter amsl (above mean sea level). Note that the gray contour lines are at an elevation level of below 300 meters and are very long.
To validate that the proposed algorithm can delineate the SWPs with the incomplete contour lines, we compare the obtained results with the SWPs from the complete contour lines. First, we select an area which has the SWPs flowing down from a top hill to the Bang Wad Dam. The elevation levels of this area are presented with contour lines from 500 to 340 meters. Second, we determine the contour lines' resolution that is appropriate for our proposed SWP delineation algorithm. The appropriate resolution will make the delineation fast and accurate. By trial-and-error, a contour line should contain the discrete elevation data at every 20 meters, and a contour interval between two adjacent contour lines should be 20 meters or below. For the sake of computation time and efficiency, we choose a contour line interval of 20 meters as depicted with green lines in Fig. 13. Third, some contour lines are converted into incomplete contour lines. The contour lines that are not crossed over by the SWPs flowing down from the top hill to the Bang Wad Dam are converted into incomplete contour lines. Fig. 14 (a) presents the complete contour lines which are the input data of the proposed algorithm depicted with the green lines. To delineate the SWPs from these complete contour lines, we employ Algorithm 1, which determines the SWP direction in steepest descent with θ min , and let every SWP segment have an equal length of r 0 = 0.5 meter. In Fig. 14 (b), we choose the extracted SWPs' starting points ⃗ P 1 , which are labeled on the red line between the contour lines of 500 and 480 meters. Then, the proposed algorithm delineates the continuous SWPs, some of which converge (labeled with "C1" and "C2") while the others diverge (labeled with "D").
Next, we convert the complete 20-m contour lines in Fig. 14 (a) into the incomplete contour lines in two different types, namely type-1 and type-2 which are plotted in Fig. 14  (c) and (e), respectively. Since we focus on the SWPs located over the right-hand side of the hill, we do not make any change to these contour lines and shorten the length of the contour lines on the left-hand side of the hill. The type-1 incomplete contour lines in Fig. 14 (c) comprise only two incomplete contour lines at the elevation of 360 * and 340 * meters. The type-2 incomplete contour lines in Fig. 14 (e) comprise all incomplete contour lines at the elevation between 500 and 340 meters. Then, we use previous configuration in Algorithm 1 to delineate the SWPs in cases of the type-1 and type-2 incomplete contour lines, which are shown in Fig. 14 (d) and (f), respectively. The obtained results point out that the SWPs' positions in Fig. 14 (d) are exactly the same as the reference which are the extracted SWPs from the complete contour lines in Fig. 14 (b). Moreover, the SWPs are converging at the same points, which are labeled with "C1" and "C2", and are diverging at the same points, which are labeled with "D". Likewise, in Fig. 14 (f) the converging SWPs (which are labeled with "C1" and "C2") and diverging SWPs (which are labeled with "D") are exactly the same as the delineated SWPs from the complete contour lines. Therefore, the proposed algorithm maintains the accuracy of the delineated SWPs even though the input contour lines are incomplete. To demonstrate the advantage of the proposed algorithm in terms of computation time, we measure the computation time of each mentioned case. The computation time depends on i) the number of zones, ii) the number of the contour lines, iii) the number of the sampled points along a contour line, and iv) the length of the SWP segment. To analyze the performance of the proposed algorithm, we compare computation times in three cases, namely complete contour lines, type-1 incomplete contour lines and type-2 incomplete contour lines, which are shown in Fig. 14 (a), (c) and (e), respectively. In each case, we designate a catchment area of nine zones and a position as the starting point for all cases. The calculation is in the environment with Intel(R) Xeon(R) processor (E3-1225 v5) 4-core 3.7GHz/8MB cache, 8GB memory and Fedora Linux (Workstation) version 35. The used compiler is GCC 11.2.1 (Red Hat 11.2.1-9). In Table  3, the proposed algorithm can reduce the total computation time from 21 to 3 minutes when using the type-2 incomplete contour lines.
All along, the SWPs are delineated from upper contour line to lower contour line. In fact, the proposed algorithm can also delineate SWPs in the reverse direction, that is, from lower contour line to upper contour line as shown in the following example. First, the input data are chosen. We choose the type-2 incomplete contour lines as the input data to demonstrate the robustness of the proposed algorithm as shown in Fig. 15. Second, we configure Algorithm 1 with θ max instead of θ min to guide the SWP direction in steepest ascent. Third, we put the starting points along the twelve sink mouths which are depicted along the purple line in Fig. 15. Executing the Algorithm 1, we obtain the SWPs which are depicted with the blue lines in Fig. 15. Note that the sink mouths depicted with "C1" and "C2" have the SWPs which are surrounding wider area than the sink mouths labeled with VOLUME 4, 2016 of delineating reverse SWPs is to identify the boundary of the catchment area. The size of the catchment area is an essential parameter in the hydrological model [42]. It is used to calculate the water flow rate and evaluate the flood risk level.

VII. CONCLUSION
This paper introduces a new contour-based algorithm for delineating the SWPs from a DEM image. Unlike conventional contour-based algorithms, which demands computation resources, the proposed algorithm requires less computation resources while maintaining the position accuracy of the delineated SWPs. The semi-analytical solution of Laplace's partial differential equation is used to identify the SWP direction in either the steepest downslope or steepest upslope. With the formulated semi-analytical solution based on BEM, the problem domain is discretized along the contour lines instead of the area surrounding the contour lines. Therefore, the domain is parameterized by one variable instead of two variables. Also, the possible SWP direction ranges from 0 to 360 degrees. To validate the proposed algorithm, the method to measure the SWP position error, compared with the theoretical SWPs, in [4] is used here. The standard synthetic surfaces, where theoretical SWPs exist, namely ellipsoid, inverse ellipsoid, and inclined plane, were used for validation. The position error E(%) versus the distance of the SWP over the ellipsoid, the inverse ellipsoid, and the inclined plane are all below 0.0002, and are less than the error of any grid-based algorithms, which are in the ranges from 16.3 to 75.2. Moreover, the semi-analytical solution is robust to BEM with incomplete contour lines, which demand even less computation resources. The real terrain with dam and mountains, where incomplete contour lines exist due to the very long contour lines at foothills, were used to show the robustness. In doing that, we compare the SWPs that are delineated from incomplete contour lines, which are shortened differently.
All obtained SWPs are smooth all the way from the upper contour line to the lower. Most importantly, they converge identically when the length of the SWP segment is set at or less than 0.5 meter, and the contour line interval is 20 meters or less, and the spacing between adjacent discrete elevation data is 20 meters or narrower. Also, they diverge identically.
One of the main practical advantages, the proposed algorithm can also delineate SWPs from the lower contour line to the upper. In other words, SWPs are delineated in the reverse direction of water flow. The delineated SWPs can be used to identify the catchment area's boundary that is beneficial for the catchment's estimation in the hydrological model to analyze the water flow rate. However, we found that, in the phase of determining the direction of the SWP segment, the criteria for choosing the direction of the steepest ascent is less accurate when the starting points of the SWP segment are designated inside the inappropriate area; for example, the center of the area where the spatial distance between adjacent contour lines exceeds 50 meters. Therefore, we should be aware of this limitation when applying the proposed algorithm to estimate the catchment area. .

APPENDIX A APPROXIMATION OF GRADIENT
The numerical solution of ∇ P (h( ⃗ P )) is where: The F ′ αj , F ′ βj , F ′ γj , F ′ ηj , F ′ ζj and F ′ ξj are geometric terms, as known as shape function. Their values do not depend on the elevation function and its gradient. They depend only on the shape of the boundary and the position of the point ⃗ P , and can be calculated after the boundary and point ⃗ P were specified.

APPENDIX B APPROXIMATION OF GRADIENT OF GRADIENT
The numerical solution of ∇ 2 P (h( ⃗ P )) is The F ′′ αj , F ′′ βj , F ′′ γj , F ′′ ηj , F ′′ ζj and F ′′ ξj are geometric terms, as known as shape functions. Their values do not depend on the elevation function nor its gradient. They depend only on the shape of the boundary and the position of the point ⃗ P and can be calculated once the boundary is specified and the point ⃗ P is chosen.

APPENDIX C APPROXIMATION OF SURFACE
If the length of the SWP segment is short enough, the acceptable approximation can be obtained using only the constant, linear, and quadratic terms: h(r, θ) ≈ a 0 + r(a 1 cos θ + b 1 sin θ) + r 2 (a 2 cos 2θ + b 2 sin 2θ); r < r 0 (61) Since the longest length of the SWP segment that still gives an acceptable approximation depends on the shape of the surface, r 0 is chosen by trial-and-error.
To fit the quadratic approximation in (61), we select the value of the coefficients a 0 , a 1 , a 2 , b 1 and b 2 for r < r 0 . Usually, these coefficients are not unique since they depend on the criterion that we adopt to measure the difference between the surface and the approximation. Here, we adopt the criterion, where the approximation error at r = 0 is zero.
Since the approximated surface h(0, θ) equals the elevation h(x, y) at point ⃗ P and the series expansion has only a 0 left when r = 0, we obtain h(0, θ) = a 0 = h( ⃗ P (x, y)).
where the polar coordinates r and θ can be converted to the Cartesian coordinates x and y by using the trigonometric functions, namely sine and cosine: x = r cos θ and y = r sin θ, and h( ⃗ P ) can be determined from (27). Next, we get rid of the term with a 0 to find coefficients a 1 and b 1 by taking first derivative to (61), and get rid of the term with a 2 and b 2 by substituting r = 0, we obtain h r (0, θ) = h x ( ⃗ P )x r + h y ( ⃗ P )y r a 1 cos θ + b 1 sin θ = h x ( ⃗ P ) cos θ + h y ( ⃗ P ) sin θ, (63) where, ∂x/∂r = cos θ and ∂y/∂r = sin θ. When substituting θ = 0 and θ = π/2 into (63), we obtain the coefficients a 1 and b 1 , respectively, that is where the coefficient a 1 and b 1 are fit by some components of the gradient of the elevation ∇ P h( ⃗ P ), which can be determined from (32). Then, we take the second derivative to the elevation h(x, y) at point ⃗ P with respect to r. We obtain h rr (0, θ) = (h x ( ⃗ P )x r + h y ( ⃗ P )y r ) r = h xr x r + h x x rr + h yr y r + h y y rr = (h xx x r + h xy y r )x r + h x x rr + (h yx x r + h yy y r )y r + h y y rr = h xx (x r ) 2 + 2h xy y r x r + h yy (y r ) 2 + h x x rr + h y y rr When substituting ∂ 2 x/∂r 2 = ∂ 2 y/∂r 2 = 0, we obtain h rr (0, θ) = h xx cos 2 θ + h yy sin 2 θ + 2h xy cos θ sin θ Since the second derivative of the approximated surface in (61) equals h rr in (67), we obtain 2(a 2 cos 2θ + b 2 sin 2θ) = h xx cos 2 θ + h yy sin 2 θ + 2h xy cos θ sin θ.
When substituting θ = 0 and θ = π/4 into (68), we obtain the coefficients a 2 and b 2 , respectively, that is To satisfy (1), h xx ( ⃗ P ) + h yy ( ⃗ P ) equals zero. Therefore, where the coefficient a 2 and b 2 are fit by some components of the gradient of gradient of the elevation ∇ 2 P h( ⃗ P ), which can be determined from (41). Now, we obtain the approximated surface over the neighborhood of the starting point of the SWP segment h(r, θ) ≈ h( ⃗ P ) + r h x ( ⃗ P ) cos θ + h y ( ⃗ P ) sin θ ; r < r 0 . + r 2 1 2 h xx ( ⃗ P ) cos 2θ + 1 2 h xy ( ⃗ P ) sin 2θ (71)

APPENDIX D DETERMINATION OF THE STEEPEST DESCENT
To find the θ min , we differentiate the trace on the circle of the approximated surface, that is dh dθ (r 0 , θ) = r 0 [−h x sin θ + h y cos θ] + r 2 0 [−h xx sin 2θ + h xy cos 2θ] .
Then, we equate the derivative to zero, that is h x sin θ + r 0 h xx sin 2θ = h y cos θ + r 0 h xy cos 2θ. (73) To solve this equation, trigonometric identities are applied to convert the equation into a fourth order (quartic) equation with either cos θ term or sin θ term. First, substituting sin 2θ = 2 sin θ cos θ and cos 2θ = 2 cos 2 θ − 1 into (73), as well as replacing cos θ by ν and sin θ by ω, we obtain (h x + 2h xx r 0 ν) ω = h y ν + h xy r 0 2ν 2 − 1 . (74) Next, since ν 2 + ω 2 = 1, ω is replaced by √ 1 − ν 2 , that is, (h x + 2h xx r 0 ν) 1 − ν 2 = h y ν + h xy r 0 2ν 2 − 1 . (75) Last, taking square to both sides of the equation to get rid of the square root, we manipulate into where the coefficients c i equal The standard techniques can be used to solve the roots of the quartic equation in (76), and the four (possible complex) roots: ν 1 , ν 2 , ν 3 , ν 4 are obtained. Next, substituting these roots into (74) gives the four corresponding values: ω 1 , ω 2 , ω 3 , ω 4 . Finally, θ min can be calculated using the ν, ω pairs and the arc-tangent function. Note that θ min are real number when the ν, ω pairs are real numbers with magnitudes ≤ 1. Pairs, which are complex numbers or have magnitudes greater than 1, should be discarded.