GIS-based watershed characterization and morphometric analysis in Bohol Watersheds, Philippines

ABSTRACT

The study presents a GIS-based morphometric analysis characterizing the watersheds in Bohol, Philippines. An IfSAR-derived DEM data were processed and used to delineate the watersheds and evaluate morphometric parameters in QGIS/QSWAT. The analysis revealed that Loboc Watershed (LW) and Wahig-Inabanga Watershed (WIW) were sixth-order watersheds. Covering 675 km² (LW) and 610 km² (WIW), the watersheds encompassed about 34% of the mainland area. The bifurcation ratios, 3.91 (WIW) and 4.06 (LW), suggested highly dissected watersheds with high flood peaks. LW has low permeability, and high runoff and erosion susceptibility compared to WIW, as indicated by higher values for the form factor, circulatory ratio, drainage density, and stream frequency. The ruggedness numbers of 1.653 (WIW) and 1.457 (LW) indicated steep slopes in WIW than in LW, associated with high erosion and environmental degradation. The hypsometric integral of 0.22 (WIW) and 0.37 (LW) implied lesser erosional areas in WIW than in LW. The study found that both watersheds were highly susceptible to erosional processes. The two watersheds performed critical environmental functions on Bohol Island; hence, the results of this study present a piece of baseline information for resource managers to conduct further investigations that would identify hotspot areas and enable targeted soil and water conservation initiatives.

KEYWORDS:

• Landform dynamics
• erosion susceptibility
• runoff characteristics
• watershed management

I. Introduction

Philippine watersheds provide environmental services to downstream communities (Francisco & Rola, 2004) by supplying water for crop production, domestic use, hydroelectric power generation, or a mix of these. However, adversities also occur in watersheds in response to climate extremes and tectonic activities. These natural responses of a hydrologically defined landform are described and better understood by analyzing surface features in the study of watershed morphology (Beven et al., 1988; Kabite & Gessesse, 2018; Zavoianu, 1985).

Watershed morphometric analysis is crucial in providing insights into the connections between landform dynamics and the hydrologic responses within the unit. A watershed’s morphologic feature is a function of its topography and the form and extent of its stream system or drainage network (Horton, 2004). A morphometric analysis measures the mathematical properties of watersheds and stream networks, including the landform’s length, surface area, size, and slope. These parameters are significant determinants of water dynamics within watersheds which can give valuable knowledge for land and water resource management in regions with limited hydrological data (Meshram & Sharma, 2020). Horton’s efforts were crucial in developing the analytical method for assessing the watershed’s morphological characteristics (Strahler, 1957).

Flood forecasting and hazard studies (Alam et al., 2021; Charizopoulos et al., 2019; Taha et al., 2017) are a few of the applications of watershed morphometric analyses. Within a geographic information system (GIS) platform, spatial datasets such as the digital elevation model (DEM) are the primary input for delineating watersheds and extracting morphometric parameters. Flood zoning studies in southeastern Bangladesh (Alam et al., 2021) and in the Wadi Qena basin, Eastern Desert of Egypt (Taha et al., 2017), integrated morphometric analysis. Alam et al. combined the morphometric parameters of a watershed with local topographic effects, such as the topographic wetness index and the topographic position index, to map flash flood susceptibility areas. In a similar study by Taha et al. (2017), high flood hazard catchments have low bifurcation ratio, high stream frequency, and high drainage density values. Flood hazards were investigated and linked with tectonic faulting on Samos Island in Greece (Charizopoulos et al., 2019). The study by Charizopoulos et al. (2019) observed that where fault orientation and flow directions differ, tectonics had little effect on the evolution of the drainage systems. Despite data constraints to carry out detailed flood assessment, morphometric analysis in flood hazard studies produced valuable information to comprehend watershed flood hydrology better for developing mitigation strategies in the study areas.

