MANILA BAY WATERSHED SCORECARD: A GIS-BASED QUANTITATIVE WATERSHED HEALTH ASSESSMENT

: Watershed health refers to the maintenance of the normal status of a watershed. With the use of Geographic Information Systems (GIS), a watershed health assessment framework was developed and applied to the Manila Bay Watershed (MBW) based on the sub-indices used by the US EPA and the Minnesota DNR. These three sub-indices were used and modified: geomorphology, connectivity, and hydrology. Geomorphology sub-index accounted for the effects of soil erosion using the Unit Stream Power-based Erosion/Deposition model (USPED). Connectivity sub-index considered the connection between habitats within three environments: terrestrial, aquatic, and riparian zone. Hydrology sub-index accounted for the effects of impervious cover and urbanization on the movement of water using three factors, namely, natural cover, tree cover, and loss of hydrologic storage. Land cover map was the most used dataset in this scorecard, where the maintained natural land covers generally received high scores and built-up areas received the lowest. The overall watershed health score of MBW is 75.713 from the mean of the three sub-indices. Pampanga River Basin, which is the largest river basin within the MBW, got the highest score of 79.462 since it consists of huge portions of maintained natural land cover. Manila River Basin, known to have dense built-up areas, got the lowest average of 60.773. On the provincial level, the province of Nueva Ecija got the highest score, and the National Capital Region (NCR) got the lowest. The developed framework successfully quantified a relative health score which can be used to rank and prioritize subwatersheds, and to measure in totality the improvement or degradation of subwatershed


INTRODUCTION
Watershed health refers to the maintenance of the "normal" status of a watershed's complex adaptive system, and a healthy watershed can provide services essential to humans ranging from water supply to cultural benefits and ecological functions (Mosaffaie, et al., 2021).Although the concept of quantifying the health of an ecosystem has been realized by researchers since the 1940's, the process of quantifying the ecosystem health still does not follow any direct methods (Banarjee, et al., 2017).In this research, the Manila Bay Watershed (MBW) Scorecard was derived from the frameworks developed by the US Environmental Protection Agency (US EPA) and the Minnesota Department of Natural Resources (Minnesota DNR) which commonly used certain sub-indices to calculate for the score.The three sub-indices used in this scorecard are as follows: Geomorphology, Connectivity, and Hydrology, where each index corresponds to a specific category for watershed health assessment.The Watershed Health Score is the mean of means of the three sub-indices.

Watershed Scorecard
The US EPA developed the Healthy Watersheds Program based on six sets of indices: Landscape Condition, Habitat, Geomorphology, Hydrology, Biological Condition, and Water Quality.These are sets of ecological attributes that are measurable, comparable, and consistent across the area of assessment.(US EPA, n.d.)On the other hand, the Minnesota DNR developed a similar framework called Watershed Health Assessment Framework (WHAF) defined as a structured and science-based approach to grow a common understanding of Minnesota's complex natural resource systems (Minnesota DNR, n.d.a).WHAF uses health scores organized into five ecological components, namely, Hydrology, Geomorphology, Biology, Connectivity, and Water Quality, to compare the health of ecological systems in the state of Minnesota (Minnesota DNR, n.d.b).Considering the data availability for the assessment of the MBW, the developed scorecard reflects the same ecological attributes used in the programs developed by the US EPA and the Minnesota DNR.These attributes are organized into three subindices: Geomorphology, Connectivity, and Hydrology.

Geomorphology
Geomorphology is the study of landforms, which are produced when rocks and sediments are eroded, transported, and deposited from one area to another.(British Society of Geomorphology, n.d.)These landforms change due to natural and anthropogenic processes shaping the land which rivers eventually flow.The shape of each stream reacts to certain variables in predictable and measurable ways.Healthy streams are said to be able to carry a certain amount of sediment over time in a sustainable balance.The relationship between discharge and sediment transport is important because it determines whether the stream channel is stable, aggrading, or degrading.(Minnesota DNR, n.d.c)

Unit Stream Power-based Erosion/Deposition (USPED) model
USPED is a simple model designed for complex terrain, soil, and cover conditions which predicts the spatial distribution of erosion and deposition rates for a steady state overland flow with uniform rainfall excess conditions (Mitas & Mitasova, 1998).It is a 2dimensional soil erosion model which assumes that soil erosion and deposition depend mainly on the sediment transport capacity of the surface runoff.This assumes that when soil particles are detached by rain, but if there is not enough runoff to transport the soil particles due to the terrain shape or land cover, the amount of erosion will be significantly reduced.(Liu, et al., 2007)

