Determining Conservation Priorities in the Urban Citarum Watershed, West Java: An Ecosystem Services Approach

2.1 Study area

The study was conducted in the Citarum Watershed, Indonesia (Figure 1).

1.png

Figure 1. Research location of the Citarum Watershed, West Java Province, Indonesia

It is located between latitudes 106°51′ - 107°51′E and longitudes 7°19′ - 6°24′S, which spreads throughout eight regions in West Java Province and covers a total area of 690,916 hectares. The area has a tri-monthly dry climate with 2,358 mm of annual rainfall on average. Three large dams (Saguling, Cirata, and Jatiluhur) bridge the Citarum River, which flows through the watershed and is an essential source of fresh water, agriculture, and electricity for many West Javanese [25, 26].

Three landform-based areas define the Citarum Watershed: upstream, middle, and downstream. The Bandung Basin, which is the upstream portion, is situated at an elevation of 625 and 2,600 meters above sea level. The primary components of its geological composition are tuff, lava, lapilli, and breccia. The average annual rainfall in the highland and mountainous regions of the upper region is 4,000 mm, with a minimum temperature of 15.3℃. The upper catchment contains various soil types, including latosol (35.7%), andosol (30.76%), alluvial (24.75%), red-yellow podzolic (7.72%), and regosol (0.86%) [27].

2.2 Data sources

This study used primary data from Landsat-8 OLI imagery (2020) and secondary data from previous research on WY [19], SC [20], and CS [21]. Software like Google Earth Engine (GEE), the InVEST model, and ArcGIS 10.8 were employed for data processing [9, 28]. The InVEST model was utilized as a spatial analysis tool to evaluate ES, including annual WY, SC, and CS. All data were processed using the WGS84 datum with a 30-meter spatial resolution. To date, more consideration is placed on the significant role of ES in decision-making processes for sustainable natural resource management. Therefore, regional development planning should account for ES dynamics resulting from changes in LULC [29].

Changes in LULC can lead to land degradation, hinder the provision of ES in a specific area, and affect the efforts to develop sustainable ecosystems. These changes include modifications to land use types and adjustments in intensity and spatial patterns [1].

The use of 2020 data was part of the evaluation of the West Java Province Long-Term Development Plan (PJP) from 2000 to 2020 [30]. This analysis aims to identify trends over the 20-year period and project changes for the next 20 years, including potential risks, impacts, and trends related to anticipated environmental dynamics [31].

2.3 Methods

The research was executed in four phases, namely (1) analyzing Landsat-8 OLI images; (2) mapping ES (WY, SC, and CS); (3) mapping ES hot spots; and (4) identifying priority areas for conservation. The research workflow is illustrated in Figure 2.

2.png

Figure 2. Research framework

2.3.1 Deriving LULC maps

The GEE platform was utilized to generate cloud-free satellite images. The Landsat-8 OLI datasets included atmospherically corrected surface reflectance that was processed using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [28]. Supervised classification with a Random Forest (RF) algorithm on the GEE platform was used to categorize Landsat pixels into LULC 11 classes.

2.3.2 Mapping various ES

The InVEST tool was employed to measure and visualize ES [9]. The circumstances of the region and the accessibility of geographic and spatial data led to the selection of three particular ES: WY, SC, and CS.

WY. Based on the idea of water balance, this module determines each pixel's yearly WY ($W Y_x$) in landscape x by dividing the average annual precipitation by the actual annual evapotranspiration [9] using Eq. (1).

$W Y_x=\left(1-A E T_x / P_x\right) \times P_x$                  (1)

where, $P_x$ is the yearly precipitation on pixel x and $A E T_x$ is the yearly actual evapotranspiration for pixel x. The InVEST user manual contains a detailed computation procedure [9]. The data required for the WY module contains the volumetric plant accessible water content, root limiting layer depth, average annual potential evapotranspiration, plant evapotranspiration coefficient, and LULC map.

CS. The "pool" sizes of four CS—soil organic matter, dead organic matter, belowground biomass, and aboveground biomass—have a major impact on how much carbon is stored in ecosystems. Quantification of CS was done using the carbon module of the InVEST model. The CS in these pools is combined based on the LULC maps in this model. To quantify CS, this study derived data from previous research [32, 33]. Thus, the total CS_x for each pixel on the LULC type x was calculated using Eq. (2).