Morphometric parameters are proxies of the erodibility status of a watershed. Several studies demonstrate the application of morphometric analysis to identify and prioritize areas for land and water conservation initiatives and other natural resources management efforts (Benzougagh et al., 2022; Bharath et al., 2021; Fenta et al., 2017; Puno & Puno, 2019). Watershed morphometric parameters, derived and analyzed using a digital elevation model in a GIS environment, were grouped into linear, aerial, and relief aspects. The linear and relief parameters directly relate to erodibility, while the aerial parameters indirectly associate erodibility (Nookaratnam et al., 2005). The ranking of parameter values, classification and prioritization of watersheds and sub-watersheds, as presented in Nookaratnam et al. (2005), is commonly adopted in several studies. Other methods, such as hierarchical cluster analysis and principal component analysis, were presented in the study of Kulimushi et al. (2021). The hierarchical cluster analysis grouped the watersheds and sub-watersheds based on morphometric characteristics. The principal component analysis reduced the number of morphometric parameters by identifying several factors that explain most of the variance observed in the parameters considered. Puno and Puno (2019) included landuse/landcover information to arrive at compound values in the final prioritization of areas. In the studies mentioned, morphometric analysis in erosion assessment across watersheds enables the development of maps showing spatial erosion variability and priority areas. These maps serve as a basis for improved planning, management, and decision-making to ensure the sustainable use of watershed resources.

Morphometric analysis studies also helped comprehend the evolution of river basins in response to climatic and tectonic activities (Mahadev et al., 1947; Singh et al., 2021). A study by Mahadev et al. (1947) combined morphometric analysis with OSL dating of sediments collected from rivers to help understand climatic and tectonic activities with the evolution of watersheds and drainage systems. Morphometric analysis and luminescence chronology established channel shift time frames in response to tectonic activity (Mahadev et al., 1947). The study of Singh et al. (2021) quantified various aspects of a river basin employing a morphometric analysis.

As demonstrated in the cited studies, the morphometric analysis provided a rapid but scientific way of assessing watershed flood hazards and erodibility status and the connections between tectonics and climate factors. Several studies have shown that with GIS’s advanced capability and spatial datasets’ availability, watershed characterization using morphometric analysis produced datasets critical for decision-makers and resource managers. Hence, in the present study, a morphometric analysis is carried out to assess the watersheds of Bohol Island in the central Philippines, where such an approach has not been applied. This study is limited to using a digital elevation model to extract watershed morphometric parameters and decipher surface processes as initial steps to watershed characterization activities. The information generated will contribute to a better understanding of watershed hydrology in support of sustainable resource management efforts in the study area.

II. Materials and methods

Description of the study area

Figure 1 depicts the Visayas region of the Philippine archipelago, where Bohol Island is located (A). Bohol is the tenth-largest island in the country, covering about 4,110 km², where 66% of the area is devoted to agriculture-related activities (PPDO, 2006). The overall land area has an existing forest cover that accounts for 24.6% of the total.

Figure 1. The geographic location of the study area. A) The major watersheds, drainage networks, and fault lines overlay on an elevation map. B) The mouth of Loboc Watershed. C) The midstream agricultural area of Wahig – Inabanga Watershed.

The island lies on several Cretaceous oceanic lithosphere formed and deposited in various geological settings and patterns (Dimalanta et al., 2006). The geology of Bohol island (JICA, 1985) is composed of volcanic rock in the eastern and western parts and a sizable area of sedimentary rock in the middle of the island known as the Ilihan formation. The formation is a relatively broad sequence of well-bedded sandstone, shale, and mudstone/siltstone. These units often consist of volcanic buildup formed by weathering nearby country rocks. About 30% of the province’s rocks are intrusive and extrusive igneous rocks, particularly in the island’s northeastern region. Uplifted reef limestone facies of significant height surrounded the perimeter of mainland Bohol. On the island, sedimentary rocks make up roughly 70% of the terrain.

The island has two major drainage networks: the Loboc Watershed (LW) and the Wahig-Inabanga Watershed (WIW). The significant economic activities in these watersheds, primarily tourism and agriculture, highlight the importance of studying landform dynamics concerning resource conservation and sustainable use. LW hosts the Loboc River (Figure 1B), where the cruise is a tourist attraction. At the same time, WIW supports crop production activity in the island province (Figure 1C). LW and WIW, encompassing 675 km² and 610 km², respectively, cover about 34% of the mainland area. The two drainage systems, with contrasting geometry, drain the island in different directions, as shown in Figure 1A.

