Spatially Explicit Systematic Conservation Planning Based on Multi-Ecosystem Services in an Urban Tropical Watershed: A Case Study of the Citarum Watershed, Indonesia

2.1 Study area

This study focuses on spatial modeling of ecosystem service-based conservation areas in the Citarum River Basin in West Java Province, as shown in Figure 1. The Citarum River Basin in West Java Province, covering approximately 690,000 ha, is undergoing rapid change and faces significant development pressures.

The watershed is located between 5º54'11" and 7º15'3" S latitude, and between 106º56'31" and 107º58'33" E longitude. Its pronounced elevation gradient influences climatic variability, hydrological processes, and land-use patterns across the watershed. The region’s topography is varied, with hills and volcanic formations among its topography. Variations in slope can be found at the base (5 to 15%), on mountain slopes (15 to 30%), and at the peaks (30 to 90%)

The area experiences a dry climate every three months, with an average annual rainfall of 2,358 millimeters. Three major dams, Saguling, Cirata, and Jatiluhur, are located along the Citarum River, which flows through the watershed and plays an essential role in freshwater supply, agriculture, and electricity generation in West Java. The region has varied topography, including hills and volcanic formations. Slope variations ranging from 5-15% in the lower areas, 15-30% on mountain slopes, and 30-90% at the peaks.

tu_pian_1.png

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

The Citarum Watershed is divided into three geomorphological regions: upstream, midstream, and downstream. The upstream region consists of mountainous volcanic terrain at elevations ranging from 750 to 2,600 m above sea level [14]. This sub-watershed is dominated by Latosol, Andosol, and Regosol soils [15], with an average annual rainfall of approximately 4,000 mm and a minimum temperature of 15.3 ℃. The midstream region is also located within volcanic formations, with elevations ranging from 200 to 800 m above sea level, annual rainfall between 1,000 and 4,000 mm, and temperatures ranging from 15.3 ℃ to 27 ℃. The downstream region is a flat lowland area at elevations of 1 to 200 m above sea level. This area is characterized by Alluvial, Entisol, and Inceptisol soils [15], with average annual rainfall of around 1,000 mm and a minimum temperature of 27 ℃.

2.2 Data sources

This study used ecosystem service maps, land use and land cover (LULC) data, climatic variables, soil data, and conservation-area information to assess ecosystem service dynamics and support conservation planning in the Citarum Watershed. Multi-temporal LULC maps for 2010 and 2020 were derived from Landsat 5 TM and Landsat 8 OLI/TIRS imagery (path/row 122/65; 30 m spatial resolution) obtained from the USGS EarthExplorer platform. The images were classified using a supervised Random Forest algorithm to generate LULC maps at a spatial resolution of 30 m × 30 m.

Topographic information was obtained from the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) provided by the United States Geological Survey (USGS). Watershed boundaries were compiled from maps provided by the College of Forestry, Environment and Resources Management, Ministry of Environment and Forestry of the Republic of Indonesia, and the Citarum Ciliwung River Basin Center, and were further processed through digital watershed delineation using the DEM.

Climate data, including rainfall and temperature records, were obtained from the Indonesian Agency for Meteorology, Climatology and Geophysics (BMKG), and the Citarum Ciliwung River Basin Center. Spatial climate surfaces were generated using spline interpolation. Evapotranspiration data were obtained from the WorldClim database and processed using the same interpolation approach. Soil data, including soil texture, organic matter content, and effective rooting depth, were obtained from the Citarum Ciliwung River Basin Center and converted from polygon to raster format through extraction and resampling procedures.

These datasets were used to model WY, SR, and CS. CS was estimated using Indonesian land-cover carbon pool values. Detailed information on the datasets, processing procedures, and data formats is provided Appendix 1. Data processing, and spatial analysis were conducted using several software tools, including GeoDa, ArcGIS 10.8 and Marxan [16-18].

2.3 Research methodology

The InVEST model was used to simulate three ES, namely WY, SR, and CS. Spatial relationship pattern analysis was then applied to identify the spatial distribution, trade-offs, and synergies among the ES. SCP was subsequently employed to delineate conservation zones in the Citarum Watershed. The conceptual framework of Ecosystem Service-based SCP is presented in Figure 2.

tu_pian_2.png

Figure 2. Conceptual framework of ecosystem service-based Systematic Conservation Planning (SCP) for the Citarum Watershed

