Assessment of the Impacts of Climate and LULC Changes on the Water Yield in the Citarum River Basin, West Java Province, Indonesia

Abstract

Changes in climate and land use land cover (LULC) are important factors that affect water yield (WY). This study explores which factors have more significant impact on changes in WY, spatially and temporally, within the Citarum River Basin Unit (RBU), West Java Province, Indonesia with an area of ±11.317 km². The climate in the area of Citarum RBU belongs to the Am climate type, which is characterized by the presence of one or more dry months. The objectives of the study were: (1) To estimate a water yield model using integrated valuation of ecosystem services and tradeoffs (InVEST), and (2) to test the sensitivity of water yield (WY) to changes in climate variables (rainfall and evapotranspiration) and in LULC. The integration of remote sensing (RS), geographic information system (GIS), and the integrated valuation of ecosystem services and tradeoffs (InVEST) approach were used in this study. InVEST is a suite of models used to map and value the goods and services from nature that sustain and fulfill human life. The parameters used for determining the WY are LULC, precipitation, average annual potential evapotranspiration, soil depth, and plant available water content (PAWC). The results showed that the WY within the territory of Citarum RBU was 12.17 billion m3/year, with mean WY (MWY) of 935.26 mm/year. The results also show that the magnitude of MWY in Citarum RBU is lower than the results obtained in Lake Rawa Pening Catchment Areas, Semarang Regency and Salatiga City, Central Java (1.137 mm/year) and in the Patuha Mountain region, Bandung Regency, West Java (2.163 mm/year), which have the same climatic conditions. The WY volume decreased from 2006, to 2012, and 2018. Based on the results of the simulation, climatic parameters played a major role affecting WY compared to changes in LULC in the Citarum RBU. This model also shows that the effect of changes in rainfall (14.06–27.53%) is more dominant followed by the effect of evapotranspiration (10.97–23.86%) and LULC (10.29–12.96%). The InVEST model is very effective and robust for estimating WY in Citarum RBU, which was indicated by high coefficient of determination (R²) 0.9942 and the RSME value of 0.70.

1. Introduction

Water yield (WY) is an indicator of the health of a watershed. In a healthy watershed, water fluctuations between the rainy and dry seasons tends to be small [1]. The guaranteed quantity, quality, and continuity of water in a watershed is essential to the concept of water security, which directly or indirectly supports national food and energy security [2]. Quantitative evaluation and visualization of WY is valuable for understanding trends in the function of water supply in an ecosystem. Understanding WY is very useful for water resource managers to determine the effect of human activities on water resources [3].

In the last 20 years, environmental conditions and water quality along the Citarum River have declined significantly. Urban areas are hotspots that drive environmental change at multiple scales. Rapid urbanization, as a result of accelerated development, is linearly proportional to industrial activity, high rates of population growth, expansion of residential areas, and the conversion of land to built-up areas [4]. Various negative impacts arise as cumulative compensation for the imbalance between rapid economic development activity and environmental preservation [5].

The Citarum Watershed Pollution and Damage Control Team has been established in an effort to improve pollution and damage in the Citarum watershed in compliance with Presidential Regulation No.15 of 2018. Accurate estimation and calculation of the elements that affect WY are critical to determine the appropriate means to protect ecosystem services, such as revegetation techniques, and to meet water demand for socio-economic systems [6].

Previous research [7] stated that WY, especially those controlled by rainfall and evapotranspiration (ET), as well as land use change caused by humans, may indirectly affect WY. According to [8], changes in climate and land use/cover (LULC) resulting from human activities are the most critical factors that drive change in WY. Referring to [9], WY was greatly influenced by precipitation. The higher the precipitation, potentially the higher water yield. Meanwhile, other authors [10] have observed that there are several natural dynamic that can affect water yield: Land cover type, soil type, and land surface. Water in bare land cannot be stored effectively, and water will become a surface runoff. On sandy soil, water will easily be lost, because sandy soil has a low ability to hold water. The land surface also affects the ability of a land to hold water. Land that has a higher slope has a lower ability to hold water than flat land; whereas [11] stated the combined impact of climate change and LULC can also affect water yield.

The integrated valuation of ecosystem services and tradeoffs (InVEST) model has been widely applied globally, especially with respect to environmental service valuation [8,12,13,14]. This model, which was developed by the Natural Capital Project at Stanford University [15], can be used to calculate the WY from a watershed. InVEST is a tool that can be used to assess ecosystem/environmental services and support decision making in environmental management [12]. The WY model, one of the modules in InVEST, uses a water balance approach [16,17] and is based on the Budyko curve [18] and average annual rainfall.

The WY ecosystem modelled by the InVEST was used in this study to estimate spatial variation in WY capacity in the area of Citarum River Basin Unit (RBU), as a representative for the application of the InVEST Model in developing countries with tropical climates. The objectives of the study were: (1) To estimate a WY model using InVEST, and (2) to test the sensitivity of WY to changes in climate variables (rainfall and evapotranspiration) and in LULC.

2. Materials and Methods

2.1. Study Site

This study was conducted in the Citarum RBU, which covers 19 watersheds. The Citarum RBU lies between longitude 106°51′36″–107°51′ E and latitude 7°19′–6°24′ S, covering an area of ±11.317 km². The climate in the area of Citarum RBU belongs to the Am climate type, which is characterized by the presence of one or more dry months. Administratively, Citarum RBU extends over 13 Regency/City administrative areas in the West Java Province. The Citarum RBU is bordered by the Java Sea to the north, Cianjur Regency and Bandung Regency to the south, Garut Regency, Indramayu Regency, and Sumedang Regency to the east, and Sukabumi Regency, Bogor Regency, and Bekasi Regency to the west (Figure 1).