Water resources in LW are indispensable to tourism and energy production. In WIW, the resource is vital to irrigation for agriculture. However, natural catastrophes such as flooding, erosion, and drought are commonplace. These are significant concerns to the local government units and the people. Typhoons and other weather disturbances frequently visit the country, with an average of 20 storms per year. Generally, according to the modified Coronas classification, the island’s climate is under Type IV, characterized as having rainfall more or less evenly distributed over the whole year. The mean annual rainfall from 2000 to 2012 is approximately 1653 mm, as recorded in the provincial capital of Tagbilaran City. The average monthly precipitation records show that April and December are the driest and wettest months, with 78 mm and 185 mm of precipitation, respectively. The island’s temperature ranged between 24°C and 32°C, while the relative humidity averaged 83%.

Bohol is a tectonically active island. Figure 1A shows two prominent fault lines running through the province, the East Bohol Fault and the newly discovered North Bohol Fault (Kobayashi, 2014). The faults’ presence raises the risk of tectonic hazards such as earthquakes and earthquake-induced hazards (e.g., tsunamis and landslides). A combined tectonic movement and eventual land erosion shaped the island in the last 150 million years. Although these features provide vital information in the watershed characterization process, the current study is limited to topographic analysis using a digital elevation model as a preliminary step.

Data sources and processing

The morphometric analysis reported in this study uses an IfSAR-derived digital elevation model (DEM) from the National Mapping and Resource Information Authority (NAMRIA), the mapping agency of the Philippines. The DEM has a spatial resolution of 5 m with 1.0-m root-mean-square error vertical accuracy and a 0.625-meter resolution ORI (orthorectified radar image) at 2.0-m horizontal accuracy (NAMRIA, 2021, page 12). The study also adopts the watershed delineation process of the QSWAT model to extract the drainage network, the watersheds, and the geometric attributes. Figure 2 shows the process flow of the analysis.

Figure 2. The process flow of morphometric analysis in Bohol watersheds.