$C S_x=C_{x a}+C_{x b}+C_{x s}+C_{x d}$                 (2)

where, the carbon concentrations in soil (Mg C ha^-1), dead matter (Mg C ha^-1), aboveground biomass (Mg C ha^-1), and belowground biomass (Mg C ha^-1) were represented by C_xa, C_xb, C_xs, and C_xd, respectively, for each pixel with LULC type x.

SC. This study evaluated SC using the "sediment delivery ratio" module within the InVEST model. It calculates SC by measuring the difference between potential and actual soil loss before applying the Revised Universal Soil Loss Equation (RUSLE) to estimate soil loss per pixel [9]. The SC for a given pixel x was calculated using Eqs. (3) to (5).

$S C_x=R K L S_x-U S L E_x$                  (3)

$R K L S_x=R_x \cdot K_x \cdot L S_x$                 (4)

$U S L E_x=R_x \cdot K_x \cdot L S_x \cdot C_x \cdot P_x$                   (5)

where, $U S L E_x$ is the actual soil loss for pixel x (ton·ha^-1·year^-1), $R_x$ is the rainfall erosion factor for pixel x (MJ·mm·ha^-1·h^-1·year^-1), $K_x$ is the soil erosion factor for pixel x (ton·ha·h·ha-1·MJ-1·mm-1), and $S C_x$ is the SC (ton·ha^-1·year^-1).

$L S_x$ is the slope length and gradient factor (dimensionless) while $C_x$ and $P_x$ represent the vegetation cover and support practice factors (dimensionless, ranging from 0 to 1), respectively. Both R_x and K_x can be estimated by methods defined in the literature [34]. This study assigned C and P values according to existing literature [9, 33, 35].

2.3.3 Mapping ES hot spots

The Getis-Ord statistics (Gi*) were used in this study to determine the hot spots and cold spots for ES. This technique is available in ArcGIS 10.6 and assesses every raster pixel through its surrounding pixels [36, 37]. Such analysis produces a new feature class with matching p-values, z-scores, and confidence levels. Low p-values and high z-scores indicate statistically significant hot spots, whereas low p-values and high negative z-scores suggest statistically significant cold spots.

A spatial analytical method called hot spot analysis explores the features and patterns of high and low ES values. It evaluates statistical significance through the standardization of z-values. Eqs. (6) to (8) were applied for calculating county unit I's (Gi).

$G_i=\frac{\sum_{j=1}^n w_{i j} x_j}{\sum_{j=1}^n x_j}$               (6)

$\left(G_i^*\right)=\frac{\sum_{j=1}^n w_{i j} x_{j-} \bar{x} \sum_{j=1}^n w_{i j}}{S \sqrt{\frac{\left[n \sum_{j=1}^n w_{i j}^2-\left(\sum_{j=1}^n w_{i j}\right)\right]}{n-1}}}$              (7)

$S=\sqrt{\frac{\sum_{j=1}^n x_j^2}{n-1}}-(\bar{x})^2$                (8)

where, Gi^* is the county unit I's spatial agglomeration index, W_ij is the weight matrix between county units I and J, x_i  and x_j are the attribute values for the i-th and j-th county units, n is the total number of counties in the research area, x is the average attribute value across all counties, and S is the standard deviation of attribute values across all counties.

2.3.4 Identify priority areas for conservation

The study area was divided into a total of 7,239 planning units, each measuring 1.000 m × 1.000 m. Following the Indonesian National Standard SNI 7645-1:2014 on the classification of medium land cover and [38], medium scale is defined as the classification of land cover at scales of 1:250,000, 1:50,000, and/or 1:25,000. This study falls within the scale of 1:100,000, thus complying with SNI 7645-1:2014.

Based on planning unit (PU), the average ES values were divided into a 1,000 m × 1,000 m grid using ArcGIS, thus creating the planning units. Zonal statistics were calculated for the map of the three ES within the planning units, allowing each grid to obtain values for the three services. To ensure consistent representation, the planning unit values were standardized using Eq. (9).

$S_i=\frac{E_i-\min \left(E_i\right)}{\max \left(E_i\right)-\min \left(E_i\right)}$                   (9)

where, E_i is the i-th ES value, $\max \left(E_i\right)$ and $\min \left(E_i\right)$ are the i-th ES maximum and minimum values, and $S_i$ is the standardized value.