The topography of the Citarum watershed, described morphologically, can be grouped into three parts, namely the upstream zone, middle zone, and downstream zone. The upper Citarum RBU is a large basin and is known as the Bandung Basin, with an elevation range of 625–2600 m msl. The morphology of the central part of the Citarum watershed varies between plains (elevation of 250–400 m msl), weak wavy hills (elevation of 200–800 m msl), steep hills (elevation of 1400–2400 m msl), and volcanic bodies. The downstream part of the Citarum watershed is dominated by plains and weak and steep wavy hills with various elevations between 200 and 1200 m above msl. All rivers within the area of the Citarum RBU flow from south to north, upstream from Mt. Burangrang, Bukit Tunggul, and Canggah, and downstream to the north coast of the Java Sea [19].

2.2. Research Data and Tools

The integration of remote sensing (RS), geographic information system (GIS), and the integrated valuation of ecosystem services and tradeoffs (InVEST) approach were used in this study. The data used in this study (spatial and non-spatial data) were collected from relevant agencies; this included catchment and sub-catchment boundaries, LULC maps, precipitation (in mm), average annual potential evapotranspiration (in mm), soil depth (in mm), and plant available water content (PAWC) (percentage), in addition to LULC attributes.

The software used in this study included ArcGIS 10, InVEST, and SPSS. InVEST modelling is spatial based modelling. Raster format data with spatial resolution of 30 m × 30 m and WGS 84 system coordinates were used in the InVEST model. In general, the research was divided into two stages: Data preparation and data analysis. The flowchart has been modified from the flowchart developed by [20] (Figure 2).

2.2.1. Land Use/Land Cover (LULC)

LULC maps (2006, 2012, and 2018) were employed in this study. LULC data for 2006 and 2012 in shapefile format were collected from existing data developed by [21]. However, the LULC 2018 data were generated from analysis of processed images from Landsat 8 OLI (22 March 2018). The selection in initial year (2012) is based on the declaration of the Citarum SWS management with a new paradigm and relatively complete data availability. Considering that the selection of one period is 6 years, changes of forest in Landsat imagery (woody plants) can be easily detected. This is in line with the research conducted by [22] where the analysis of land cover using Landsat imagery was over a period of 6 years. All LULC data were assigned to 12 classes (

Table 1. Classes and total area of land use/cover (LULC) in Citarum River Basin Unit (RBU) in 2006, 2012, and 2018.

Code	LULC	2006		2012		2018
		Ha	%	Ha	%	Ha	%
1	Virgin Forest	36.70	3.24	29.07	2.57	27.98	2.47
2	Plantation Forest	118.98	10.51	121.39	10.73	116.14	10.26
3	Shrub	18.84	1.66	19.61	1.73	7.39	0.65
4	Estate Crops Plantation	57.17	5.05	56.73	5.01	51.02	4.51
5	Settlement Area	95.57	8.44	112.61	9.95	108.80	9.61
6	Bare land	10.34	0.91	9.98	0.88	7.80	0.69
7	Lake	16.42	1.45	16.39	1.45	16.40	1.45
8	Pure Dry Agriculture	125.88	11.12	317.55	28.06	151.56	13.39
9	Mixed Dry Agriculture	118.46	10.47	144.80	12.79	132.99	11.75
10	Paddy Filed	499.14	44.10	269.37	23.80	477.02	42.15
11	Fishpond	34.07	3.01	34.07	3.01	34.47	3.05
12	Airport	0.19	0.02	0.19	0.02	0.19	0.02
	Total	1131.75	100.00	1131.75	100.0	1131.75	100.00

).

The InVEST model requires a biophysical table containing information on LULC along with an appropriate code, crop coefficient (Kc), and root depth. The InVEST model does not use root depth information for land in classes without vegetation cover/use [15], so any value can be entered (in this study, a value of 1 was used). Vegetated LULC classes have a value of 1 and LULC classes lacking vegetation (freshwater, buildings, and settlements) have a value of 0.

2.2.2. Rainfall

Rainfall data for 2000–2018 were obtained from the Meteorological, Climatological, and Geophysical Agency (BMKG) [23], the River Basin Territory Organization or the Balai Besar Wilayah Sungai (BBWS) of Citarum Ciliwung [24], and PT Jasa Tirta II. Rainfall data were obtained from 35 rainfall observation stations located in Bandung, Cicalengka, Bandung Geophysics, Pamanukan, Cikao Bandung, Cikarang Dam, Bekasi Dam, Jatiasih, and Karang (Figure 3). Rainfall analysis was carried out for three periods: 2000–2006, 2006–2012, and 2012–2018. The average annual rainfall in the Citarum RBU ranged from 676 to 3894 mm/year for the 2000–2006 period, between 817 and 3446 mm/year for the 2006–2012 period, and between 789 and 3284 mm/year for the 2012–2018 period.

The average annual rainfall data for each period were used to create a rainfall map for 2006, 2012, and 2018. The spline spatial interpolation technique in ArcGIS was chosen to create a monthly rainfall map and annual rainfall map. The monthly rainfall map was used for calculating monthly reference evapotranspiration, while the annual rainfall map was used in the InVEST model analysis.

2.2.3. Annual Reference Evapotranspiration

The annual reference evapotranspiration map was compiled from analysis of extra-terrestrial solar radiation, minimum air temperature, maximum air temperature, and monthly rainfall. Daily extra-terrestrial solar radiation for each rainfall station was calculated using Microsoft Excel, then the values were combined to obtain the monthly value. Monthly extra-terrestrial solar radiation maps in raster format were generated using spline techniques. The amount of air temperature data from local meteorological stations that could be obtained was not sufficient to generate an air temperature map. Therefore, this study used the average minimum and maximum air temperature obtained from global climate data [25] at a spatial resolution of 1 km × 1 km. Then, the data were resampled to 30 m × 30 m, which is developed by [26].

The calculations to obtain monthly and annual reference evapotranspiration maps were performed using the ArcGIS raster calculator. Annual reference evapotranspiration values ranged from 792 to 1921 mm/year for 2006, from 764 to 1749 mm/year for 2012, and from 794 to 2039 mm/year for 2018.