Multiple ecosystem service modeling. This stage focused on the quantification and spatial mapping of ecosystem services across the Citarum Watershed. The InVEST model was used to simulate three key ES, namely WY, SR, and CS, for 2010 and 2020. These services represent important hydrological, geomorphological, and climate-regulating functions related to watershed sustainability. The outputs were generated as spatially explicit raster layers, which provided a biophysical basis for subsequent spatial analysis and conservation planning. This stage established the ecological foundation of the study by identifying the magnitude and spatial distribution of ecosystem service provision.

Although the InVEST model provides a practical framework for ecosystem service assessment, the outputs remain subject to uncertainty due to input data limitations, parameter assumptions, and model simplifications [17]. Therefore, the consistency and plausibility of the modeled ecosystem service outputs were assessed using several complementary approaches.

For WY, modeled outputs were qualitatively compared with available hydrological observations from runoff monitoring stations located at the watershed outlet as a general plausibility check. However, due to the limited availability of consistent long-term field measurements, comprehensive field-based validation could not be fully conducted. For SR and CS, direct field validation was not possible because in situ observational datasets were unavailable.

To provide a quantitative assessment of model consistency, a benchmarking exercise was conducted by comparing the modeled ecosystem service values with previously published estimates for the same watershed reported by the study [19]. This benchmarking analysis was performed at the sub-watershed level using the following statistical indicators:

1. Root Mean Square Error (RMSE), calculated as the square root of the average squared differences between the modeled values and the study values [19], expressed in the same units as each ecosystem service;

2. Coefficient of determination (R²), used to evaluate the degree of linear correspondence between the modeled outputs and the reference dataset;

3. Percent deviation, calculated as: Percent Deviation = (|Modeled − Reference| / Reference) × 100%.

It is important to note that this benchmarking exercise represents a model-to-model comparison and consistency assessment rather than true field validation. Therefore, the resulting statistical metrics primarily reflect the degree of agreement between datasets generated using similar methodological frameworks and input sources. Future studies should incorporate field-based measurements and independent observational datasets to improve the robustness and empirical reliability of ecosystem service assessments.

The modeling was carried out using InVEST version 3.14.8, specifically the Annual Water Yield (AWY), Sediment Delivery Ratio (SDR), and CS modules. The WY module estimates the spatial distribution of WY based on the water balance principle. In this approach, WY is defined as the portion of precipitation remaining after plant transpiration and surface evaporation within each grid cell. The model assumes that the remaining water eventually reaches the watershed outlet through surface and subsurface runoff. The WY estimation incorporates several factors, including precipitation, potential evapotranspiration, root system characteristics, and soil depth. The model outputs were further adjusted using observed runoff data from hydrological stations located at the basin outlet. The final WY values for each grid were calculated using the algorithms outlined in Eqs. (1) and (2) [17].

$Y_{x j}=\left(1-A E T_{x j} / P_x\right) \times P_x$                            (1)

$\frac{A E T(x)}{P(x)}=\frac{1+w_x R_{x j}}{1+w_x R_{x j}+1 / R_{x j}}$                            (2)

where, $Y_{x j}$ represents the WY of grid cell $x$ with land-use type $j\left(\mathrm{~m}^3 \cdot \mathrm{hm}^{-2}\right), A E T_{x j}$ represents the annual actual evapotranspiration of grid cell $x$ with land-use type $j(\mathrm{~mm}), P_x$ represents the annual average precipitation of grid cell $x(\mathrm{~mm})$, $R_{x j}$ represents the dryness index of grid cell $x$ with land-use type $j$, and $w_x$ represents the available water content for vegetation.

SR was assessed using the SDR module in the InVEST model. In this approach, SR was calculated as the difference between potential soil erosion and sediment export for each pixel. Potential soil erosion was calculated using the RKLS factor, while actual soil erosion was estimated using the Universal Soil Loss Equation (USLE). The calculations are presented in Eqs. (3) to (6) [17].

$R K L S_i=R_i \times K_i \times L S_i$                     (3)

$U S L E_i=R_i \times K_i \times L S_i \times C_i \times P_i$                     (4)

$S R_i=R K L S_i-S E D_i$                    (5)

$S R_i=R_i \times K_i \times L S_i-\left(U S L E_i \times S D R_i\right)$                        (6)

where, $S R_i$ represents SR at pixel $i, R K L S_i$ denotes potential soil erosion, calculated from rainfall erosivity ($R_i$), soil erodibility ($K_i$), and slope length-steepness factors ($L S_i$); $U S L E_i$ represents actual soil erosion after considering landcover management $\left(C_i\right)$ and conservation practice factors $\left(P_i\right)$.