The hot spot and cold spot maps were reclassified using the Getis analysis. All hot spots and cold spots with a 90% confidence level were deemed insignificant and only those with 95% and 99% confidence levels were highlighted [39]. A composite map was generated by overlaying individual ES hot spot map. On each ES map, pixels designated as hot spots and cold spots were given a value of 1 and -1, respectively. Meanwhile, a value of 0 was assigned to every other pixel. After the three layers were integrated, a map representing the number of overlapping hot spots and cold spots was created, with pixel values ranging from -3 to 3 [40].

On the aggregate hot spot map, grids with higher values indicate higher priority for multiple ES protections to determine priority areas. The number of ES hot spots was used to designate priority regions for conservation [41]. In this study, all ESs were assumed to provide equal benefits, so each ES was given similar weight.

In this study, conservation area planning was conducted on a regional scale using the Getis analysis, which provides valuable insights for initial conservation area assessments. Employing Systematic Conservation Planning (SCP) methodologies like Marxan is essential for more detailed and effective regional-scale conservation planning when financial resources are limited.

The Getis analysis method is well-regarded for effectively detecting spatial clusters and hot spots, making it useful for informing conservation planning. It has been applied in several studies to evaluate spatial patterns and support conservation-related decision-making [42].

SCP frameworks utilize algorithms to evaluate multiple spatial solutions for selecting conservation areas, ensuring that critical biodiversity components are represented [43]. New frameworks, such as those developed for municipal ecological conservation in China, incorporate multiple indicators and decision-making scenarios, thus balancing conservation benefits with economic feasibility. The InVEST model, alongside Marxan, exemplifies the integration of ecological and socioeconomic factors in prioritization efforts [44].

By integrating cost considerations, SCP approaches help planners allocate resources efficiently, ensuring that conservation investments target the most effective areas for protecting biodiversity [45]. It identifies regions with high positive z-scores (greater than 0) as hot spots, thus representing significant ecological importance or vulnerability. In contrast, regions with low negative z-scores (less than 0) are cold spots, reflecting low ecological activity or concern. To assess the impact of specific variables, they are excluded one at a time from the analysis before recalculating the Gi* statistics. Comparing these results will provide valuable information on the variables that significantly influence the hot spot and cold spot distributions. This process helps pinpoint the key factors affecting the conservation objectives [46].

As economic activities and human populations grow, effective LULC management through the ES approach becomes vital for balancing conservation and economic interests [2]. Changes in LULC significantly impact ES, including energy exchange, soil erosion, and water cycles, leading to both direct and indirect effects on ecosystem functionality [47, 48]. Forests and grazing lands typically provide a higher supply of ES, while agricultural and developed areas often exhibit deficiencies [49].

Such models are especially valuable in conflict-prone regions, where economic pressures threaten the delicate balance of natural systems. Shifts in LULC, especially conversion of vegetated to developed or agricultural land, disrupt regional ecological functions by reducing biodiversity, carbon storage, and groundwater supply while increasing soil erosion and habitat fragmentation [50].

LULC management serves as the foundation for the sustainable utilization of ES. Without it, the delicate balance of our environment will become increasingly vulnerable [51]. Spatial analysis is evolving beyond a mere tool; in this context, it is essential for effective decision-making. By allowing for sophisticated, strategic choices informed by the present and mindful of the future, spatial analysis empowers us to unravel the complex web of ecological relationships [52].

In the study area, hot spots are primarily distributed in forest and grassland areas, whereas agricultural and built-up areas are dominated by cold spots (see Table 5). The loss of forest—as seen in cases like the Marghazar Valley where forest cover decreased by 33% while urban areas expanded by 38.4%—clearly demonstrates the pressures of urbanization on natural resources [53].

Although the ES framework offers a path to balance economic and conservation priorities, it is still a challenge to align immediate economic benefits with long-term ecological health. This balance is essential to identify areas at risk from LULC changes for water regulation, biodiversity, and other services [49]. Understanding ES distribution is fundamental for effective LULC management due to its significant implications for human well-being and ecosystem functionality [47, 48].

In conservation planning, areas with greater spatial overlap are more likely to succeed due to the strong spatial connectivity of ES hot spots [40]. Ecosystems that feature interconnected ES areas tend to maintain higher resilience and support stronger service networks [54]. Integrating spatial dynamics improves conservation efforts by targeting key services, such as WY and SC [55]. Although LULC changes typically reduce ES, strategic land-use approaches can enhance resilience and maintain service provision by reinforcing ecosystem adaptability.