Connectivity
Connectivity is defined as the extent an organism moves through a landscape to access resources for its survival (Merriam, 1984).Two important elements are considered when talking about connectivitythe physical structure of the landscape and the biology and behavior of the organisms present (Taylor, 2000).
Landscape connectivity is the extent to which movement is facilitated or impeded through different patch types across the landscape (Taylor, et al., 1993).Fragmentation of landscapes changes the movement of animals, plants, and water.Connected landscapes allow the movement of organisms to reproduce and maintain genetic diversity (Minnesota DNR, n.d.d).On the other hand, water resources and transportation infrastructures such as dams, bridges, and culverts provide countless benefits including water security, transportation, power, and flood control.However, they also induce certain ecological concerns such as decreasing the connectivity within the river ecosystem for the movement of organisms, sediments, and water itself (McKay, et al., 2016).
Between the aquatic and terrestrial ecosystems lies the riparian zones which act like boundaries affecting the food web dynamics and ecosystem function of both habitats.Riparian connectivity is directly affected by human activities.Urbanization, for example, directly affects the terrestrial communities and in turn alters the connectivity between terrestrial and its adjacent aquatic ecosystems (Kupilas, et al., 2021).

Hydrology
Urbanization has caused changes in watershed hydrology resulting to the alteration of watershed sediment and solute export.Construction of dams and impoundments, for example, caused the decline of the natural filtering capacity of river systems and the regulation of water flow.Urbanization has also increased the spread of the total impervious areas (TIA) in the form of parking lots, roads, lawns, and rooftops.This reduces infiltration and surface storage of rainfall which causes surface runoff to increase.Aside from TIA, vegetation also affects the stream's hydrological response to precipitation.Through simulations, it was found that increased tree cover could reduce peak flows by up to 12%, and trees in urban catchments could intercept up to 41% of the precipitation during storms.(O'Driscoll, et al., 2010) The presence and densification of canals and drainage systems also influence water storage and water quality.Direct connection between the stormwater drainage systems and surface waters, and bypassing the riparian vegetation, reduces the ability of the riparian zone to decrease nitrogen concentration as water enters the stream network.(O'Driscoll, et al., 2010)

METHODOLOGY
The MBW was divided into subwatersheds using built-in tools in ArcGIS.Pour points were manually chosen considering upslope contributing areas.Delineated subwatersheds were then split by municipality to assess the health of the subwatershed within the jurisdiction of each.These subwatersheds were then evaluated using the three sub-indices.All health analyses were done per subwatershed split by municipality.

Geomorphology Sub-index
The Geomorphology Sub-index is the normalized inverse of the Sediment Transport Capacity (T) which was computed using the Unit Stream Power-based Erosion/Deposition (USPED) model with the following five factors: terrain factor, soil erodibility factor (K), cover management factor (C), support practice factor (P), and rainfall-runoff erosivity factor (R). Fig. 2 shows the workflow for calculating this sub-index.

Sediment Transport Capacity using USPED
Based on the USPED model, T of a given land pixel can be determined from the total water flow through this given land pixel.The equation for the sediment transport capacity is: where R = rainfall-runoff erosivity factor K = soil erodibility factor C = cover management factor P = support practice factor b = slope A = upslope contributing area R-factor approximates the uniform rainfall intensity.On the other hand, factors K, C, and P approximate the transportability coefficient Kt.The expression   • (sin )  approximates the terrain factor, which is described in the next section.

Terrain Factor (LS)
=   • (sin ) (2) The Terrain Factor was calculated using equation ( 2) where A is the upslope contributing factor derived using the Flow Accumulation tool in ArcGIS.b, on the other hand, is the slope derived from the DTM using the slope tool also in ArcGIS.m and n are empirical exponents such that m = 1.6 and n = 1.3 for prevailing rill erosion, while m = n = 1 for prevailing sheet erosion.(Liu, et al., 2007;Blanco & Nadaoka, 2006).In this research, m = n = 1 was used.

Soil Erodibility (K-Factor)
Soil Erodibility Factor (K-Factor) is a quantitative description of the susceptibility to erosion of a particular soil.Every soil has its own rate of erosion based on its physical characteristics such as its organic matter content and soil texture classification (Wajid, et al., 2020).The soil map was acquired from the FAO Digital Soil Map of the World (DSMW).The dominant topsoil texture for each classification in the FAO Soil Map was acquired from the Harmonized World Soil Database (HWSD).Based on the dominant topsoil texture present in the MBW, the soil map raster was reclassified with the K-factor values obtained from Stewart, et al. (1975).

Crop Management (C-Factor) and Support Practice (P-Factor)
C-Factor represents the ratio of soil loss from land cropped under specific conditions to the corresponding loss from a tilled, continuous fallow condition (Teng, et al., 2016).For the MBW, C-factor was derived by reclassifying the land cover map using the values obtained from Yang, et al. (2003) and Bakker, et al. (2008).On the other hand, P-factor is the ratio between the soil erosion with a specific support practice and the corresponding soil loss with straight-row upslope and downslope tillage.The P-factor also accounts for the control practices that reduce the erosion potential of the runoff by their influence on drainage patterns, runoff concentration, runoff velocity and hydraulic forces exerted by runoff on the soil (Renard, et al., 1997).For the MBW, P-Factor was derived by reclassifying the land cover map using the values obtained in Yang, et al. (2003).R-Factor is the sum of individual storm EI-values for a year averaged over a long period of time.However, other research studies provide a different approach to obtain the R-Factor.
In this research, the average monthly and annual precipitation values were derived from daily precipitation data from CHIRPS v2.0 package (de Sousa, et.al, 2020).

Connectivity Sub-index
The Connectivity Sub-index was computed from the mean of three components: Terrestrial, Aquatic, and Riparian Connectivity.Each component provides different parameters concerning the condition of connectivity for each habitat.The workflow is shown in Fig. 3.

Terrestrial Connectivity
Terrestrial Connectivity was calculated based on the percentage of like adjacencies of patches of the same land cover.This index was computed using Fragstats, a spatial pattern analysis program used to quantify landscape structure.Percentage of Like Adjacencies (PLADJ), one of the landscape aggregation metrics available on Fragstats, uses this equation: With the land cover map as input raster, Fragstats initially identified the patches of land cover types and created an adjacency matrix.Then, PLADJ was calculated by computing for the sum of like adjacencies between patches of the same land cover type and divided it by the total number of adjacencies.PLADJ is higher for landscapes with greater aggregation of patch types like those with larger patches and compact shapes as compared to landscapes with disaggregated patch types or those with smaller patches and complex shapes.(McGarigal, et al., 2012)

Aquatic Connectivity
Aquatic Connectivity was computed from the mean of three factors: Road-Stream Crossing Density, Canals and Drains Density, and Stream Patch Size.These factors are significant in assessing the connectivity of streams within the subwatersheds.Thus, subwatersheds without streams and those without available stream data are not included in this analysis.
Road-Stream Crossing Density was derived from the number of road-stream points or the overlap between the road and stream vector polylines.The density was computed by getting the ratio between the number of road-stream points within a sub-watershed and the area of the sub-watershed they are in.The score is the inverse of the density which gave higher points for areas with lower road-stream overlap per area ratio.These areas are considered to have the highest connectivity considering the effects of bridges and culverts.Canals and Drains Density was computed by getting the ratio of the total length of canals and drains divided by the total flow line length of streams within the sub-watershed.The score is also the inverse of the density, giving higher points for areas with lower canals and drains length per total flow line length ratio.Stream Patch Size considers the discontinuity of streams caused by dams and reservoirs.With that, portions of the stream network line vectors that overlaps with the reservoir and dam polygons were erased.Ideally, Stream Patch Size is computed considering the total length of connected streams within the subwatershed.However, since the subwatersheds in this research were also split by municipal boundaries, the whole stream network may not be within the area of interest.With this, the length of stream considered was only the total length of streams within the area of interest but are connected to each other at some point upstream or downstream.

𝑆𝑡𝑟𝑒𝑎𝑚 𝑃𝑎𝑡𝑐ℎ 𝑆𝑖𝑧𝑒
Considering the equation above, Stream Patch Size is higher for areas where streams are connected, and smaller for disconnected streams.This means that scores are lower if the disconnected streams passing through the area are almost equal in length.

Riparian Connectivity
Riparian Connectivity was computed using the land cover map and the stream polyline.Riparian areas were identified by creating a 200m buffer zone from the streams.The land cover map was also reclassified to differentiate natural and impervious land covers.Scores were calculated based on the percentage of the natural land cover of the riparian area within each subwatershed.Since Riparian Connectivity also involves stream data, subwatersheds without streams or those without stream data are not included in this analysis.

Hydrology Sub-index
Hydrology Sub-index was computed from the mean of three components -Natural Cover, Tree Cover, and Loss of Hydrologic Storage (Fig. 4).All of which provides certain emphasis on hydrological conditions affecting the health of the watershed.Natural Cover was computed from the land cover map.Areas identified to have croplands, wetlands, aquaculture, forestland, mangroves, paddy rice, shrubland, and water were classified as natural land covers while areas with urban and barren land as land covers were classified as impervious.Scores were computed from the percentage of classified natural land covers within each subwatershed.Meanwhile, existing data was used for the Tree Cover and scores were proportional to the percentage of tree cover per subwatershed.On the other hand, scores for the Loss of Hydrologic Storage are just the same as the scores acquired for the Canals and Drains Density of the Aquatic Connectivity Index.

Geomorphology Sub-index
The Sediment Transport Capacity was computed using the USPED model.High sediment transport capacity means that a higher sediment load can be transported at a given flow rate therefore sediments are more susceptible to erosion.On the other hand, the Geomorphology Sub-Index is the normalized inverse of the Sediment Transport Capacity such that an area with relatively high erosivity will have a low Geomorphology score.The mean score is 87.096 with a standard deviation of 18.183.Most of the subwatersheds have scores between 90-100.Areas with a score of 100 are mostly those with slope equal to zero since terrain factor alone yields a score of zero for Sediment Transport Capacity.On the other hand, areas that have relatively low scores are mostly where rivers are located.Fig. 5 shows that significantly low scores can be seen along Pampanga River.

Terrestrial Connectivity
The mean score for Terrestrial Connectivity is 81.359 with a standard deviation of 18.484.Score distribution for Terrestrial Connectivity is highest between 80-90.High scores were given in areas where dominant land cover types can be easily identified.Most of which are in the mountainous areas of the watershed where land cover is dominated by forests and grasslands.High scores were also observed in the flat lands where paddy rice is the dominant land cover type while relatively low scores are mostly observed in areas where urban and annual crop lands are scattered in smaller patches.

Aquatic Connectivity
Aquatic Connectivity mean value is 73.4 with a standard deviation of 23.657.The score distribution is highest between 60-70.However, subwatersheds with scores between 90-100 is also noticeably high.Subwatersheds with a score of 100 are those considered to have streams with no discontinuity as it passes through these subwatersheds.Most of which are in the mountainous areas in the eastern and western-most portions of the watershed.On the other hand, areas with high stream discontinuity or those with relatively low scores can be observed in the areas of Metro Manila and nearby areas of Cavite, along the Pampanga River, and near the Pantabangan reservoir.

Riparian Connectivity
Riparian Connectivity mean score of 66.59 with a standard deviation of 32.269 was obtained for the MBW.Most of the scores are within 90-100.Since Riparian Connectivity considers the effect of manmade structures, areas near Metro Manila obtained the lowest scores.High scores for Riparian Connectivity were obtained for areas in the mountainous regions mostly on the eastern part of the watershed.

Connectivity Sub-index
Connectivity Score was calculated from the mean of the three connectivity components previously discussed.A high score is an indication that habitats within that certain subwatershed are connected with each other, thus making them healthier than subwatersheds with lower scores.Most subwatersheds obtained relatively high scores with a mean score of 73.269 with a standard deviation of 19.964 (see Fig. 6).However, it is significant to note that some areas may not have all three connectivity components since not all subwatersheds have streams and data may not be available.

Natural Cover
Metro Manila is a well-known highly urbanized region in the Philippines.Thus, the percentage of natural cover left is at the lowest compared to nearby regions within the watershed.Other regions outside Metro Manila obtained a relatively high score and majority of which have scores between 90-100 producing a mean of 80.89 with a standard deviation of 29.195.

Tree Cover
The scores for the Tree Cover obtained the lowest mean score of 43.197 with a standard deviation of 22.593.Majority of the subwatersheds have scores falling under the range of 30-40.It is also important to note that the lowest values are in areas classified as aquaculture in the land cover map since trees are not present at all in these areas of the watershed.These can mostly be observed on the areas surrounding Manila Bay.

Loss of Hydrologic Storage
The data used for the Loss of Hydrologic Storage component is the same with the data used for the Canals and Drains Density.The mean score for the whole subwatershed is 85.051 with a standard deviation of 35.599.High value for the standard deviation is due to the significant imbalance in the distribution of the scores where most of the scores fall under the range of 0-10 and 90-100 only.This means that most of the subwatersheds either have a lot of manmade canals and drains or none at all.

Hydrology Sub-index
Hydrology Score was computed from the mean of the three components where the mean score is 67.066 with a standard deviation of 21.730.As seen in Fig. 7, most of the low values can be seen in the Metro Manila area.Provided that most of the components used for the calculation of this sub-index deals with the prominence of the natural cover, this also means that the hydrological condition of the watershed in highly urbanized areas are relatively unhealthy.

Watershed Health Score
Landcover map was the most used dataset in this scorecard.From the results of the indices, natural land covers such as croplands and forests lands generally received high scores and built-up areas received the lowest.Most of the watershed health scores fall between the range of 80-90.High scores indicate a healthier watershed condition, while low scores indicate a relatively unhealthier condition.Subwatersheds with a score of 100 are those in the healthiest relative condition, while those with zero are the unhealthiest.The mean score for the whole MBW is 75.713 with a standard deviation of 14.551.(See Fig. 8) Pampanga River Basin is the largest river basin within the MBW.It has the highest computed watershed health score of 79.462.Although   From Table 2, the highest mean score per province was calculated for the province of Nueva Ecija.Additionally, high scores were also calculated for the mountainous areas in the western portion of Bulacan and in the eastern-most portion of Bataan and Pampanga (see Fig. 9).Meanwhile, the lowest scores were obtained in the NCR Districts.The lowest mean score was calculated for the first district of NCR, where the capital of the Philippines, the City of Manila, is located.Not considering the rest of the NCR Districts, the province of Cavite came in second to the last and the province of Rizal follows.Relatively low scores can be observed for the subwatersheds along the Pampanga River in the provinces of Bulacan, Nueva Ecija, and Pampanga.Similarly, low scores were calculated for some areas surrounding the Laguna Lake in the provinces of Rizal and Laguna.
On the municipal level, the municipality of Dupax del Sur got the highest score of 91.163 followed by San Luis (89.333), in the province of Aurora, and Alfonso Castañeda (89.035) in the province of Nueva Vizcaya.On the other hand, the Port Area of the City of Manila got the lowest score of 29.791.Followed by Cavite City (31.917) in Cavite, and Tondo in the City of Manila (36.226).

CONCLUSIONS AND RECOMMENDATIONS
The methodology is a framework created to provide a means to quantify watershed health.Through the various sub-indices and components, the spatial variation of watershed health across the MBW was elucidated.The results were only able to quantify a relative health score which can be used to rank and prioritize subwatersheds.This score cannot quantify the health of a subwatershed alone.However, the score card can be used to somehow measure in totality the improvement or degradation of subwatershed/s over time because of anthropogenic and natural factors.
Although it successfully provided quantifiable data for watershed health assessment, certain aspects must still be considered to come up with a more specific and refined scoring system.However, the most limiting aspect of this research is data availability.Significant data may not be available or not all locations might have the data needed for the refinement of the scoring system.If data is openly available, using additional indices and sub-indices could help refine the scoring system.Water quality, for example, considering the nitrogen and phosphorus content, as well as turbidity and pollution sources may be used.This will provide a scoring system more specific to the quality of water which is a significant aspect in watershed health assessment.Another would be the biological aspect considering the abundance of species, and the introduction of endemic species through artificial links which affects the health of a watershed.

Figure 1 .
Figure 1.Map of Manila Bay Watershed and all river basins.
monthly precipitation P = average annual precipitationUsing the Modified Fournier's Index in equation (3),Arnoldus (1977) obtained a relation between the mean annual rainfall and the mean annual R-value: −  = 0.2641.50

Figure 2 .
Figure 2. Workflow for the Geomorphology Sub-index.

Figure 3 .
Figure 3. Workflow for the Connectivity Sub-index.

Figure 4 .
Figure 4. Workflow for the Hydrology Sub-index.

Figure 5 .
Figure 5. Map of the Geomorphology Sub-Index.

Figure 6 .
Figure 6.Map of the Connectivity Sub-index.

Figure 7 .
Figure 7. Map of the Hydrology Sub-index.

Figure 8 .
Figure 8. Map of the Watershed Health Score of the subwatersheds within the Manila Bay Watershed.itcontains cities with high densities of built-up areas, a huge portion of the basin has a natural land cover including the mountainous areas on the west and the agricultural areas in the central plane regions.This thus gave the basin a relatively high score indicating that the Pampanga River Basin is the healthiest among all river basins within the MBW.On the other hand, is the Manila River Basin.As the location of the country's capital, it is known to have a large density of built-up areas where natural covered lands have already been developed as urban hubs.This thus resulted to the lowest average of 60.773 indicating that the Manila River Basin is the unhealthiest within the MBW.

Figure 9 .
Watershed health per province of subwatersheds within the Manila Bay Watershed.

Table 1 .
Watershed health score of the major river basin in the Manila Bay Watershed.

Table 2 .
Mean watershed health score of provinces within the Manila Bay Watershed.