$S D R_i$ represents the sediment delivery ratio, indicating the proportion of eroded soil transported from each pixel to the stream network, which reflects landscape connectivity an topographic conditions. Meanwhile, $S E D_i$ refers to sedimen export, calculated as $U S L E_i \times S D R_i$, representing the amoun of soil effectively delivered to the river system.

CS was estimated using the InVEST CS module. CS within the landscape was calculated based on land-use data from different periods and their corresponding carbon pool values. Total carbon stock was computed as the sum of aboveground, belowground, soil, and dead organic matter carbon, as shown in Eq. (7) [20, 21]:

$C S_{\text {total}}=C S_{\text {above}}+C S_{\text {below}}+C S_{\text {soil}}+C S_{\text {dead}}$                      (7)

where, total CS ($\mathrm{CS}_{\text {total}}$) consists of aboveground carbon ($\mathrm{CS}_{\text {above}}$), belowground carbon $\left(\mathrm{CS}_{\text {below}}\right)$, soil carbon $\left(\mathrm{CS}_{\text {soil}}\right)$, and carbon stored in dead organic matter $\left(\mathrm{CS}_{\text {dead}}\right)$. All values are expressed in ton$\cdot \mathrm{ha}^{-1}$. The analytical procedures followed the guidelines described in the research [17].

The InVEST model outputs were validated by comparing the results with those reported in a previous study [19]. Due to the limited field observation data, model performance was evaluated using the RMSE, which measures the deviation between modeled and reference values. Lower RMSE values indicate better model performance.

Spatial autocorrelation of multiple ecosystem services. This stage examined the spatial structure and interaction patterns among ESs. Spatial autocorrelation analysis was applied to identify clustering patterns, spatial dependence, and potential synergies or trade-offs among ESs. Using spatial statistical approaches, the analysis assessed whether ESs were spatially clustered, dispersed, or randomly distributed. This step helped identify areas with overlapping or concentrated ES and provided important information for conservation prioritization.

The ecosystem service outputs, initially generated in raster format, were converted into vector format with a spatial resolution of 500 m × 500 m. The study area was then divided into planning units (PUs), resulting in a total of 28,302 PUs. Ecosystem service values for each PU were extracted using the zonal statistics tool in ArcGIS [22]. The values were subsequently standardized using Eq. (8):

$E S_i=\frac{E S_{i, o b s}-E S_{i, \min}}{E S_{i, \max}-E S_{i, \min}}$                          (8)

where, $E S_{i, \text { obs}}$ represents the observed ecosystem service value, while $E S_{i, \text { max}}$ and $E S_{i, \text { min}}$ denote the maximum and minimum values, respectively. $E S_i$ refers to the standardized value obtained from the normalization process.