In areas with natural vegetation and relatively limited economic activities, hot spots were encountered more frequently (with frequencies of 3 and 2 times) compared to areas with high activity, such as residential zones and shrublands, where hot spots were found only once or were even dominated by cold spots [56].

The interaction between natural vegetation, economic activities, and climate factors plays a crucial role in shaping the distribution of ecological hot spots and cold spots across different landscapes. Areas dominated by natural vegetation with minimal economic activity showed more frequent hot spot occurrences, recorded at three and two instances, respectively. Conversely, regions with intensive activities, such as residential areas and shrublands, exhibited fewer hot spots, often appearing only once or being dominated by cold spots.

This pattern underscores how land use and economic activity influence the spatial distribution of ecological values. Research suggests that urban green spaces, such as Central Park in Helsinki, function as intricate social-ecological systems. These spaces offer a range of ES that can be spatially mapped to identify hot spots and cold spots related to landscape values and visitor activity. Previous findings revealed a low overlap between landscape value hot spots and visitor use, indicating that ecological quality does not always coincide with patterns of human activity [57].

The variability between ES hot spots and cold spots depends largely on the selected ES type. Notably, regions of high ES provision may lack conservation focus, especially if they face intense human activity that degrades their support for regulating services. Conservation funding limitations necessitate prioritization of hot spots with low edge-to-area ratios, enabling focused resource application [12].

Urbanization, agriculture, and agricultural development processes gradually move away from critical ES locations, gradually reducing their ability to provide essential habitats and retain water. According to previous study [58], ecological harmony is disrupted when these critical parts are damaged. Past study reported that prioritizing hot spots with minimal edges increases ecosystem connectivity [59].

Human actions have a significant effect on ES and often diminish their ecological value as ecosystems shift from multifunctional to single-service dominance, which makes conservation efforts more difficult [60]. To address this, multi ES planning aligns conservation strategies with available resources while enhancing ecosystem resilience [61].

The Getis-Ord method is one of the powerful tools that identify high and low ES clusters to explain the spatial patterns of ES distribution [62]. In the Sichuan Basin, spatial modeling has provided important insights by mapping carbon sequestration and soil conservation patterns, conversely revealing how local drivers influence ES distribution and availability [63]. Moreover, monitoring temporal shifts in ES is no longer a luxury; it is a necessity. Tracking these changes can provide a holistic approach to resource management, which is critical for securing water and food resources in a changing world [64].

Hot spots are vital for ES like soil conservation, which can be spatially mapped to identify areas of high service provision. For instance, in Shaanxi Province, hot spots provided 59.7% of total soil conservation services despite occupying only 29.6% of the area. Getis-Ord Gi* statistics can assist ES conservation efforts by identifying multi-functional priority areas for conserving multiple ES and biodiversity [65].

Hot spots are not only abundant in biodiversity but also offer vital ES like carbon sequestration, water purification, and soil fertility. While the spatial overlap between regions rich in biodiversity and those delivering essential ES is typically minimal, focused conservation initiatives can safeguard both [66]. Meanwhile, hot spots are where high biodiversity, ES, and threat overlap cover 0.1 to 7.1% of the area and targeting them will support reactive conservation strategies [67].

Ecologically rich areas, where biodiversity and resilience converge in a delicate balance, are often labeled as overlapping hot spots. These areas are vibrant with life and serve as the important, irreplaceable pillars for the health of ES rests. It demands more resources, strategic management, and attention. Therefore, strengthening the networks that bind ecosystems together is crucial. Prioritizing these areas can boost their resilience to the ever-looming threat of climate change and fortifies the stability of our environment as a whole [68].

The Getis analysis results identified hot spots and cold spots with confidence levels of 90%, 95%, and 99%. In hot spot and cold spot analysis, results with 90% confidence are considered insignificant. Furthermore, the overlay of three ESs revealed that the distribution of hot spots, with frequencies of 3, 2, and 1, corresponds to high, medium, and low conservation priorities, respectively. Higher frequency indicates a greater priority for designating the area as a conservation zone. This is because such areas typically maintain good ES conditions, which implies lower conservation cost allocation [69]. The Getis-Ord Gi* statistics are a crucial tool in this analysis as it evaluates spatial clustering by comparing local sums of feature values to their neighbors, ultimately determining statistical significance through z-scores and p-values [70].