2.2.4. Depth of Soil Solum and Plant Available Water Content

Soil maps in shapefile format were obtained from the Citarum BBWS. Based on these maps, there were 10 soil types in the study area: (1) Alluvial, (2) Regosol, (3) Latosol, (4) Andosol, (5) Grumusol, (6) Litosol, (7) Mediterranean, (8) Podzolic, (9) Resin, and (10) Gley humus. The depth of the soil solum was obtained from the land system map. For WY analysis, the soil solum depth map was converted into raster format and into mm.

The PAWC was obtained based on data of soil type and texture (percentage of clay/loam and sand fractions) using soil water characteristics software developed by the Agricultural Research Service USDA in collaboration with the Department of Biological Systems Engineering, Washington State University [27]. Due to the limited amount of data available for soil characteristics, default values were used for other parameters (e.g., percentage of organic matter and salinity). The map of available water capacity for plants was converted into raster format.

2.2.5. Watershed Boundaries

The data regarding the boundaries of watersheds in the Citarum River Basin Territory were obtained from the Citarum-Ciliwung BBWS in shapefile format. The Citarum River Basin Territory comprises 19 watersheds. This watershed map was entered as input to the InVEST model.

In this study, the separation of watershed boundaries was also carried out based on the Regulation of the Minister of Public Works of the Republic Indonesia Number 4/PRT/M/2015 concerning Criteria and Designation of River Basin, namely upstream zone, middle zone, and downstream zone river basin.

2.3. Data Analysis

2.3.1. Water Yield

Reference evapotranspiration (mm.day⁻¹), $E{T}_0{}_{\left(x\right)}$, was calculated using the modified Hargreaves equation [28], which gives better results than the Penman–Monteith method when the data required is given as follows:

$$E{T}_{0_{\left(x\right)}}=0.0013\times 0.408\times Ra\times \left[\left(\frac{T_{max}+T_{min}}{2}\right)+17\right]\times {\left[\left(T_{max}-T_{min}\right)-0.0123P\right]}^{0.76}$$

(1)

where:

• $Ra$ = extra-terrestrial solar radiation (MJ.m⁻².day⁻¹)
• $T_{max}$ = average maximum daily air temperature (°C)
• $T_{min}$= average minimum daily air temperature (°C)
• $P$ = monthly rainfall (mm.day⁻¹)

For non-vegetation LULC (e.g., water bodies or settlements), actual evapotranspiration was calculated directly from reference evapotranspiration, $E{T}_0{}_{\left(x\right)}$, and it has an upper limit determined by rainfall as follows:

$$AE{T}_{\left(x\right)}=Min\left(K_{c_{\left(x\right)}}\times E{T}_{0_{\left(x\right)}}, P_{\left(x\right)}\right)$$

(2)

$AW{C}_{\left(x\right)}$ determines the amount of water stored in the soil and released for use by plants. This parameter was estimated using PAWC, the minimum root restricting layer depth, and vegetation rooting depth as follows:

$$AW{C}_{\left(x\right)}=Min\left(Rest.layer.depth, root.depth\right)\times PAWC$$

(3)

$E{T}_0{}_{\left(x\right)}$ reflects local climatic conditions based on the evapotranspiration of reference plants at that location. Whereas $K_c{}_{\left(x\right)}$ is mainly determined by the vegetation characteristics of land use/cover at each pixel [29]. The coefficient of plant available water capacity at each pixel, ${\omega }_{\left(x\right)}$ [30], was calculated as follows:

$${\omega }_{\left(x\right)}=Z\times \frac{AW{C}_{\left(x\right)}}{P_{\left(x\right)}}+1.25$$

(4)

where:

• $AW{C}_{\left(x\right)}$ = the volume (mm) of plant available water capacity.
• $Z$ = an empirical constant (sometimes referred to as the seasonality factor/Zhang coefficient), reflecting the local precipitation pattern and additional hydrogeological characteristics. In this study, the Z value used was 4, which is the recommended value for watersheds in tropical areas [31].

The non-physical parameter that characterize the natural climatic-soil properties are both detailed below. Potential evapotranspiration, PET(x), was calculated as follows:

$$PE{T}_{\left(x\right)}=E{T}_0{}_{\left(x\right)}\times K_{c_{\left(x\right)}}$$

(5)

where:

• $E{T}_0{}_{\left(x\right)}$= the reference evapotranspiration at pixel $x$
• $K_c{}_{\left(x\right)}$ = the plant (vegetation) evapotranspiration coefficient at pixel x associated with its LULC.

For vegetated type of LULC [16], $\frac{AE{T}_{\left(x\right)}}{P_{\left(x\right)}}$ is estimated in a spatially explicit way on pixelx, that is:

$$\frac{AE{T}_{\left(x\right)}}{P_{\left(x\right)}}=1+\frac{PE{T}_{\left(x\right)}}{P_{\left(x\right)}}-{\left[1+{\left(\frac{PE{T}_{\left(x\right)}}{P_{\left(x\right)}}\right)}^{\omega }\right]}^{\frac{1}{\omega }}$$

(6)

where:

• $PE{T}_{\left(x\right)}$ = potential evapotranspiration for pixel $x$
• $\omega$ = non-physical parameters that characterize the correlation between climate and soil properties also called the coefficient of available water capacity for plants [12,32].

In this study, the annual WY for each pixel, $Y_{\left(x\right)}$, for a given LULC was determined as follows [18]:

$$Y_{\left(x\right)}=\left(1-\frac{AE{T}_{\left(x\right)}}{P_{\left(x\right)}}\right)\times P_{\left(x\right)}$$

(7)

where:

• $AE{T}_{\left(x\right)}$ = annual actual evapotranspiration for the pixel x
• $P_{\left(x\right)}$ = annual precipitation at pixel x.