The DEM preprocessing consists of projecting the data to a local projection system and filling the sinks to remove spurious areas to produce sink-free elevation data. Watershed delineation uses the QSWAT software (https://swat.tamu.edu/software/qswat) and follows the process presented in the QSWAT User manual available on the SWAT website. Briefly, watershed delineation involves defining a threshold area to create streams and channels. Once the channels and streams are defined, QSWAT delineates the watershed boundaries using user-identified outlet points or drainage locations. The delineation process also includes assigning a stream order to channel and stream, calculating the stream and watershed geometry, and landscape analysis. A text file output stores the attributes calculated in the landscape analysis. Following the equations in Table 1, the morphometric analysis proceeds using a spreadsheet and visualization in QGIS software.

Table 1. List of morphometric parameters used in the study, equations, and the source of information.

Morphometric parameters (Notation)[Unit]	Equation	References/Source
Linear aspects
Stream order (U)	Hierarchical rank	(Strahler, 1952)
Stream number (Nu)	Nu = N1+ N2 + …Nn	(Horton, 1932)
Stream length (LT)	LT = L1+ L2 + …Ln	(Horton, 1932)
Mean stream length (Lm)	Lm = Lu/Nu	(Strahler, 1957)
Stream length ratio (RI)	RI = Lu/L(u-1)	(Horton, 1932)
Bifurcation ratio (Rb)	Rb = Nu/N(u+1)	(Strahler, 1957)
Mean bifurcation ratio (Rbm)	Rbm = average of Rb	(Horton, 1932; Schumm, 1956)
Watershed length (Lb) [km]	GIS
Main stream length (LM) [km]	GIS
Areal aspects
Watershed area (A) [km^2]	GIS
Perimeter (P) [km]	GIS
Form factor (Rf)	Rf = A/(Lb)^2	(Horton, 1932)
Elongation ratio (Re)	Re = 2/Lb*(A/π)^1/2	(Schumm, 1956)
Circulatory ratio (Rc)	Rc = 4πA/P^2	(Schumm, 1956)
Drainage density (Dd)	Dd = Lu/A	(Horton, 2004)
Drainage texture (Dt)	Dt = Nu/P	(Horton, 1932)
Stream frequency (Fs) [km^−2]	Fs = Nu/A	(Horton, 2004)
Compactness coefficient (Cc)	Cc = 0.2841 × (P/A)^0.5	(Horton, 1932)
Relief aspects
Maximum (Emax)	GIS
Minimum (Emin)	GIS
Mean (Emean)	GIS
Relief (R)	R = (Emax – Emin)	(Strahler, 1952)
Relief ratio (Rr) [m/km]	Rr = H/Lb	(Schumm, 1956)
Relative relief (Rhp)	Rhp = (Emax × 100/ P)	(Schumm, 1956)
Ruggedness number (Rn)	Rn = Dd* R/1000	(Asfaw & Workineh, 2021)
Hypsometric curve	GIS
Hypsometric integral (HI)	HI = (Emean – Emin)/(Emax – Emin)	(Strahler, 1952)

Evaluation of morphometric parameters

The morphometric parameters, grouped into linear, areal, and relief aspects (Asfaw & Workineh, 2021; Bharath et al., 2021), are given in Table 1 with the corresponding equations to quantify the values. These parameters directly or indirectly affect runoff characteristics and erosion processes, so they can be used to characterize watersheds.

III. Results and discussion

Linear aspects of LW and WIW

A watershed’s linear aspects comprise the streams and the stream network and are typically one-dimensional quantities (Asfaw & Workineh, 2021). Table 2 lists the linear parameter values of the study watersheds.

Table 2. Stream network characteristics and the linear aspects of LW and WIW.

Parameters	LW								WIW
U	1	2	3	4	5	6	-	1	2	3	4	5	6	-
Nu	951	304	56	11	3	1	1362	852	194	48	10	3	1	1108
Lu	0.7	1.1	2.9	7.0	17.6	39.2	-	0.6	1.4	3.4	6.6	7.7	85.0	-
LT	666	347	160	77	53	39	1342	557	264	163	66	23	85	1158
Rl,u	-	0.52	0.46	0.48	0.69	0.74		-	0.47	0.62	0.40	0.35	3.69
Rl,m	0.58							1.11
Rb	3.13	5.43	5.09	3.67	3.00	-	-	4.39	4.04	4.80	3.33	3.00	-	-
Rbm	4.06							3.91
Lw	38.6							47
LM	88							105

U – Stream order; N_u – Stream number; L_u – stream length [km]; L_T – Total stream length [km]; R_l,u – Stream length ratio; R_l,m – mean stream length ratio; R_b – Bifurcation ratio; R_bm – Mean bifurcation ratio; L_w – Watershed length [km]; L_M – Main stream length [km]; LW – Loboc Watershed; WIW – Wahig–Inabanga watershed

The stream order (U) is a metric used to describe the position of a stream within its tributary hierarchy. According to Horton (1932), U is a quantitative basis for comparing the degree of development of the drainage networks of watersheds of similar size. Strahler’s (1957) stream ordering scheme applied in the study assigns an order to the streams in this way: The first orders are the smallest finger type and non-branched tributaries; the second orders streams formed when two first stream orders meet; the third orders formed when the second orders meet, and so on. The analysis found the two watersheds as sixth order and generally of dendritic patterns, as shown in Figure 3. A drainage network with a dendritic pattern suggests textural homogeneity and a lack of structural control (Christopher et al., 2010).

Figure 3. Drainage network organization and the stream orders in the watersheds under study.

The stream number (N_u) refers to each order’s total count of stream segments, providing information on surface runoff characteristics. According to Horton (1932), N_u of a given stream order U forms an inverse geometric series: N_u decreases with increasing U. The analysis reveals that LW has 12% more first-order streams than WIW. Further, LW has about 20% more streams on the total count of streams than WIW, implying a better drainage system but higher chances of erosion and a lower potential for groundwater recharge. As Saha et al. (2022) noted, a watershed with numerous first-order streams has a high probability of flash floods during heavy rainfall than a watershed with fewer first-order streams.

When other factors are held constant, the time required for the water to flow a given distance is directly proportional to the stream length (Langbein, 1947). Stream length (L_u) refers to the mean stream length in a given stream order U. As noted in Table 2, L_u increases as U increases which means that the L_u of the first-order stream is lower than the L_u of the next higher-order stream. Also noted is the considerably more extended length of the sixth-order stream in WIW, which is more than twice the length of LW’s sixth order stream. Watersheds with longer stream lengths tend to be flatter, while watersheds with short stream lengths are steeper (Christopher et al., 2010; Kabite & Gessesse, 2018). Table 2 also revealed that stream orders 2, 3, 4, and 6 of LW are shorter than the length of the same stream orders of WIW, suggesting steeper and faster runoff in LW. Longer stream lengths tend to have higher chances of infiltration and slower runoff velocity, thus reducing surface erosion. The sum of all stream lengths (L_T) in a given order is inversely proportional to U: L_T decreases with increasing U. Asfaw and Workineh (2021) observed the same pattern in their empirical study.

The stream-length ratio (R_l,u) is the ratio of the mean stream length of a given stream order (L_u) to the mean stream length of the next lower stream order (L_u-1; Horton, 1932). The mean stream-length ratio (R_l,m) is a watershed’s average R_l,u.R_l,u changes with U, as shown in Table 3, where such change, according to Mahala (2020), is related to geomorphic development stages and the irregularity of hydrological behaviour. The R_l,m of LW is lower than the R_l,m of WIW, suggesting that LW is at an earlier stage of geomorphic development than WIW, as implied in the study of Mahala (2020).

Table 3. Results of the areal aspect analysis in LW and WIW.

Watersheds	A	P	Rf	Re	Rc	Dd	Dt	Fs	Cc
LW	676	164	0.454	0.552	0.315	1.99	7.77	1.89	1.78
WIW	610	203	0.276	0.556	0.186	1.90	5.57	1.85	2.32

A – area [km²]; P – perimeter [km]; R_f – form factor; R_e – elongation ratio; R_c – circulatory ratio; D_d – drainage density; D_t – drainage texture; F_s – stream frequency [km⁻²); C_c – compactness coefficient

Horton’s bifurcation ratio (R_b) is the ratio of the stream number of any order to the stream number of the next higher order (Strahler, 1957). In Horton’s (1932) study, R_b values range from about 2 for flat or rolling and 3–4 for mountainous or highly dissected drainage basins. Hence, R_b is higher for hilly and well-dissected watersheds than for rolling watersheds. Further, Fitzpatrick (2017) indicated that R_b values falling between 2 and 4 indicate runoff peakedness. Low ratios are associated with flashy hydrographs with high flood peaks for small watersheds (Fitzpatrick, 2017). From the R_b values in Table 3, both LW and WIW can be described as highly dissected and mountainous and could potentially experience high flood peaks. Within the R_b values of 3.0 to 5.0, according to Strahler (1957), the geological structures of LW and WIW cannot be associated with stream pattern development.

Watershed length (L_b) is the dimension that runs parallel to the primary drainage channel. The maximum watershed length is close to twice the diameter of a circle having the same area as the watershed for circular-shaped catchments (Schumm, 1956). L_b reveals the watershed’s main channel, which carries most water, and is commonly used to compute the time element of water flow over a watershed. Longer L_b indicates more time for surface water to reach the mouth of the watershed. Using GIS, this study measures the watershed length as the distance between the mouth and the most distant point in the watershed. Using this approach, the length of WIW is 8.4 km longer than LW. Hence, runoff water from WIW would travel longer than the runoff in LW.

The mainstream length (L_M) is the longest flow path in the watershed. The analysis revealed that WIW longest flow path is about thrice of the L_M of LW. One of the consequences of a longer runoff distance is the high tendency of runoff to infiltrate the ground, recharging groundwater sources because of the longer travel time.

Areal aspects of the LW and WIW

The areal aspects or shape indices of watersheds, given in Table 3, are associated with drainage characteristics and hydrological properties. The watershed area (A) is the upslope area that contributes flow to an identified reference point. In the study, the reference point refers to the outlets or mouths of the watersheds. A watershed’s area directly reflects the water volume that will flow or pass through a reference point. According to Langbein (1947), the discharge volume in a humid climate varies directly with the size of the tributary drainage area. As depicted in Table 3, LW has about 10% more watershed area than WIW, indicating more runoff from LW.

A watershed’s perimeter (P) refers to the outer boundary’s total length. P indirectly influences the drainage texture, relative relief, and circulatory ratio and directly affects the compactness coefficient. As P increases with increasing catchment area, the geometry of the enclosed area dictates its magnitude. As shown in Table 3, LW’s area is 10% bigger than that of the WIW, but its perimeter is about 24% shorter than that of the WIW.

The form factor, elongation ratio, and circulatory ratio are indicators of watershed shape, volume, and the duration of runoff. The form factor (R_f; Horton, 2004) is the ratio of the watershed area divided by the square of the length. R_f explains a watershed’s flow property (Bogale, 2021) as affected by shape. In Table 3, the F_f values of the watersheds suggest the elongated shape of WIW and the round shape of LW. Low F_f values reflect an elongated shape, while high F_f values indicate rounded-shaped watersheds. High values of F_f are associated with increased runoff and short-time concentration and are highly sensitive to flooding (Bogale, 2021). The elongation ratio (R_e) is the ratio of the diameter of a circle having the same area as the watershed to the watershed length, which varies from 0.6 to 1.0 (Fitzpatrick, 2017). R_e values between 0.6 and 0.8 are associated with high relief and steep ground slopes (Schumm, 1956). R_e value closer to 0.6 is a character of an elongated basin with slower runoff and less flashy floods, while a ratio of 1.0 is a round watershed (Fitzpatrick, 2017). Table 3 shows similar R_e values of the two watersheds, which indicate similar behaviour. The circulatory ratio (R_c) is the ratio of the watershed area to the area of the circle having the same circumference as the perimeter of the watershed (Miller, 1953). R_c for elongated watersheds approaches zero, while R_c for circular watersheds approaches 1. R_c values also indicate the shape of a watershed where a circular to oval shape allows quick disposal of runoff resulting in a high-peaked and narrow hydrograph. In contrast, the elongated shape allows slow water removal, resulting in a broad and low-peaked hydrograph (Pandey et al., 2005). The circular shape of LW with R_c of 0.315 suggests quick and high-peak runoff, which encourages soil erosion, while an elongated shape WIW (R_c = 0.186) implies slow and low peaked runoff and less erosion.

Horton (1932) defines drainage density (D_d) as the total length of channels per unit watershed area. D_d reflects channel development, according to Langbein (1947) and indicates the physical-environmental conditions that favour erosion in the watershed, as suggested by Rodrigues et al. (2019). A low drainage density watershed sits on surfaces highly resistant to erosion or has a highly permeable substratum (Strahler, 1952), forming a relief with gentle slopes covered by dense vegetation (Zavoianu, 1985). D_d implies runoff characteristics (Fitzpatrick, 2017), whereby a dense stream network watershed will have flashy high peak flows (Horton, 1932). Drainage density seldom exceeds 3.0 and is commonly 1.0 to 2.0 in humid regions where soil erosion is active (Horton, 1932). Results of the analysis show that the D_d of the two watersheds is close, about 2 km per square kilometre of the watershed area, indicating that the runoff characteristics do not differ across the two catchments.

The drainage texture (D_t) represents the total number of streams per unit perimeter of a watershed (Horton, 1932) and reflects the topography’s fundamental lithology, infiltration potential, and relief (Rai et al., 2018). Smith (1950) classified D_t into coarse (≤4), moderate (4–10), and fine (>10). The drainage textures of LW and WIW fall between 5 and 8, described as moderately textured drainage, which, according to Fenta et al. (2017), is associated with the dominance of low permeability soft rock with low erosion resistance.

The stream frequency (F_s) is the total number of streams in a watershed divided by area. F_s may not be comparable for watersheds of different sizes (Horton, 1932). However, the study results show that the F_s of the two watersheds are almost the same at about two streams per square kilometre watershed area.

The compactness coefficient (C_c) represents the ratio of the actual basin perimeter to the perimeter of a circle of equal area (Horton, 1932). C_c equals unity when the watershed shape is a perfect circle, increases to 1.128 in the case of a square, and may exceed 3 for a very elongated basin (Zavoianu, 1985). C_c does not depend on watershed size but is highly slope-dependent (Pareta & Pareta, 2011; Verma et al., 2022). Table 3 shows that C_c of WIW is higher than that of LW, which suggests that WIW is more elongated than LW.

Relief aspects of LW and WIW

Relief analysis is critical for developing a thorough understanding of the spatial organization of landforms. Table 4 shows the result of the relief aspect analysis in the study watersheds. The relief aspects are related to the elevation differences between reference points in a watershed.

Table 4. Relief aspect analysis and the computed values in LW and WIW.

Watershed	Emax	Emin	Emean	R	Rr	Rhp	Rn	HI
LW	732	0	271	732	14	0.004	1.457	0.37
WIW	869	0	195	869	17	0.004	1.653	0.22

E_max – Elevation, maximum [m]; E_min – Elevation, minimum [m]; E_mean – Elevation, mean [m]; R – Relief [m]; R_r – Relief ratio [m/km]; R_hp – Relative relief; R_n – Ruggedness number; HI – Hypsometric integral

LW and WIW drain into the sea on opposite sides of the island (see, Figure 1), and their minimum elevations (E_min) are at zero. The island’s highest point (E_max) is inside WIW, but WIW’s mean elevation (E_mean) is 76 m lower than LW. As shown in Figure 4, the river profiles of the two watersheds revealed that a large portion of the LW lies between 100 and 400 m in elevation. The elevation values are critical in the hypsometric analysis, which examines a catchment’s geological development stages.

Figure 4. Elevation profile of the LW and WIW main rivers from the headwaters to their outlets. The main trunk’s length and peak elevations are 89 km and 732 m, and 105 km and 869 m for LW and WIW, respectively.

A watershed’s total relief (R) is the elevation difference between the highest and lowest points. High R is associated with high erosion susceptibility. R has an essential role in landform development, drainage development, surface and subsurface water flow, permeability, and erosional properties of the terrain (Magesh et al., 2011). Between the two watersheds, the WIW has a higher R than LW, which indicates that WIW has a higher potential for erosion.

The relief ratio (R_r; Schumm, 1956) is the elevation difference per unit of watershed length expressed in m/km. Compared to R, the parameter R_r provides a better representation of the relief aspects of a watershed. One practical application of R_r is estimating sediment loss (Schumm, 1956). A higher relief ratio is associated with increased susceptibility to erosion. WIW is more prone to erosion between the two watersheds because of its higher R_r than LW.

The relative relief (R_hp) is the watershed’s total relief per unit perimeter (Schumm, 1956). Across the two study watersheds, R_hp shows identical values. Similar to other relief aspects, R_hp is also directly affecting erosion susceptibility. With this metric, the two watersheds are on an equivalent level of erosion susceptibility.

The ruggedness number (R_n) is the product of the watershed’s drainage density and relief. R_n represents the intricacy of the landform structure and the average slope steepness along its length. In the study of Guidolini et al. (2020), R_n stands as an environmental parameter that directs potential land use regarding soil capability, where higher R_n values connote higher erosion and environmental degradation. Rodrigues et al. (2019) suggested that areas with high R_n are for forest use because of steep slopes and dense drainage networks, whereas areas with low R_n are appropriate for natural cropland use. Table 4 shows that WIW has a higher R_n value than LW, associating WIW with a higher erosion potential than LW.

Figure 5 shows the study watershed’s hypsometric curve and hypsometric integral (HI). The HI values indicated that only 35% and 28% of the total landmass remained after erosion for LW and WIW, respectively. Both watersheds are considered geologically matured landscapes. However, WIW is at an advanced maturity stage because of a lower HI than LW. A previous study on the island landscape revealed a spatial variation of HI ranging from 0.1 to 0.6, indicating young to mature geologic development stages (Torrefranca et al., 2022a).

Figure 5. Hypsometric curve and hypsometric integral (HI) of the watersheds under study.

Implications for soil and water conservation

LW and WIW watersheds offer a range of environmental services, such as the provision of water flow for several economic activities. However, these watersheds are threatened by severe resource degradation, as soil erosion is a serious environmental problem in many Philippines watersheds (Puno & Puno, 2019) due to changing climatic conditions, tectonic settings, and anthropogenic factors. For example, erosion studies in the upper WIW found that areas cultivated for crop production generated a high soil loss of 42.5 t/ha-yr compared to forested and grassland areas with only 0.2 t/ha-yr and 6.4 t/ha-yr, respectively (Genson, 2006). Further, the study also found that high rainfall intensities resulted in high sediment concentrations in monitored rivers.

Land and water resources in the island landscape are also exposed to tectonic activities owing to geographic location. In 2013, a 7.2 magnitude earthquake hit the island and caused landslides that spoiled agricultural areas, cut transportation networks, and blocked irrigation infrastructures, to name a few. Landslides are commonplace on the island, so resource use must consider the long-term impact of landform evolution. One study on the landscape incorporated landform dynamics in identifying potential small-scale hydropower sites (Torrefranca et al., 2022b). Although the study was an initial table assessment, it incorporated proxies of tectonic activity and showed the spatial distribution and the geohazards associated with each potential hydropower site. The information provided insights to data-limited resource managers and energy planners in targeting possible locations and minimizing field investigation costs and time.

Soil and water resource degradation in the watershed are critical issues. Identifying a watershed’s erosion vulnerability status is a precondition for taking appropriate and relevant actions to tackle the situation (Asfaw & Workineh, 2021). However, the lack of hydrology-related datasets and field observations to proceed with watershed characterization is another constraint to hurdle in the watersheds. As a result, many studies employ morphometric analysis as a rapid approach to characterize hydrologic responses of watersheds using topographic data and GIS (Benzougagh et al., 2022; Negese et al., 2021; Puno & Puno, 2019) to proceed with planning for soil and water conservation actions.

Results of the morphometric analysis indicated that both LW and WIW were susceptible to erosion. Hence, there is a need for a more thorough evaluation and assessment of watersheds to better understand resource degradation. As noted in many studies, it is advantageous to disaggregate a large watershed into smaller and manageable units in identifying and prioritizing areas of urgent concern (Benzougagh et al., 2022; Najia et al., (2018).; Negese et al., 2021; Puno & Puno, 2019).

IV. Conclusion

The morphometric analysis conducted in LW and WIW characterizes the watersheds in terms of their susceptibility to erosion using a digital elevation dataset and GIS. It links the watershed’s linear, areal, and relief aspects to hydrologic responses and surface erosional processes. The study made the following conclusions:

1. The linear aspects indicated that LW and WIW are sixth-order watersheds dominated by lower-order streams. The bifurcation ratio of around 4 for both watersheds suggested highly dissected landscapes with a high potential for flood peaks.

2. LW’s areal aspects, such as form factor, circulatory ratio, and drainage texture, implied that its nearly circular shape is associated with runoff’s short concentration time and quick disposal, which encourages flash flooding and erosion.

3. The greater relief characteristics of WIW suggested steep slopes in the watershed, implying increased soil loss generation compared to LW.

4. Both watersheds are susceptible to erosion.

The study recommends further investigations on the two watersheds to identify erosion hotspot areas and enable targeted conservation programs and initiatives. Future studies should disaggregate the watersheds into smaller and more manageable units to identify and prioritize areas for conservation and protection initiatives.

Acknowledgments

The main author is grateful for the Graduate Scholarship Grant from DOST-ERDT at the School of Engineering at the University of San Carlos – TC, Cebu City, and the Bohol Island State University – Bilar Campus support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by a Graduate Scholarship Grant from the Philippines Department of Science and Technology - Engineering Research and Development for Technology (DOST-ERDT) at the University of San Carlos - TC, Cebu City, Philippines.