Hot spots deliver vital ES like water filtration, climate regulation, and carbon storage, which are essential for human well-being and survival. Incorporating advanced techniques, such as species distribution modeling and remote sensing, can further improve the identification and management of these critical regions [71]. Preserving these areas ensures the ongoing provision of these services, benefiting both ecosystems and human communities. Directing conservation efforts toward identified hot spots enables the protection of numerous species within a concentrated area, making resource use more efficient [72].

While hot spots are key focal points, it is important to acknowledge that biodiversity also exists beyond these areas. Effective conservation strategies should address not just hot spots but also the broader ecological networks to achieve comprehensive biodiversity protection. Therefore, prioritizing hot spots can help sustain global biodiversity and support local communities reliant on these ecosystems [73].

To avoid oversimplified conclusions that could undermine conservation outcomes, ES mapping requires careful, nuanced interpretation. By prioritizing on overlapping ES hot spots, resource-limited conservation efforts can achieve maximum efficiency across well-connected landscapes [74]. Modeling fidelity, which involves measuring discrepancies between actual spatial ES patterns and model outputs, underscores the importance of environment-specific adjustments. Furthermore, adaptation and adjustment are essential to address the complexities inherent in diverse ecological and urban contexts [75]. Discrepancies in models highlight the importance of environment-specific adjustments to accommodate complex ecological and urban contexts [76].

Identifying cumulative impacts on ES through mapping is essential for pinpointing high-impact areas and informing conservation priorities. This technique reveals regions that are critical for sustaining ES flows and should be incorporated into management strategies [77].

An ES hot spot can denote areas with high levels of a single service or regions providing multiple services. Variations in spatial configurations between hot spot identification methods can introduce uncertainties in decision-making. These differences can also impact the analysis of spatial overlaps between multiple service hot spots and between ES and biodiversity [12].

Conservation strategies highlight the need to prioritize biodiversity hot spots to achieve the best outcomes, especially under resource constraints. This targeted approach ensures that vital ecosystems receive the necessary protection [78].

The policy for implementing conservation area planning prioritizes areas: conservation area planning prioritizes high-value areas with significant ES to minimize costs and maximize ecological benefits. Integrating biodiversity and ES goals enhances efficiency and effectiveness as targeting both can provide substantial benefits without significantly increasing biodiversity losses. This approach emphasizes the importance of focusing on high-priority areas to achieve optimal ecological and economic outcomes [79].

This study introduces a novel approach by integrating advanced spatial techniques to assess and prioritize ES conservation at a finer spatial resolution. Unlike previous studies that mainly focus on individual ES or broad-scale analyses, this research examines the compactness and efficiency of overlapping ES hot spots and cold spots, specifically for water yield, soil retention, and CS. It provides a comprehensive framework to identify spatial co-occurrence and variations in ES provision through a combination of the InVEST model, Getis-Ord (G*) statistics, and ArcGIS overlay technique.

The results showed that: (1) The study area is dominated by the following LULC types: paddy field covering 216,141 ha (31.78%), dry farming covering 205,412 ha (30.20%), plantation forest covering 76,854 ha (11.30%), and settlement areas covering 72,264 ha (10.62%). The image analysis achieved a kappa coefficient of 0.84 and an overall accuracy of 88%. (2) The total water yield amounts to 77.93 × 10⁸ m³ year^-1, SC amounts to 31.02 × 10³ tons year^-1, and CS amounts to 31.65 × 10⁵ tons year^-1. (3) WY, SC, and CS are statistically significant with confidence levels exceeding 95%, covering 137,078 ha (19.84%), 86,433.59 ha (12.51%), and 41,455 ha (6.00%) of the study area, respectively. Based on the ES levels, conservation priority areas were classified as high, medium, and low priority, encompassing 14,124 ha (2.04%), 38,268 ha (5.54%), and 137,133.47 ha (19.85%) of the study area, respectively. Based on the overlay of two conservation area maps using the Three ES and the TES approaches in the Citarum Watershed, both methods showed similarity in results over an area of 620,202 hectares (89.77%) and a difference in results covering 70,714 hectares (10.23%). This study provides a framework for setting conservation priorities using environmental criteria, aiming to guide decision-makers in identifying areas that require protection.