Validation of the InVEST model was carried out on the total WY. The data were obtained from the Geospatial Information Agency [33]. Linear regression analysis was carried out between the actual observed data and the modelling results’ estimation data. Based on the analysis results, the coefficient of correlation (R²) and Pearson correlation (r), and root mean square error (RMSE) were obtained to determine the validation of the model. All statistical analyses were performed using SPSS software.

2.3.2. Impact of Changes in Climate and LULC on the WY

Changes in LULC and climate are the major drivers of changes in water yield. Climate change limits changes in rainfall and reference evapotranspiration. To find out which factors have the most significant impact on WY, four simulation models were run using the WY InVEST model (

Table 2. Parameter change for simulation model.

Parameter Change	Scenario Name	Description
Climate		Climate change simulation: It was assumed that there will be a climate change, while the LULC data used is data of 2006, 2012, and 2018.
	Climate 1	It was assumed that there was no change in climate inputs. The climate data inputs are rainfall data of in 2006, 2012, and 2018 and reference evapotranspiration data of 2006.
	Climate 2	It was assumed that rainfall increased by 10% and evapotranspiration were constant.
	Climate 3	It was assumed that there was an increase in rainfall and evapotranspiration by 10%.
	Climate 4	It was assumed that rainfall decreased by 10%, and evapotranspiration was constant.
	Climate 5	It was assumed that there was a decrease in constant rainfall and evapotranspiration by 10%.
LULC		Simulation of changes in LULC: There was a change in LULC. Meanwhile, the climate data remains unchanged.
	LULC 1	LULC was assumed to be unchanged, so the input data of LULC in 2006, 2012, and 2018 was LULC data in 2006.
	LULC 2	LULC was assumed to have changed. All industrial plantations were converted into open land (an increase in the area of open land by 124% (in 2006), by 129% (2012), and by 122% (in 2018). The average increase in the open area was 125%.
	LULC 3	LULC was assumed to have changed. All industrial plantations have been converted into open land (increase in open land area by 1217% (in 2006), by 3182% (2012), and by 1944% (in 2018)). The average increase in the open area was 2114%.
	LULC 4	LULC was assumed to have changed. All paddy fields have been converted into open land (increase in open land area by 4827% (in 2006), by 2699% (2012) and by 4548% (in 2018)). The average increase in the open area was 4548%.
	LULC 5	LULC was assumed to have changed. Shrubs have been converted into open land (an increase in open land area by 1217% (2006), 196% (2012), and 95% (2018)). The average increase in the open area was 158%.

).

3. Results

3.1. Water Yield in Citarum RBU

Based on the results of the analysis of the InVEST model, the volume of WY at RBU Citarum is around 12.17 × 109 m³/year, where this value is obtained from the average thickness of water 935.26 mm/year multiplied by the total area of the entire RBU (

Table 3. The WY in the Watershed and River Basin Territory in the Citarum River in 2018.

Watershed Area/WS (River Basin Territory)	Area		Mean WY	Total WY		WY Coefficient
	103 Ha	Percentage	Mm	109 m3/Year	Percentage
Catchment Area
Citarum	659.50	58.25	1220.70	8.05	66.18	0.63
Cipunara	128.06	11.31	1126.94	1.44	11.86	0.62
Ciasem	73.19	6.46	969.28	0.71	5.83	0.60
Cimalaya	52.06	4.60	974.30	0.51	4.17	0.58
Cikarokrok	36.33	3.21	656.34	0.24	1.96	0.52
Others	183.06	16.17	663.99	1.22	9.99	0.48
Total RBU	1132.20	100.00	935.26	12.17	100.00
WS (River Basin Territory)
Upstream	24.40	21.56	951.71	2.32	19.07	0.60
Middle	69.51	61.39	1084.48	7.53	61.89	0.60
Downstream	19.31	17.05	1201.63	2.32	19.05	0.62
Total RBU	1132.20	100.00	1079.27	12.17	100.00

). This WY reflects natural river flow, not taking into account the use of water in human activity, such as by households, industry, and agriculture [10,13]. The Citarum Watershed has the highest annual rainfall (1.994 mm/year) with the lowest potential evapotranspiration (1.291 mm/year). With actual evapotranspiration of 649 mm/year, the Citarum Watershed produces the highest WY at 1.220 mm/year (Table 3). The producer of the second-largest WY is the Cipunara Watershed (1.126 mm/year), the third-largest is contributed by the Cimalaya Watershed (974 mm/year), and the fourth-largest is from the Ciasem Watershed (969 mm/year).

In terms of volume, the Citarum Watershed, which is the widest watershed, makes the largest contribution to WY, with 8.05 × 10⁹ m³/year (66.18%), followed by the Cipunara Watershed with 1.44 × 10⁹ m³/year (11.89%), the Ciasem Watershed at 0.71 × 10⁹ m³/year (5.83%), and the Cimalaya Watershed with 0.51 × 10⁹ m³/year (4.17%).

With regard to distribution, the middle river basin produces the largest amount of water as it covers an area of 69.51 thousand ha (61.39%) with a total WY of 7.53 × 10⁹ m³/year (61.89%). The remainder of the yield contribution is evenly distributed between upstream and downstream. Meanwhile, on average, the highest WY comes from the downstream area at 1201 mm/year, the middle area provides 1,084 mm/ year, and the upstream area yields 951 mm/year. The spatial pattern of WY and rainfall is shown in Figure 4. Meanwhile, changes in WY from 2006 to 2018 are shown in

Table 4. Changes in WY for the period 2006–2018 in Citarum River in the period 2006-2018.

Catchment Area	Total WY 109 m3/Year			Change in Total WY 109 m3/year (% Class of Change)
	2006	2012	2018	2006–2012	2012–2018	2006–2018
Citarum	9.54	10.10	8.05	5.80 (LI)	(20.28) (MD)	(15.65) (LD)
Cipunara	1.69	1.79	1.44	5.92 (LI)	(19.59) (LD)	(14.83) (LD)
Ciasem	1.03	0.94	0.71	(8.83) (LD)	(24.81) (MD)	(31.45) (MD)
Cimalaya	0.94	0.73	0.51	(22.82) (MD)	(30.17) (MD)	(46.11) (HD)
Cikarokrok	0.26	0.36	0.24	37.68 (MI)	(34.27) (MD)	(9.50) (LD)
Others	2.38	1.81	1.22	(23.82) (MD)	(32.92) (MD)	(48.90) (HD)
SWS Citarum	15.86	15.74	12.16	(0.75) (LD)	(22.71) (MD)	(23.29) (DM)

Remark: LD = Low decrease (0–20%), MD = moderate decrease (20–40%) and high decrease (more than 40%), LI = low increase (0–20), MI = moderate increase (20–40%) and HI = high increase (more than 40%).

.

Figure 4 shows the average annual WY in Citarum RBU ranges from 470 to 1220 mm with an average value of 763 mm year⁻¹. About 50% of the watersheds in Citarum RBU have WY lower than the average WY, while the remaining watersheds (50%) have an average WY greater than the average WY. The WY in SWS Citarum in 2012 was relatively the same as in 2006. Overall, the WY in 2018 decreased by 23.23% However, when viewed from the conditions of each watershed, there were various changes. In 2006–2012, there has been an increase in WY in some watersheds, while in 2012–2018, there was a decrease in WY.

The increase in WY occurred in Citarum, Cipuana and Cikarokrok Watersheds in 2006–2012 period, while in the 2012–2018 there was a decrease. Meanwhile, the Cimalaya and Ciasem watersheds have always experienced a decrease in WY in the period 2006-2012 to the period 2012–2018. This condition indicates poor number of vegetation cover and is reflected by the changes in land cover. Results of land cover in the period 2006–2018 show that there has been a decrease in rice field area of 22,120 ha (4%), natural forest area of 8720 ha (23.75%), and plantations of 2830 ha (2.38%). There was an increase in residential area by 13,220 ha (13.84%), dry land agriculture covering 25,680 ha (20.40%), and mixed dry land farming covering an area of 14,530 (12.26%).

Spatially, the WY distribution pattern in the Citarum River Basin Territory follows the annual rainfall distribution (Figure 5).

The spatial pattern of WY and rainfall have a linear correlation, where lower rainfall is associated with a smaller WY, and vice versa, as shown in Figure 5c. According to [34], this value refers to the coefficient WY, which is the ratio between WY and precipitation per hectare for each type of LULC. The WY coefficient represents the amount of WY converted from precipitation, taking into account evapotranspiration, degree of saturation, and infiltration.

The WY coefficient ranged from 0.39 to 0.64, and varied between watersheds. Watersheds that produced high WYs (Citarum Watershed, Cipunara Watershed, Ciasem Watershed, and Cimalaya Watershed) had WY coefficients above 0.57. Meanwhile, the WY coefficient for the other watersheds was less than 0.57.

LULC in the study location was dominated by paddy field in 2018 (42.15% of the total area), pure dry agriculture (13.39%), settlements (9.61%), and plantation forest (10.46%).

Analysis of changes in LULC for the period 2006–2018 showed a reduction of 22,120 ha (4%) in the extent of the rice fields, a reduction of 8720 ha (23.75%) in virgin forest, and of 2830 ha (2.38%) in land used for plantations. In contrast, settlement areas had grown by 13,220 ha (13.84%), dryland agriculture areas by 25,680 ha (20.40%), and mixed dryland agriculture by 14,530 ha (12.26%). The agricultural and plantation land cover classes in the study sites are a mixture of timber crops (forestry plants and fruits) and agricultural crops. WY by land cover type (

Table 5. WY by LULC in the Citarum River Basin Territory, 2018.

Classes of LULC	Area		MWY	Total WY		Coefficient of WY
	Ha	Percentage	Mm	109 m3/Year	Percentage
Paddy Filed	479.41	42.36	985.01	4.72	38.83	0.59
Pure Dry Agriculture	151.85	13.42	1348.35	2.05	16.84	0.69
Mixed Dry Agriculture	133.98	11.84	1219.10	1.63	13.43	0.67
Plantation Forest	114.69	10.13	1226.80	1.41	11.57	0.61
Settlement Area	109.63	9.69	1018.73	1.12	9.18	0.64
Estate Crop	51.12	4.52	1196.45	0.61	5.03	0.61
Virgin Forest	26.67	2.36	1444.36	0.39	3.17	0.66
Shrubs	7.30	0.65	1516.25	0.11	0.91	0.78
Lake	16.51	1.46	549.50	0.09	0.75	0.24
Bare Land	7.85	0.69	424.08	0.03	0.27	0.22
Airport	0.19	0.02	898.61	0.00	0.01	0.61
Fish Pond	32.55	2.88	3.21	0.00	0.01	0.10
Total	1131.75	100.00	1075.20	12.16	100.00	0.54

).

Based on LULC, shrubs produce the highest average WY of 1516 mm/year (Table 5). The second-largest WY is from virgin forests (1444 mm/year), followed by pure dry agriculture (1348 mm/year), and plantation forests (1226 mm/year).

With respect to volume, paddy fields, which cover the most extensive area (42.36%), are the largest water producers with a total yield of 4.72 × 10⁹ m³/year (38.83%), followed by pure dry agriculture with 2.05 × 10⁹ m³/year (16.84%), mixed dry agriculture with 1.63 × 10⁹ m³/year (13.43%), and plantations at 1.41 10⁹ m³/year (11.57%).

For vegetated land cover, the WY coefficient ranges from 0.57 to 0.76. Paddy fields have the lowest ratio (0.57) and shrubs have the highest ratio (0.76). Meanwhile, the WY coefficient for non-vegetated land use ranges from 0 to 0.2. The difference inWY between vegetation types was analyzed by extracting raster pixel values MWY and mean annual precipitation (MAP) for the land use types virgin forest (VF), shrubs (SH), tea plantation (TP), settlement (ST), bare land (BA), agriculture (AG), and paddy field (PF). The relationship MAP with MWY values is shown in Figure 6.

Figure 6 illustrates the linear correlation between MAP and MWY in general terms. A high MAP will produce high and reversed MWY. Anomalous conditions occur in the land use type BA, where rainfall is relatively high but a low MWY is found. In order to explore the correlation between rainfall and WY, a paired linear regression analysis was performed between the MWY and MAP pixel values. Correlation and significance for MWY and MAP are shown in

Table 6. Correlation and significance for MWY and mean annual precipitation (MAP).

Type of Land Cover	R2	p-Value	Significance
Forest	0.97	0.00	0.00
Bare land	0.97	0.00	0.00
Tea plantation	0.94	0.00	0.00
Agriculture	0.93	0.00	0.00
Paddy field	0.79	0.00	0.00
Shrubs	0.68	0.00	0.00
Settlement	0.40	0.00	0.00

.

Referring to Table 6, WY and rainfall are positively correlated, with a correlation coefficient of greater than 0.79, except for land use types SH and ST. It can be concluded that land use types VF, BA, TP, and AG demonstrate a reasonably healthy correlation, while for SH and ST, the correlation between rainfall and WY is weak.

Data for WY from the Citarum RBU published by [33] were used to validate the InVEST model. By pairing this WY data with the WY modelling data, the correlation coefficient (R²) and Pearson correlation can be obtained, using the equation Y = 0.8682x + 0.2798, (R² = 0.9885). The RSME was also calculated.

3.2. Impact of Changes in Climate and LULC on theWY

The simulation results for changes in climate and LULC (

Table 7. The WY with climate scenarios.

Year	Normal	Climate 1		Climate 2		Climate 3		Climate 4		Climate 5
	109 m3	109 m3	Percen-tage	109 m3	Percentage	109 m3	Percen-tage	109 m3	Percentage	109 m3	Percentage
2006	15.85	15.85	0	18.08	14.06	17.59	10.97	13.65	(13.92)	14.12	(2.21)
2012	15.70	16.30	3.82	18.67	18.90	17.40	11.05	13.64	(13.16)	14.05	(2.07)
2018	12.17	14.93	22.68	15.52	27.57	15.07	23.86	10.96	(9.89)	11.47	(1.20)

and

Table 8. The WY with LULC scenarios.

Year	Normal	LULC1		LULC2		LULC3		LULC4		LULC5
	109 m3	109 m3	%	109 m3	%	109 m3	%	109 m 3	%	109 m 3	%
2006	15.85	15.85	0	16.67	5.18	16.27	3.25	14.97	(5.60)	16.92	6.71
2012	15.70	15.36	(2.23)	15.34	(2.33)	13.81	(12.04)	14.62	(6.90)	15.58	(0.80)
2018	12.17	12.10	(0.56)	11.47	(5.73)	10.91	(10.29)	8.66	(28.81)	12.09	(0.60)

). The WY under normal conditions was used as a reference for the amount of change in WY. Under normal conditions, the WY in 2012 is relatively similar to that in 2006. However, the 2018WY, compared to 2006 and 2012, decreased by 23.23%. In general, the 2018WY decreased compared to the initial condition (in 2006).

Water yield is driven by LULC. Changes in LULC for the period 2006–2018 showed a reduction in paddy field area by 22,120 ha (4%), in areas of natural forest by 8720 ha (23.75%), in plantation forest by 2830 ha (2.38%), and an increase in the area occupied by settlement of 13,220 ha (13.84%), by dryland agriculture of 25,680 ha (20.40%), and by mixed dryland agriculture by 14,530 ha (12.26%). This transformation contributes to changes in the WY.

The results of the Climate 2 and Climate 3 scenarios show an increase in WY by 10.96–27.57%, compared to the initial conditions. However, in the Climate 3 scenario WY does not increase significantly (relatively the same as Climate 2). The reduced climatic factor, as shown in the results of Climate 4 and Climate 5, has an impact on the decreasing WY. The Climate 5 scenario shows that reference evapotranspiration has no significant impact on decreasing WY.

In general, the condition of normal WY in 2018 has fallen compared to the initial condition (2006). The WY in the LULC1 scenario (without changes in LULC) demonstrated the same pattern as that under normal conditions, where WY decreased in 2012 and 2018. The decrease in WY ranged 0.56–2.23%. However, the simulation which included the addition of open areas of 1.58–45.46% (from vegetated land converted to open land) shows a decrease in WY of around 0.06–28.81%. The correlation between changes in area of types of LULC (%) and changes in WY (%), based on type of land cover, is shown in

Table 9. Correlation and significance for LULC areal Change for annual WY variations.

Type of LULC	R2	p-Value	Significance
Forest	0.9871	0.0000	0.000
Bare land	0.9963	0.0000	0.000
Tea plantation	0.9879	0.0000	0.000
Agriculture	0.9922	0.0000	0.000
Paddy field	0.9 856	0.0000	0.000
Shrubs	0.9819	0.0000	0.000
Settlement	0.4614	0.0000	0.000

.

Table 9 shows a significant positive correlation between change to areas of forest and annual WY (p < 0.000) in the study area. The correlation output matched the results shown in Table 6. The highest positive correlation is that between changes in WY and changes in area of bare land, with a correlation coefficient of 0.9963.

Conversely, the lowest correlation coefficient is that between change in WY and the area occupied by settlements. With a Pearson correlation value of 0.9819–0.9963, a p value of 0.0000, and significance of 0.0000, it can be concluded that there is a close correlation between changes in area of land cover and WY values, except when the land cover type is residential area.

4. Discussion

4.1. Water Yield in Citarum RBU

The WY in Citarum RBU is 12.17 × 10⁹ m³/year. The MWY is 935 mm/year, which is in line with the published data [33,35] that show the value of WY at Citarum RBU is 12.95 × 10⁹ m³/year and MWY by 994 mm/year. The magnitude and spatial distribution pattern of the WY from the modelling results are comparable to the [33] data, except that the WY value is smaller.

According to [33,35], the potential water resource in Citarum RBU is very large, amounting to 12.95 billion m³/year, supported by high rainfall of 2000–4000 mm/year [33]. The water required from the Citarum RBU yield is only by 7.65 × 10⁹ m³/year. The water requirement comprises irrigation, 6.63 × 10⁹ m³/year (86.7%), clean water, 0.46 × 10⁹ m³/year (6%), industry, 0.15 × 10⁹ m³/year (2%), urban water (municipal), 0.02 × 10⁹ m³/year (0.3%), and maintenance, 0.38 × 10⁹ m³/year (5%). The remaining 5.3 × 10⁹ m³/year is accounted for by potential water that has not been utilized (wasted at sea). The main irrigation facilities and infrastructure, in the form of dams, that currently function to supply water needs in the Citarum RBU are the Saguling Reservoir (for irrigation/hydropower), Cileunca Reservoir (hydropower), Cipanjuang Reservoir (hydropower), Cirata Reservoir (hydropower), and Jatiluhur (irrigation/hydropower).

Figure 5c and Table 6 show the WY of each watershed and the type of LULC, which is closely related to rainfall. The Pearson correlation value is R² = 0.79–0.95 with a p-value of 0.000, which indicates that there is a correlation between WY and rainfall. The findings for the Citarum RBU were in accordance with the results of [36], who found that when rainfall changes over time, the WY of a watershed also changes significantly, and there is a high correlation (R² = 0.954) between the two variables. Other research [37] reported that there is a simple linear correlation between annual rainfall data as an independent variable and the estimated volume of the WY in the same period as the dependent variable, with R² being 0.9992–0.9999. According to [38], the gradient distribution of MWY is consistent with MAP, i.e., MWY increases with MAP, in almost all ranges ofmean annual temperature (MAT).

The WY coefficient values in Citarum RBU ranged from 0.39–0.64 for each watershed and from 0.00–0.76 based on the type of land use. This coefficient is very similar to the Liang study results in the Qinghai Watershed, China, which are 0.00–0.82, but there are differences for each type of land cover. The yield coefficient value in the Citarum RBU area for vegetated areas is higher than bare land and built-up areas. However, other findings [36] in the Qinghai Lake region, China and [39] in the Jing-Jin-Ji region, China found that the WY coefficients for areas with buildings and for bare land are higher than the WH coefficient for vegetated areas. Built-up land is normally covered with asphalt, cement, and concrete, forming an impermeable layer that may reduce infiltration time and concentration. An increase in built-up land leads to an increase in WY. This condition is more appropriate in describing the WH coefficient for surface water (surface runoff), not the cumulative WY coefficient. The authors of the findings regarding WY in Qinghai and Jing-Jin-Ji, China, although they do not explicitly state that this is surface water, imply that the increase in WY is more accurately defined as surface water, as confirmed by [40]. Runoff is usually considered to be the quantity of water in a hydrological system that represents the movement of water on the earth’s surface.

The limitations of the InVEST model means that it is unable to distinguish between surface water and groundwater and is less sensitive to natural variability, so the hydrological cycle cannot be interpreted correctly, which allows the InVEST model to give different results. An application of the InVEST Model in South Ecuador by [41] estimated the annual sub-watershed and watershed surface runoff.

The problem of findings of different WH coefficients was identified by [42], who found that the correlation between LULC and water resources is complicated and difficult to predict due to the natural variability of watersheds, difficulties in controlling changes in LULC, and studies based on catchment areas where control is limited. According to [43], LULC is the key parameters to be considered in the study of WY, besides soil texture, surface runoff depth, stakeholders’ priorities, and stream order.

These findings are consistent with the results of the [1] study, with states that WY is an indicator of the health of a watershed, and a healthy watershed tends to have small water fluctuations. The decrease in WY in a period of 10 years is by 0–20% (low decrease), 20–40% (moderate decrease), and more than 40% (high decrease). Low-level decreased WY indicates good watershed management, moderately decreased WY indicates moderate watershed management, while highly decreased WY indicates poor watershed management.

Referring to [1], the management of Citarum, Cipuana, and Cikarokrok Watersheds are classified as good, the management for Ciasem Waterdheds is classified as moderate, and the management of the Cimalaya watershed and others are poor. Overall, the management of SWS Citarum is moderate.

The results of the validation of the modelling results with the observation data showed a Pearson Correlation value of 0.9885, an RSME value of 0.70, a p value of 0.0005, and significance of 0.0000. These results indicate that InVEST modelling of WY can be used to predict the actual WY conditions in Citarum RBU. Based on the spatial pattern, the distribution of WY from each watershed was relatively similar. The four largest watersheds also provided the largest WY. This finding is consistent with research by [44], who found that the total WY in each basin is influenced by MWY and area, while MWY distribution is closely related to rainfall.

4.2. Impact of Changes in Climate and LULC on theWY

Referring to Table 7, climate change is linearly correlated with WY, so an increase in climate variables will affect the increase in WY and vice versa. Increased rainfall and evapotranspiration gave relatively similar results compared to an increase in rainfall only, suggesting that evapotranspiration has a relatively small effect. This shows that the reference evapotranspiration factor has no significant effect on WY increase or decrease.

The results of a simulation with a 10% increase in rainfall predicted an increase in WY in Citarum RBU by 14.05–27.57%. However, the impact of the increase or decrease on evapotranspiration is not clearly visible, producing only a small effect (1.30–2.21%) as shown by simulation Climate 5 scenario, so this result can be ignored. This finding is consistent with the findings of [45], who found that WYs from 22 watersheds in the UK are very sensitive to changes in rainfall, a 10% increase in rainfall resulted in an increase in WY by 11–27%, but WY is not sensitive to variation in evapotranspiration.

The simulation showed that WY decreased by 0.56–2.23% under stable land cover conditions. This shows that under stable land cover conditions, the changes in WY are small. However, WY under stable climate conditions increased by 3.82–22.68%. The impact on WY under stable land cover conditions is relatively small when compared to stable climate conditions. This finding is consistent with [36] who stated that compared to land use/cover changes, rainfall has a bigger impact on WY. The study in China’s Qinghai Lake Watershed shows that the impact of land use change are much smaller than the impact of climate change (rainfall). According to [44], the magnitude of the change in WY depends on the size of the area converted and the type of land cover. The total WY of each basin is influenced by area, MWY, and MWY distribution, which is closely related to rainfall (MAP).

In general, the addition of open areas by 16,316–470,325 ha (158–4548%) (agricultural land converted into open land, plantations converted into open land, rice fields converted into open land, and shrubby areas converted into open land) caused a decrease in WY of 0.80 × 10⁹ m³ year − 13.50 × 10⁹ m³ per year (0.60–28.81%).

The results of the simulation of land cover change on WY in Citarum RBU showed that deforestation affected WY reduction. Previous authors [46] found that a decrease in WY is caused by deforestation, which is in accordance with the results of research in the Citarum watershed, where the impact of deforestation include decreased discharge, a substantially increased runoff coefficient, and reduced low flow during the dry season. According to [47], loss of forest causes a decrease in its function, marked by increased fluctuations in river flow in the dry season and the rainy season or heavier currents, increased flooding, and decreased reservoir capacity due to elevated sedimentation.

The results of deforestation and urbanization were clarified by [48] in the Brantas Watershed, Java (claiming to represent regional land change patterns in developing countries in Southeast Asia). The modelling results produced from the Soil and Water Assessment Tool (SWAT) show that reduction in forest cover and increased urbanization results in moderate changes to long-term water runoff (+8%), WY (+0.28%), and reduction in groundwater (−1.8%), and evapotranspiration (−1.15%). The SWAT Model can compensate for the InVEST model’s limitations, in that it cannot differentiate between surface water, subsurface water, and bottom flow. Therefore, according to [48], the impact of changes in land use on WY can only be seen in the long term.

The same change may have different impacts on WY, according to [49], whose research revealed that several findings do not show natural variability consistently. For example, deforestation increases WY [36,39,50] in China, in Indonesia [51], in Australia, the USA, Africa, and Germany [52], in Pakistan [53], in Malaysia [54], in Africa [55,56], and in Ecuador [41].

When WY is estimated using the InVEST model [36,39,55,56], in our opinion, the finding that WY has increased is misleading; in reality, the increase is in the quantity of surface water, not the total WY. The findings in Qanghai and Jing-Jin-JI, China, clearly stated that built-up land is normally covered with asphalt, cement, and concrete that form an impermeable layer, reducing infiltration times and concentrations, which are closely related to increased surface water. An increase in built-up land leads to an increase in WY, indicating an increase in surface water. This is made clear by the research of [56] who assessed WY in Africa using the InVEST model and clearly stated that an increase in surface water will eventually cause water scarcity and food insecurity. This was also supported by [41], who reported that the application of the InVEST Model in South Ecuador can estimate the annual runoff WY. This finding is in accordance with research [57] that the disconnection of the runoff and sediment delivery was confirmed by the reduction in the runoff delivery at plot scale due to the control of the length of the plot (slope) on the runoff and sediment delivery.

The application of the InVEST model in Jing-Jin-Ji, China, used a scale of 100 m × 100 m (especially on a scale of 1 km × 1 km), which is a medium and global scale where one pixel represents a vast area (1–100 ha). This produces different results because a region’s representation becomes more globalized; for example, a very remote station will generate relatively inaccurate interpolation. The model is unable to adequately record the hydrological process so the water results obtained are different.

Meanwhile, deforestation in Pakistan, Malaysia, and Indonesia correlates with the homogeneous forest type. Certain types of vegetation (coniferous forest in Pakistan and mercury pine forest in Kedungbulus, Central Java) can absorb more or less groundwater [51,52]. The absorption capacity of groundwater by forests is more determined by forest density, coniferous forest type, and terrain, so that an increase in forest cover (area density and vegetation) can be associated with a cumulative decrease in WY.

5. Conclusions

In this study, the WY was analyzed using the InVEST Model. LULC and climate changes have been simulated to determine the main factors driving changes in WY. InVEST modelling results show that the WY volume within the area of Citarum RBU in 2006 was 15.85 × 10⁹ m³/year, in 2012 it was 15.70 × 10⁹ m³/year, and in 2018 was around 12.17 × 10⁹ m³/year. Based on the validation of the InVEST model, these results indicate that the InVEST model can be used in estimating WY in the Citarum RBU as well as rule model for the application of the InVEST Model in developing countries with tropical climates. The climate in the area of Citarum RBU belongs to the Am climate type, which is characterized by the presence of one or more dry months with rain intensity <60 mm in 1 month. Referring to the magnitude of changes in WY in one period, it shows that the Citarum, Cipuana and Cikarokrok watersheds are classified as good management, the Ciasem watershed is classified as moderate management, and the Cimalaya watershed and others are poorly managed. Overall, the management of SWS Citarum is moderate.

The study also found that LULC for the period 2006–2018 showed a reduction of 22,120 ha (4%) in the extent of the rice fields, a reduction of 8,720 ha (23.75%) in virgin forest, and of 2,830 ha (2.38%) in land used for plantations. In contrast, settlement areas had grown by 13,220 ha (13.84%), dryland agriculture areas by 25,680 ha (20.40%), and mixed dryland agriculture by 14,530 ha (12.26%).

Based on the results of the InVEST simulation, it is known that climate change is a major factor affecting WY compared to changes in LULC in the Citarum watershed. This model also shows that the effect of changes in rainfall (14.06–27.53%) is more dominant followed by the effect of evapotranspiration (10.97–23.86%) and LULC (10.29–12.96%).