Spatial autocorrelation analysis was conducted using GeoDa software. Univariate LISA was used to identify spatial patterns within individual ES, while bivariate analysis was applied to examine relationships between ES based on standardized values in each planning unit. Following the principles of Exploratory Spatial Data Analysis (ESDA), the analysis was used to detect spatial dependence, heterogeneity, and distribution patterns in the spatial data [23]. LISA and multivariate Local Moran’s I were further implemented to highlight areas with significant spatial autocorrelation. Global Moran’s I and Local Moran’s I were calculated for the entire study area using GeoDa 1.18 (http://geodacenter.github.io).

A first-order queen contiguity spatial weights matrix was used, where planning units sharing either a boundary or a vertex were considered neighbors. Binary weights were assigned and subsequently row-standardized so that the weights for each planning unit summed to one, minimizing the influence of unequal numbers of neighbors. The spatial weights matrix was generated in GeoDa 1.18. No isolated planning units (islands) were identified after planning-unit generation; therefore, no special treatment for zero-neighbor observations was required.

Interactions among ESs in the Citarum Watershed, including trade-offs and synergies, were examined using multivariate Global Moran’s I to assess the presence and strength of spatial autocorrelation at the grid level [24]. The formulas for Global Moran’s I and Local Moran’s I are presented in Eqs. (9) and (10):

Global Moran’s I:

$I=\frac{n \sum_{i=1}^n \sum_{j=1}^n W_{i j}\left(x_i-\bar{x}\right)\left(x_j-\bar{x}\right)}{\sum_{i=1}^n \sum_{j=1}^n W_{i j}\left(x_i-\bar{x}\right)^2}$                   (9)

Local Moran’s I:

$I=\frac{n^2}{\sum_{i_i} \sum_j W_{i j}} \cdot \frac{\left(x_i-\bar{x}\right) \sum_j W_{i j}\left(x_j-\bar{x}\right)}{\sum_j\left(x_j-\bar{x}\right)^2}$                        (10)

In these equations, $n$ represents the total number of spatial unit samples within the study area. The variables $\mathrm{x}_{\mathrm{i}}$ and $\mathrm{x}_{\mathrm{j}}$ denote the attribute values of spatial units $i$ and $j$, respectively, while $\overline{\mathrm{x}}$ refers to the mean attribute value across all units. The spatial weight matrix $\mathrm{W}_{\mathrm{ij}}$ represents the spatial relationship between unit $i$ and unit $j$, reflecting their relative spatial proximity. Moran's I and LISA were applied to evaluate both global and local autocorrelation of ecosystem service values across the planning units.

The Z-score and p-value were used to assess whether global autocorrelation was positive or negative. The sign of Moran’s I further indicates whether the spatial pattern is clustered or dispersed. A positive Moran’s I value indicates positive spatial autocorrelation, meaning that similar values tend to cluster spatially, with higher values indicating stronger spatial relationship and clearer global clustering patterns. In contrast, a Moran’s I value close to 0 suggests a random spatial distribution without a clear pattern. Local spatial autocorrelation analysis was conducted in GeoDa using Anselin Local Moran’s I cluster and outlier analysis to identify spatial relationships among ES values across the study area.

To further examine trade-offs and synergies among ES, global and local bivariate Moran’s I analyses (bivariate LISA) were applied. Global bivariate Moran’s I measures the overall spatial correlation between two ES across the study area, whereas local bivariate Moran’s I identifies spatial correlations within individual spatial units. The equations used in this study are presented in Eqs. (11) and (12) [25].

$I_{e u}=\frac{N \sum_i^N \sum_{\neq i}^{n N} W_{i j} z_i^e z_j^u}{(N-1) \sum_i^N \sum_{j \neq i}^N W_{i j}}$                            (11)

$I_{e u}^{\prime}=z^e \sum_{j=1}^N W_{i j}$                          (12)

where, $I_{e u}$ and $I_{e u}^{\prime}$ represent the global and local bivariate Moran's I values for ES, respectively; $N$ is the total number of spatial units; $W_{i j}$ is the spatial weight matrix describing the spatial relationship between spatial units $i$ and $j$. The matrix was generated using queen contiguity weights to first-order neighboring units [25]. The variable $z_i^e$ and $z_j^u$ represent the standardized values of ES $\mathrm{ES}_{\mathrm{i}}$ and $\mathrm{ES}_{\mathrm{j}}$.

The values of $I_{e u}$ and $I_{e u}^{\prime}$ range from -1 to 1. Positive values indicate positive spatial correlation, meaning that areas with high ecosystem service values tend to be surrounded by areas with similarly high values. Conversely, negative values indicate negative spatial correlation, where areas with high values are surrounded by areas with low values. Larger absolute values indicate stronger spatial relationships between ESs. Statistical significance was assessed using a permutation test with 999 permutations [26], and spatial correlations were considered significant at p < 0.05.

Bivariate LISA was used to visualize local spatial relationships through Moran's scatterplots, cluster maps, and significance maps. The four quadrants of bivariate LISA scatterplot represent four types of local spatial autocorrelation: high-high (HH), high-low (HL), low-high (LH), and low-low (LL). HH clusters indicate synergy, where high values of one ecosystem service are associated with high values of another service. HL and LH clusters represent trade-offs, where high values of one service are associated with low values of another. Meanwhile, LL clusters indicate areas where both ES are low, reflecting zones of ecological degradation or low ecosystem performance. Identifying these spatial interaction patterns helps distinguish areas with strong conservation synergies (HH) from areas with potential land-use conflicts (HL and LH).

Spatial modeling for watershed conservation planning. At this stage, ES information was translated into spatially explicit conservation priorities. The mapped ES layers were incorporated into a Marxan-based SCP framework as conservation features. Cost layers representing land-use constraints were also included to reflect implementation feasibility. Conservation targets were defined for each ecosystem service, and spatial optimization was performed to identify priority conservation areas that were both cost-efficient and spatially coherent. This process enabled ecosystem service information to be operationalized into conservation zoning within the Citarum Watershed. Three ecosystem service layers, namely WY, SR, and CS, were incorporated into the Marxan model as individual conservation features. In addition, the three layers were combined into a single integrated feature representing multi-ES, referred to as Aggregated Ecosystem Services (AES).

The conservation cost surface was developed based on the level of human disturbance associated with different land-use types and the total economic value of ES, which were converted into relative conservation costs. Conservation status was derived from the forest area map [27]. Detailed classifications of conservation costs and conservation status categories are presented in Tables 1 and 2.

Parameter calibration was conducted using the QMarxan platform. The primary parameters adjusted in the analysis included the Species Penalty Factor (SPF), Boundary Length Modifier (BLM), and the number of iterations and runs. Calibration of SPF and BLM is essential to balance conservation target achievement and the spatial compactness of the resulting conservation network. This approach helps generate conservation solutions that are not only cost-efficient but also spatially connected and ecologically effective as expressed in Eqs. (13) and (14) [18, 28]:

Table 1. Conservation cost feature

Table 2. Conservation status feature

$Total$ $Cost$ = $basic$ $cost$ + $fragmentation$ $penalty$ + $unmet$ $target$ $penalty$                  (13)

$\begin{aligned} & \text {Total Cost}=\sum_{i=1}^n \text {Cost}+(\text {BLM} \times \left.\sum \text {Boundary }\right) \sum_{i=1}^n(\text {SPF} \text {× penalty})\end{aligned}$              (14)

where, $\sum_{i=1}^n \operatorname{Cost}$ represents the total costs of all planning unit derived from area or economic value. The term BLM $x \sum$ Boundary represents the boundary length modifier component, which regulates the importance of maintaining compact conservation areas. Higher BLM values impose stronger penalties on fragmented solutions (those with long boundaries), thereby encouraging the formation of contiguous conservation areas. Meanwhile, $\sum_{i=1}^n$ (SPFxPenalty) represents the penalty applied when conservation targets are not achieved. Higher SPF values indicate greater conservation priority for specific features.

During the data input process, penalties were assigned proportionally based on boundary length and the additional costs required to achieve unmet conservation targets. Increasing SPF values generally improve the likelihood of achieving conservation targets [29]. However, previous studies have noted that no fixed theoretical basis for determining optimal SPF values [30]. In practice, SPF values greater than 1 are commonly recommended, and some studies have used values up to 100 to ensure target achievement [29, 31].

To improve the robustness of the Marxan optimization, sensitivity analysis was conducted for key parameters, including SPF and BLM. Multiple parameter combinations were tested to evaluate their effects on spatial configuration, compactness, and conservation target achievement [32]. The results showed that conservation priorities remained relatively stable across a reasonable range of parameter values. In this study, the BLM value was set to 1. The ecosystem service protection targets evaluated in this study ranged from 20% to 50%. The 50% ecosystem service protection target was not intended to represent a direct land-area conservation target, but rather a precautionary conservation threshold used to identify spatial priorities capable of maintaining ecosystem service resilience under future environmental pressures. The final conservation configuration generated under the 50% ecosystem service target produced approximately 28–30% conservation area coverage, which aligns with the global “30 by 30” conservation target under the Kunming–Montreal Global Biodiversity Framework. The use of a relatively high ecosystem service target was intended to increase spatial flexibility and improve the robustness of conservation prioritization in anticipation of future land-use change and ecological degradation

Marxan employs a simulated annealing optimization algorithm, which increases the likelihood of generating near-optimal spatial configurations compared with conventional heuristic approaches [33]. The model identifies planning units with the high selection frequencies that can be prioritized for conservation, as well as areas that are consistently excluded from the solution [18]. Spatial weighting is an important component of the zoning process because it can reduce management costs and minimize potential land-use conflicts by reducing spatial overlap among competing land uses [34].

Each ES raster was incorporated into the Marxan model as a conservation feature. A conservation target equivalent to 30% of the total ES supply was defined for each service, reflecting Indonesia's conservation policy targets. Two planning scenarios were evaluated: (A) incorporating existing protected areas as locked-in units, and (B) excluding existing protected areas to represent a fresh planning approach. Four Marxan runs were conducted for each scenario, and the solution with the lowest total cost was selected as the optimal conservation configuration.

Spatial modeling of conservation areas under two planning scenarios. Conservation area modeling was conducted under two planning scenarios. In this study, the impacts of LULC change and climate change were not explicitly simulated using process-based ecosystem service models. The protection targets of 20%, 30%, 40%, and 50% refer to Marxan feature-representation targets rather than percentages of watershed area. Each target specifies the minimum proportion of the total amount of each ecosystem service (WY, SR, CS, and aggregated ES) that must be represented within the selected planning units. Because Marxan identifies the minimum-cost combination of planning units that satisfies these representation targets while maximizing spatial complementarity, the resulting conservation area is substantially smaller than the feature target itself. Consequently, achieving a 50% ecosystem service representation target required approximately 29% of the watershed area rather than 50% of the land area.

To assess the effectiveness of protected-area representation, two planning approaches were applied [27]. In scenario A, existing protected areas were fully integrated into the planning process, with all formally designated protected areas locked into the Marxan simulations. This scenario was designed to identify additional priority conservation areas beyond the current protected area network while retaining existing protected areas in the final conservation configuration. Retaining these areas in the output maps also enabled comparison between newly identified conservation priorities and the formally designated protected zones. In scenario B, existing protected areas were not enforced in the planning process. As a result, all planning units had an equal opportunity to be selected as conservation areas. This approach represents a fresh planning scenario without spatial constraints imposed by the current protected area network.

After running both scenarios, the selection frequency of each planning unit across the 1,000 iterations was analyzed. Selection frequency was used as an indicator of irreplaceability, where planning units selected more frequently were considered more important for achieving conservation objectives. The classification of conservation priority levels followed [27, 35], as defined below:

• Priority I (>75%) – planning units with very high importance that form the core components of the conservation network.
• Priority II (50–75%) – planning units with considerable importance that should be considered for conservation designation.
• Priority III (25–50%) – planning units functioning as complementary areas to strengthen the conservation network.
• Priority IV (<25%) – planning units with lower conservation importance but still offering potential contributions to conservation efforts.

In this study, priority conservation areas were defined based on Priority I units, representing planning units with the highest conservation importance.

Spatial pattern similarity analysis of conservation areas. In this study, four levels of ES protection targets, namely 20%, 30%, 40%, and 50%, were evaluated, with each scenario executed using 1,000 iterations. To assess the effectiveness of protected-area representation under different planning scenario, two main approaches were applied [27, 32].

To quantify the spatial similarity between Scenario A and Scenario B, JSD was calculated from the normalized probability distributions of planning-unit selection frequencies generated by Marxan. For each planning unit, the selection frequency obtained from 1,000 Marxan runs was divided by the sum of selection frequencies across all planning units, producing a probability distribution for each scenario (P for Scenario A and Q for Scenario B). Thus, the JSD compares the complete spatial distributions of planning-unit selection probabilities rather than binary conservation maps or total conservation area. A JSD value of 0 indicates identical spatial priority distributions, whereas a value approaching 1 indicates increasingly different conservation priorities.

This method was selected because it is symmetric, produces bounded values, and can be applied to distributions containing zero values [36]. Each planning unit (PU) was assigned a relative probability value based on its ES score, as expressed in Eq. (15):

Relative Probability $=\frac{\mathrm{ES} \text {value of}\ \mathrm {PU}}{\text {Maximum ES value across PUs}} \times 100 \%$                    (15)

The probability assigned to planning unit $i$ was calculated as

$p_i=\frac{S F_i}{\sum_{j=1}^n S F_i}$                         (16)

where, $S F_i$ is the Marxan selection frequency obtained from 1,000 runs. The resulting normalized vector $\left(\mathrm{p}_{\mathrm{i}}\right)$ was used as the input probability distribution for the JSD calculation. The selection frequency of each PU across all iterations was compiled to form a spatial probability distribution representing the likelihood of that unit being selected in an optimal conservation solution. The probability distributions from the two planning scenarios were then compared using the JSD approach. The Kullback–Leibler (KL) divergence used in the calculation is expressed in Eq. (17) [37]:

$K L[P(x) \| Q(x)]=\sum_{x \in X}\left[P(x) \log \frac{P(x)}{Q(x)}\right]$                    (17)

where, $P(x)$ represents the probability distribution derived from one scenario, while $Q(x)$ represents the probability distribution from the alternative scenario. The term $X$ denotes the set of all planning units within the study area, and $\log \frac{P(x)}{Q(x)}$ represents the logarithmic ratio used to measure the relative dissimilarity between the two distributions.

The JSD was then calculated using Eq. (18) [37]:

$\begin{gathered}J S D(P(x) \| Q(x))=\frac{1}{2} K L \left(P(x) \| \frac{P(x)+Q(x)}{2}\right)+\frac{1}{2} K L\left(Q(x) \| \frac{P(x)+Q(x)}{2}\right)\end{gathered}$                      (18)

where, JSD $(\mathrm{P}(\mathrm{x}) \| \mathrm{Q}(\mathrm{x})$) represents the Jensen-Shannon divergence between distributions P and Q , while $\frac{\mathrm{P}(\mathrm{x})+\mathrm{Q}(\mathrm{x})}{2}$ represents the mixed distribution M , defined as the average of the two distributions. The factor $\frac{1}{2}$ is used because the JSD is calculated as the average of the two KL divergence terms. The two KL terms measure the divergence of each distribution relative to the mixed distribution.

The JSD values range from 0 to 1. A value close to 0 indicates that the two spatial distributions are highly similar, whereas a value close to 1 indicates substantial spatial differences. In this study, JSD values below 0.3 were interpreted as similar patterns, values between 0.3 and 0.6 as moderately different patterns, and values above 0.6 as significantly different patterns. This approach provides a quantitative assessment of the stability and consistency of conservation area selection generated by the Marxan simulations [29]. In the final stage, the conservation-area maps from both scenarios were overlaid in ArcGIS. The overlay analysis was used to quantify the percentage overlapping conserved and non-conserved areas, enabling spatial comparison between the two planning scenarios.

4.1 Spatial autocorrelation and interactions among ecosystem services

The results indicate that ecosystem service dynamics in the Citarum Watershed were strongly influenced by land-use change between 2010 and 2020. WY and CS declined by 15.08% and 7.16%, respectively, whereas SR increased by 12.46%. These contrasting trends suggest that land-use transformation altered the balance among hydrological regulation, soil conservation, and carbon sequestration functions across the watershed. Similar patterns have been reported in tropical watersheds where urban expansion and agricultural intensification reduce ecosystem service capacity and ecological resilience [19, 43].

The contrasting ecosystem service dynamics observed between 2010 and 2020 reflect differences in the ecological processes governing hydrological regulation, soil conservation, and carbon sequestration. WY declined primarily because increased vegetation cover and restoration activities enhanced evapotranspiration and infiltration, reducing surface runoff, while decreasing rainfall further limited water availability [44]. CS also decreased as a consequence of forest conversion and continued land-use intensification, which reduced aboveground biomass and soil carbon pools [4, 5]. In contrast, SR increased because improvements in vegetation cover, conservation practices, and reduced sediment connectivity enhanced the landscape's capacity to trap eroded soil before it reached the stream network. Since SR is influenced by land management, slope characteristics, and hydrological connectivity rather than biomass alone, it can improve even when CS declines. Similar contrasting responses among ES have been reported in tropical watersheds where restoration measures effectively reduce soil erosion despite ongoing losses of mature forest cover [2, 20].

In contrast, SR increased because it is controlled not only by biomass, but also by vegetation cover, land management practices, topography, and sediment connectivity. Localized improvements in vegetation cover and soil conservation measures, particularly in upstream areas, may have enhanced the landscape's capacity to retain eroded soil before it reached the river network. Consequently, SR can increase even when total CS declines, as these ES respond to different controlling factors.

Spatial autocorrelation analysis showed that WY, SR, and CS were dominated by high–high and low–low clusters, indicating strong spatial synergies among ES. High-value ecosystem service clusters were concentrated primarily in upstream and midstream areas, suggesting that these regions function as ecosystem service hotspots capable of simultaneously supporting hydrological regulation, soil conservation, and CS. Such co-occurrence reflects the influence of common biophysical drivers, including vegetation cover, topography, and landscape condition, including vegetation cover and topographic characteristics. Similar spatial synergy patterns have been reported in other watershed systems, supporting the concept of ecosystem service bundles [45].

Although synergistic relationships dominated, localized trade-offs were also identified, particularly between WY and vegetation-related ES. Areas with higher vegetation density generally exhibited improved CS and SR but may experience reduced WY due to increased evapotranspiration [46]. These results highlight the importance of spatially differentiated management strategies because conservation actions that maximize one ecosystem service may not always maximize others.

From a management perspective, the concentration of multiple ES in upstream and midstream areas suggests that conservation and restoration efforts should prioritize these regions. Protecting ecosystem service hotspots can simultaneously support several ecological functions, thereby improving conservation efficiency and long-term watershed sustainability [47, 48].

4.2 Scenario-based spatial optimization of conservation areas

The apparent difference between the 50% ecosystem service target and the resulting conservation coverage (29.21 % of the watershed) reflects the optimization principle implemented in Marxan. ES are spatially heterogeneous and highly concentrated in ecologically valuable planning units, particularly in the upstream forested landscape. Consequently, many selected planning units simultaneously satisfy multiple ecosystem service targets (feature complementarity), allowing Marxan to achieve high ecosystem service representation while minimizing total conservation area and implementation costs. This result demonstrates that ecosystem service representation targets should not be interpreted as equivalent to land-area conservation targets.

The 20–50% feature targets were evaluated as sensitivity scenarios to assess the response of conservation solutions to increasing ecosystem service representation requirements. The final analysis adopted the 50% ecosystem service representation target because it generated a stable and spatially connected conservation network while producing approximately 29.21 % watershed coverage, closely matching the international 30 × 30 conservation objective.This wording resolves the apparent contradiction and aligns the methodology with standard Marxan practice.

These findings suggest that existing protected areas provide an important spatial foundation for conservation planning. By serving as anchors within the optimization process, they enhance connectivity and reduce the amount of additional land required to meet conservation objectives. Similar observations have been reported in SCP studies emphasizing the role of protected areas in reducing fragmentation and improving conservation efficiency [29].

However, the results also indicate that the current protected-area network alone is insufficient to achieve higher ecosystem service conservation targets. Additional conservation zones are required, particularly in areas with high ecosystem service values outside the existing network. This finding supports the need for integrated conservation strategies that combine protected-area strengthening with targeted expansion into priority non-protected areas [49, 50]. The ability of Scenario A to approach the global "30 by 30" conservation target demonstrates the practical value of ecosystem service-based SCP for supporting biodiversity and sustainability objectives at the watershed scale [11, 51].

4.3 Spatial divergence and multi-ecosystem service optimization

The high JSD = 0.78 demonstrates substantial differences between the normalized planning-unit selection-frequency distributions generated under Scenarios A and B. This result indicates that incorporating existing protected areas strongly influences the spatial allocation of conservation priorities. Importantly, the JSD compares the probability distributions of planning-unit selection frequencies across all Marxan runs rather than simple overlap of the final conservation maps. Consequently, the high divergence reflects different optimization pathways and conservation priorities instead of disagreement in overall conservation targets.

Rather than indicating planning inconsistency, the observed divergence reflects the existence of multiple spatially viable conservation solutions within the watershed. Scenario A produced a more connected conservation network because existing protected areas functioned as ecological anchors that improved landscape continuity. In contrast, Scenario B identified alternative conservation opportunities outside the existing protected-area system, thereby increasing spatial flexibility but producing a more fragmented configuration [52].

The high spatial divergence also highlights the importance of spatial complementarity in SCP. Conservation efficiency is achieved not solely by selecting areas with the highest ecosystem service values, but by identifying combinations of planning units that collectively maximize ecosystem service representation while maintaining connectivity and reducing redundancy [53]. Consequently, JSD can be interpreted as an indicator of scenario sensitivity and alternative conservation pathways rather than simply a measure of spatial inconsistency [13, 54].

4.4 Implications for this study and policy

This study demonstrates how multiple ES can be integrated into a SCP framework to support watershed-scale conservation decision-making. Unlike conventional approaches that focus primarily on biodiversity representation, the proposed framework incorporates WY, SR, and CS simultaneously, providing a broader representation of ecosystem functions and socio-ecological processes. The results further show that incorporating existing protected areas improves spatial coherence and connectivity while reducing the need for extensive land conversion [55].

Despite its practical value, the framework remains predominantly biophysical and does not explicitly incorporate social, economic, governance, or climate-change dimensions. Future studies should integrate stakeholder preferences, economic valuation, and alternative land-use and climate scenarios to improve the applicability and long-term robustness of conservation planning outcomes [47].

The findings have direct implications for land-use governance and watershed management in the Citarum region. Integrating ecosystem service priorities into regional spatial planning instruments, such as RTRW and RDTR, could strengthen ecological connectivity and reduce future land-use conflicts. The identified priority conservation areas, particularly in upstream and midstream regions, provide important guidance for protecting hydrological regulation, soil conservation, and CS functions.

The ecosystem service trade-off and synergy maps generated in this study can support conservation and restoration planning. Areas characterized by strong ecosystem service synergies represent priority conservation zones, whereas areas with consistently low ecosystem service performance may be prioritized for rehabilitation and restoration interventions [45]. In addition, mechanisms such as Payments for Ecosystem Services (PES) may help encourage conservation practices on private and community-managed lands while promoting stakeholder participation and long-term sustainability.

Although developed for the Citarum Watershed, the proposed framework has potential applicability in other tropical watersheds experiencing rapid land-use change and ecosystem degradation. Future applications should emphasize local parameter calibration, higher-resolution datasets, and site-specific validation to improve predictive accuracy and support broader implementation.

This research was supported by an Internal Research Grant from Universitas Esa Unggul (UEU). The study was conducted under the cooperation agreement between UEU and the National Research and Innovation Agency (BRIN), under Agreement No. 207/I/KS/12/2023 and No. 97/MOU/R/UEU/XII/2023, dated 01 December 2023. The authors acknowledge the support of the Research Center for Limnology and Water Resources BRIN, the Institute for Research and Community Service, and the Faculty of Engineering UEU, in facilitating this collaborative research.