Groundwater recharge modeling with integration of land use/land cover and climate change projections in Surakarta City, Indonesia

Abstract

Increased population mobility in urban areas drives higher water demand and significant changes in Land Use and Land Cover (LULC), which directly impact groundwater recharge capacity. This study aims to predict LULC changes in 2030 and 2040, analyse groundwater recharge quantities for historical, current, and projected conditions, and evaluate the combined impacts of LULC and climate change. The Cellular Automata-Artificial Neural Network (CA-ANN) method was employed to predict LULC changes, using classified and interpreted land use data from Landsat 7 ETM+ (2000 and 2010) and Landsat 8 OLI (2020) imagery. The Soil and Water Assessment Tool (SWAT) model was used to simulate groundwater recharge. Input data for the SWAT model included Digital Elevation Model (DEM), soil type, LULC, slope, and climate data. Climate projections were based on five Regional Climate Models (RCMs) for two time periods, 2021–2030 and 2031–2040, under Shared Socioeconomic Pathways (SSP) scenarios 2–45 and 5–85. The results indicate a significant increase in built-up areas, accounting for 71.08% in 2030 and 71.83% in 2040. Groundwater recharge projections show a decline, with average monthly recharge decreasing from 83.85 mm/month under SSP2-45 to 78.25 mm/month under SSP5-85 in 2030, and further declining to 82.10 mm/month (SSP2-45) and 77.44 mm/month (SSP5-85) in 2040. The expansion of impervious surfaces due to urbanization is the primary factor driving this decline. This study highlights the innovative integration of CA-ANN-based LULC predictions with climate projections from RCMs, offering a robust framework for analysing urban groundwater dynamics. The findings underscore the need for sustainable urban planning and water resource management to mitigate the adverse effects of urbanization and climate change. Additionally, the methodological framework and insights gained from this research can be applied to other urban areas facing similar challenges, thus contributing to broader efforts in groundwater conservation.

Graphical abstract

Keywords

Groundwater Recharge / Climate Change / Remote Sensing / Socioeconomic Pathways / SWAT

Introduction

Among Asian countries, Indonesia has demonstrated a remarkable development record over the past 20 years is remarkable, according to the World Bank. Indonesia's population reached about 275 million in 2022 and is expected to grow to 330 million by 2050. In general, humans worldwide rely heavily on groundwater for daily water supply. Many researchers predict that freshwater resources will become scarcer in the future due to climate change and increasing human demand (Boretti and Rosa, 2019; Guo et al. 2024; Zhang et al. 2024). Groundwater resources can be assessed based on the existing storage capacity, the rate of inflow into the aquifer and the rate of outflow from the aquifer system. Groundwater aquifer systems are replenished primarily through rainfall recharge and the interaction between surface water and groundwater. Groundwater recharge processes are likely to be affected by land use and land cover change scenarios, as such changes can alter hydrological processes within the earth's surface (Adhikari et al. 2020). The absence of land use policies dramatically contributes to groundwater degradation in urban areas (Sajjad et al. 2022).

The rapid growth of the city's population has led to urbanization (Desky, 2022). This growth results increased urban density as more people migrate to cities. In addition to population growth, urban development has spatial impacts, including greater demand for land to meet the needs of the expanding population (Yunus, 2005). Consequently, land use and land cover in the city will also be affected. The expansion of built-up areas increases impermeability, limiting surface water infiltration into groundwater and thereby reducing groundwater availability.

Analysis of Land Use and Land Cover (LULC) changes using temporal and spatial data is essential for monitoring areas experiencing LULC changes by observing visual patterns and estimating changes in land area (As-syakur, 2008). This study employs Cellular Automata-Artificial Neural Network (CA-ANN) modeling to project LULC changes in Surakarta City in 2030 and 2040.

The research by Xu et al. (2019) showed that the CA-ANN model has the highest accuracy compared to Logistic Regression (LR) and Analytical Hierarchy Process (AHP). This is most likely due to the non-linear relationship inherent in urban growth, which cannot be accommodated by LR (Mustafa et al. 2018). A CA model comparison study conducted by (Ridho Fariz et al. 2020) also found that, when evaluated based on the total number of false alarms and missed detections, the CA model built by ANN outperformed LR.

Surakarta City has experienced rapid land cover change, and development activities continues to increase. The city's strategic location has led to improved accessibility and supporting infrastructure. Along with the rising population and population mobility, the demand for built-up land is also increasing. This directly affects hydrological conditions, particularly groundwater in Surakarta City. Water demand in the City has also risen significantly due to growing urbanization. In 2016, total water demand was 1,163.428 liters/second, increasing to 1,286.763 liters/second by 2022 (DLH Kota Surakarta, 2021). There is a high dependency on groundwater, yet groundwater recharge has been impacted by population growth and land use changes. The effect of land use change on groundwater recharge needs to be widely acknowledged by the public, as groundwater remains a vital source of freshwater for domestic, agricultural and commercial purposes. Without proper recognition and management, this could lead to the depletion of groundwater resources (Mengistu et al. 2022).

Spatially, Surakarta City has a critical water catchment area of 2,090 ha, accounting for 48% of the city's total area (DLH Kota Surakarta, 2022). Extensive land use changes have contributed to the deteriorating condition of infiltration areas in Surakarta City. To determine the quantity of groundwater infiltration under both current and future conditions, this study employed the SWAT (Soil and Water Assessment Tools) model. The SWAT model integrates multiple sub-models to predict river flow discharge over long periods in complex watersheds with diverse soil types, land uses, and management practices (Santhi et al. 2008). It is particularly suited for analyzing the impacts of climate change and land use change on hydrological processes within Hydrologic Response Units (HRUs).

Moreover, SWAT has been successfully applied in urban settings to estimate hydrological components, including groundwater recharge (Tamm et al. 2023; Ferreira and Uagoda, 2017; Piyumi et al. 2021; Halefom et al. 2017; Sisay et al. 2017). Previous studies have shown that even small-scale urban transformations can lead to disproportionate impacts on hydrological responses, reinforcing the necessity for city-scale modelling (Adhikari et al. 2020; Indhanu et al. 2025; Siddik et al. 2022). Applications of the SWAT model to project groundwater recharge in response to urbanization or land use changes can be found in Guo et al. (2008) and Zhang et al. (2016). Adhikari et al. (2020) also projected groundwater recharge in Ho Chi Minh City, Vietnam, based on climate change and Land Use and and Land Cover (LULC).

This research contributes to urban groundwater recharge studies in tropical regions by incorporating region-specific climate change projections into groundwater modeling. This study utilizes the CA-ANN model for LULC projection, which combines the spatial pattern recognition capability of Cellular Automata (CA) with the non-linear learning ability of Artificial Neural Networks (ANN). This integration allows for more accurate predictions of urban land use transitions by considering neighborhood effects and spatial dependencies - features not present in conventional models. Therefore, the primary objective of this study is to quantify the impact of increasing impervious built-up areas on groundwater recharge in Surakarta City. In the absence of previous research, this study offers future land cover and groundwater recharge projections using the CA-ANN and SWAT models, respectively.

1 Study area

The research location is in Surakarta City, located in Central Java Province. Geographically, it lies between 110° 45' 15" E and 110° 45' 35" E and 7° 36' S and 7° 56' S. Surakarta City covers an area of 46.72 km² and is surrounded by highlands to the West, East, and South. To the west, the city is bordered by Mount Merbabu and Mount Merapi; to the south, by the Sewu Mountains, which extend from Yogyakarta to Pacitan Regency; and to the east, by Mount Lawu. These surrounding highlands place Surakarta City within a basin zone, due to its topography being flanked by elevated terrain.

Surakarta City has an average elevation of 92 meters above sea level (masl), with a slope gradient of <20% and experiences average daily temperature ranging from 21°C to 31°C. The area lies in the humid tropics, maintaining warm temperatures throughout the year. Between 2000 and 2020, the average annual precipitation was approximately 1,936.4 mm/year, while the average annual potential evaporation was around 1,700 mm/year. According to data from the Central Bureau of Statistics (BPS), the total population of Surakarta City is 586,166, with a population density of 12,585 persons/km².

The spatial characteristics of Surakarta City, including geology, groundwater flow, elevation, and its delineation within the Bengawan Solo watershed using the SWAT model, are illustrated in Fig. 1. The city's geological featurs are covered by the Surakarta-Giritontro Sheet Geological Map, which outlines stratigraphic units composed of alluvial deposits and volcanic materials from Mount Merapi, dating to the Quaternary period. These units include:

• The Old Alluvium Formation (Q_t), consisting of conglomerate, sandstone, silt, and clay, formed during the Plistocene Period.

• The Younger Merapi Volcano Formation (Q_vm), which lies unconformably over the Old Alluvium Formation and Older Merapi volcanic Rocks in most of the western basin, and consists of volcanic breccia, lava, and tuff.

• Alluvium Formation (Q_a), which comprises loose materials such as clay, silt, sand and gravel.

The distribution of aquifers underlying Surakarta City is depicted in two cross-sections (Fig. 2): One running west to east as hydrostratigraphic unit profile A-A', and another running north to south as profile B-B'. Profile A-A', extending from Karangasem to Jebres, shows the presence of a near-surface confined aquifer around Karangasem, Mojosongo, and Jurug, with a thickness ranging from 5 m to 30 m. In the Manahan and Jebres areas, unconfined aquifers with thicknesses between 75 m and 90 m are observed, along with depressed aquifers at greater depths, ranging from 5 m to 10 m.

Profile B-B', which spans from Kadipiro to Kampung Sewu, reveals the presence of unconfined aquifers in Mojosongo, Jebres, and Kampung Sewu, with thicknesses ranging from 25 m to 90 m. Additionally, deeper depressed aquifers with thicknesses ranging from 20 m to 60 m are evident (Putranto et al. 2017).

2 Materials and methods

2.1 Data collection

The tools used in this research include RStudio for processing secondary data, QGIS 2.18 with the MOLUSCE plugin, ArcSWAT for estimating groundwater recharge, and ArcGIS 10.8 for spatial data processing and presentation. The data used in this study include physical characteristics, weather parameters, flow discharge, and demographic observation data, as detailed in Table 1. The overall methodological framework for this study, including the integration of LULC and climate change projections with the SWAT model for groundwater recharge estimation, is presented in Fig. 3.

2.2 Research methodology

2.2.1 Land Use and Land Cover (LULC) classification

Prior to the LULC classification process, satellite imagery was extracted and reprojected using the UTM WGS 1984 Zone 49S coordinate projection. The classification was conducted using a supervised classification approach. The identification of areas categorized as built-up areas, vegetation, water bodies and roads was based on visual identification of land cover, as summarized in Table 2.

2.2.2 LULC change projection

LULC change projections were analyzed using the Modules for Land-Use Change Simulation (MOLUSCE) plugin in QGIS to estimate spatiotemporal changes and calculate LULC transitions between the study intervals (2000–2010, 2010–2020 and 2000–2020). Changes in area and transition probability matrices were derived using LULC data and driving factors, which included landscape classifications from the beginning and end of each interval. A multilayer perceptron approach was employed to model the transition potential, while the driving factors included DEM, slope, distance to rivers, and distance to roads. These variables are commonly used in LULC change analyses, as they provide consistent and physically meaningful data representing both natural and anthropogenic influences on LULC dynamics (Muhammad et al. 2022).

Model validation was performed using observed LULC data from 2020, obtained from image analysis, and compared to the LULC simulated for 2020. From these two datasets, validation accuracy was calculated. Transition potential modeling was conducted using the ANN method, with a neighborhood value of 1 pixel, a learning rate of 0.01, a maximum of 100 iterations, 10 hidden layers, and a momentum of 0.05. For model validation and prediction uncertainty assessment, the MOLUSCE plugin uses the Kappa validation technique (Landis and Koch, 1977), which compares the actual LULC with the projected results.

2.2.3 Multi Model Average (MMA)

The Multi Model Average (MMA) method is a crucial approach in climate projection, aimed at minimizing uncertainty by incorporating outputs from multiple climate models. In this method, several Global Climate Models (GCMs) are run under the same emission scenarios, and their outputs are then averaged. This approach enables researchers to obtain more comprehensive and reliable climate projections, as each individual model possesses unique strengths and limitations.

By combining outputs from various models, the MMA approach helps reduce individual model biases, providing a more robust estimate of future climate changes. This method is particularly important for producing credible projections that inform climate change adaptation and mitigation planning (Tegegne et al. 2020).

2.2.4 Climate parameter projections

Before using climate parameter projection to assess groundwater infiltration with the SWAT model, bias corrections must be applied to ensure comparability between model simulations and actual climate observations. The bias correction approach used in this study is the statistical bias correction method developed by Piani et al. (2010), originally used to correct daily rainfall data from a Multi Model Average (MMA) of climate models against daily rainfall observations in Europe. The same method is adopted here to correct rainfall data from an MMA, specially the average of the five models (CMIP6), using CHIRPS-2.0 observational data. For climate projections in this research, the Shared Socioeconomic Pathways (SSP) SSP2-45 and SSP5-85 scenarios were selected from the Regional Circulation Model (RCM) MMA output.

The SSP2-45 scenario is considered the 'middle' or business-as-usual scenario. It assumes technological and social progress is ongoing, with global GHG emissions peaking around 2040 and declining thereafter. In contrast, the SSP5-85 scenario is considered the high emissions or 'worst-case' scenario, representing a world characterized by high population growth, persistent fossil fuel dependence, and a lag in green technological development, making it suitable for modeling extreme climate change impacts. MMA data was resampled to a 0.05° × 0.05° spatial resolution to describe rainfall and temperature conditions for both the baseline period (1991–2021) and future periods under the SSP2-45 and SSP5-85 scenarios: Namely period 1 (2021–2030) and period 2 (2031–2040). Daily rainfall was bias-corrected using the Quantile Mapping method, following the formula (Piani et al. 2010) (Eq. 1):

(1)$ {P}_{corr,m}={ecdf}_{obs,m}^{-1}\left({ecdf}_{raw,m}\left({P}_{raw,m}\right)\right) $

Where: ecdf = empirical cumulative distribution function, obs,m = the observation value, raw,m = raw model, and P_raw,m = the modelled rainfall. Fig. 4 presents the bias correction procedure applied to global climate model rainfall outputs. The process involves extracting model data at each station, comparing it with observed rainfall, and applying correction factors to reduce systematic errors. The corrected outputs provide improved forecasts that are more consistent with local rainfall observations

The accuracy of rainfall data from the CMIP6 model under historical conditions was first evaluated both before and after applying bias correction. The evaluation metrics used included correlation coefficient (Equation 2) and Mean Absolute Error (MAE) (Equation 3), which were employed to assess the performance of the models at a regional scale.

(2)$ {r}_{xy}=\frac{\displaystyle\sum {(x}_{i}-\bar{x}){(y}_{i}-\bar{y})}{\sqrt{\displaystyle\sum {{(x}_{i}-\bar{x})}^{2}{{(y}_{i}-\bar{y})}^{2}}} $

Model suitability was also evaluated based on the magnitude of error, with Mean Absolute Error (MAE) used to quantify this. If Y_t is the observed value at time t, F_t is the modelled value at the same time, and n is the number of observations, the error is defined as e_t = Y_t – F_t and can be calculated as follows (Willmot and Matsuura, 2005). A smaller MAE value indicates that the modelled results are closer to the observed values.

(3)$ MAE=\frac{1}{n}{\sum }_{t=1}^{n}\left|{e}_{t}\right| $

2.2.5 Calculation of groundwater recharge using the SWAT model

SWAT simulations were conducted by delineating watershed and sub-watershed boundaries using DEM data, defining Hydrologic Response Units (HRUs), processing weather generator and climate input data, and configuring and running the SWAT model. The unsaturated zone plays a critical role in groundwater recharge modeling. In this study, the characteristics of the unsaturated zone were incorporated into the SWAT framework by considering key parameters such as porosity, hydraulic conductivity, and infiltration patterns. The model simulates water movement through the unsaturated zone, accounting for both Land Use/Land Cover changes and climate change projections specific to Surakarta City.

The main hydrological processes modelled in SWAT include surface runoff, infiltration into the soil and root zone, evapotranspiration, soil and snowmelt evaporation, and baseflow (Arnold et al. 1998). The water balance for each HRU is estimated using the following equation (Equation 4) given by (Neitsch et al. 2011):

(4)$ {SW}_{t}={SW}_{0}+\sum _{i=1}^{t}\left({R}_{day}-{Q}_{surf}-{E}_{a}-{w}_{seep}-{Q}_{gw}\right) $

Where: SW_t is the final groundwater content (mm); SW₀ is the initial groundwater content on day i (mm); R_day is the rainfall on day i (mm); Q_surf is the surface runoff on day i (mm); E_a is the evaporation on day i (mm); W_seep is the percolation into the vadose zone in the soil profile on day i (mm); and Q_gw is the return flow from groundwater to the stream on day i (mm).

SWAT simulates both unconfined and confined aquifers in each sub-basin. The infiltration of water form the vadose zone into the aquifer on day i is calculated using Equation 5:

(5)$ {W}_{rchrg,i} = \left(1 - \mathrm{exp}\left[ -\frac{1}{{\delta }_{gw}}\right]\right)\cdot{w}_{seep} + \mathrm{exp}\left[ -\frac{1}{{\delta }_{gw}}\right]\cdot{w}_{rchrg,i-1} $

Where: δ_gw is the groundwater delay time in days and w_seep is the total water percolating from the bottom of the soil profile in mm (Neitsch et al. 2011). This hydrological component represents one of several results generated from water movement simulations in SWAT.

2.2.6 Calibration with SWAT-CUP

The calibration of the SWAT model involves the use of various observed and derived input datasets. Selecting appropriate parameters significantly improves calibration efficiency and accuracy, reduces projection uncertainty (Arnold et al. 2012a), and ensures a more realistic representation of the hydrological processes (Meaurio et al. 2015). To facilitate this, the SWAT Calibration and Uncertainty Program (SWAT-CUP) was developed to perform calibration, validation, and sensitivity analysis of SWAT models (Arnold et al. 2012b). SWAT-CUP supports four calibration algorithms: namely GLUE, ParaSol, MCMC, and SUFI-2. In this study, model calibration was conducted using the SWAT-CUP software, with discharge data from the Jurug Water Post in Surakarta City, covering the period from 2011 to 2020. This calibration was essential for evaluating the agreement between simulated and observed streamflow values. Model accuracy was assessed through graphical comparison of observed and simulated discharge along with the calculation of the following statistical performance metrics: Nash-Sutcliffe Efficiency (NSE), Percent Bias (PBIAS), and the coefficient of determination (R²) (Equations (6-8)).

(6)$ NSE= 1-\frac{{\displaystyle\sum }_{t=1}^{T}{\left({Q}_{m}^{t}-{Q}_{0}^{t}\right)}^{2}}{{\displaystyle\sum }_{t=1}^{T}{\left({Q}_{0}^{t}-{\bar{Q}}_{0}\right)}^{2}} $

(7)$ PBIAS= 100\frac{{\displaystyle\sum }_{t=1}^{T}{Q}_{m}^{t}-{Q}_{0}^{t}}{{\displaystyle\sum }_{t=1}^{T}{Q}_{0}^{t}} $

(8)$ {R}^{2}={\left[\frac{{\displaystyle\sum }_{t=1}^{T}\left({Q}_{0}^{t}-{\bar{Q}}_{0}\right)\left({Q}_{m}^{t}-{\bar{Q}}_{m}\right)}{\sqrt{{\displaystyle\sum }_{t=1}^{T}{\left({Q}_{0}^{t}-{\bar{Q}}_{0}\right)}^{2}}\sqrt{{\displaystyle\sum }_{t=1}^{T}{\left({Q}_{m}^{t}-{\bar{Q}}_{m}\right)}^{2}}}\right]}^{2} $

For the model performance to be considered satisfactory, NSE should be greater than 0.5, PBIAS should fall within ± 25%, and R² should be greater than or equal to 0.70 (Moriasi et al. 2007; Bucton et al. 2022). Following successful calibration and validation of the hydrological model, the impact assessment of climate and land use changes on groundwater recharge in Surakarta City was conducted under multiple projection scenarios.

3 Results and discussion

3.1 LULC change from 2000 to 2020

In general, the dominant land use types in Surakarta City are built-up areas, including industrial zones, service areas, and residential settlements. As shown in Table 3, the built-up area has experienced a sustantial expansion, accounting for 69.65% of the total area by 2020. In contrast, vegetation cover has shown a notable decline, decreasing from 9.55 km² in 2000 to 7.28 km² in 2020 (Fig. 5). Additionally, road infrastructure expanded by 1.82% over the 20-year period, reflecting ongoing urban development. Water bodies also exhibited a slight increase of 0.32% during this period.

3.2 Validation of LULC projection model

The comparison between the modeled and observed land use conditions in 2020 reveals only minor differences, as detailed in Table 4. The projected built-up area differs by 1.26% from the observed value, while vegetation shows a 0.75% difference. Water bodies and road areas exhibit even smaller percentage deviations. Following the generation of the LULC for 2020 based on 2010 data, the model's accuracy was assessed using the kappa coefficient. The resulting kappa a value of 0.82 (Table 4) indicates a strong agreement between the projected and observed land cover (Landis and Koch, 1977). This validation outcome suggests that the model has a good fit, with projection uncertainty remining below 20%. Therefore, the model is considered reliable for forecasting future land use changes using the 2010 and 2020 LULC datasets.

Fig. 6 illustrate the LULC prediction results for 2030 and 2040, generated using the ANN method, indicate that built-up areas will continue to dominate Surakarta City, expanding to 33.21 km² in 2030 and 33.56 km² in 2040. Vegetation cover is projected to increase slightly by 2030, while road areas are expected to experience a slight decrease of 0.24%. The continuous growth of built-up areas, particularly near roads, align with findings by Nabila (2023), which highlight that proximity to roads accelerates land use change.

The spatial dynamics of LULC from 2000 to 2040 show changes across all categories except water bodies, which only fluctuated around 2010. Social and physical factors remain key drivers of land use change, with built-up areas concentrating along roads and vegetation expanding in suburban zones. These trends corroborate earlier studies by Fitriana et al. (2017), who projected increased built-up and vegetation areas in Surakarta City between 2030 and 2040. Similarly, research by Wibisono et al. (2023) noted that residential development during this period is progressing southward towards Sukoharjo Regency and westward. Fig. 7 illustrate the projected LULC changes in the study area. The results show a continuous increase in built-up areas towards 2030 and 2040, accompanied by a gradual decline in vegetation and other land cover types, reflecting the trend of urban expansion.

To assess the conformity of these predictions with the regional spatial planning (RTRW) of Surakarta City, the predicted LULC results were compared with land allocation outlined in the RTRW (Tables 5 and 6). In 2030, approximately 76.43% (35.71 km²) of the predicted LULC was consistent with the RTRW, while 23.57% (11.01 km²) was not. By 2040, the conformity slightly increased to 77.39% (36.16 km²), and the non-conforming areas decreased to 22.61% (10.56 km²). These discrepancies indicate areas planned in the RTRW that are projected differently by the model, primarily involving the expansion of built-up areas. This highlights the need to evaluate and adjust spatial planning policies to ensure sustainable development, particularly in managing land-use changes that directly impact the environment and groundwater recharge capacity. More adaptive and responsive spatial planning policies are essential to address ongoing urbanization, reduce land-use conflicts, and mitigate adverse impacts on groundwater resources.

3.3 Projected future climate of Surakarta City

Climate projections were derived from an ensemble model, with spatial resolution enhanced to 0.05° × 0.05°, using the Shared Socioeconomic Pathways (SSP) climate scenario. The SSP scenarios employed in this study are SSP2-45 and SSP5-85, representing future climate conditions in Surakarta City. The bias correction method applied to rainfall data uses the gamma distribution, while other climate parameters, such as temperature, solar radiation, relative humidity, and wind speed, are corrected using the normal distribution. The choice of statistical distribution for bias correction is based on the distinct characteristics of each climate variable. The gamma distribution is well suited for rainfall data because it effectively capture its asymmetric, discontinuity, and long right-skewed tail, which reflect the complex and irregular rainfall patterns typically observed in tropical regions (Stull, 2017). Conversely, the normal distribution is appropriate for variables like temperature, solar radiation, wind speed, and relative humidity due to their generally symmetric and continuous nature, with values fluctuate regularly around the mean. These characteristics are also commonly observed in tropical climates, which are characterised by relatively stable temperature and climate patterns. This methodological choice is supported by empirical statistical analyses in meteorology and climatology, which consider the probabilistic structure, continuity, and distributional characteristics of each climate variable.

Fig. 8 presents rainfall data for historical (1991–2020) and projected (2021–2040) periods under the SSP2-45 and 5–85 scenarios. Fig. 8(a) shows the rainfall data before bias correction, based on records from the Pabelan climatological station. Prior to correction, the Multi Model Average (MMA) overestimated rainfall compared to the observed historical data, whereas the scenario-based model data tended to underestimate rainfall. In contrast, Fig 8(b) shows the rainfall data after bias correction. It is evident that both the MMA data and scenario models more accurately capture the historical rainfall patterns in Surakarta City. The projected rainfall for 2021–2040 indicates an overall increase compared to the historical period, with the SSP2-45 scenario yielding higher rainfall values than SSP5-85. To evaluate model performance, statistical analysis were conducted using the correlation coefficient (R) and Mean Absolute Error (MAE). A higher absolute R value indicates a stronger linear relationship between observed and modelled data, suggesting that the model effectively captures the overall rainfall trend. A lower MAE value reflects improved model accuracy in forecasting monthly rainfall patterns (Table 6).

3.4 SWAT model calibration and validation

The delineation of the Bengawan Solo watershed in Surakarta City using the SWAT model was performed automatically. This process utilized river network maps, DEM data, and watershed boundary delineations. The SWAT delineation resulted in 19 sub-watersheds covering a total area of 537,28 km² (Fig. 9). The SWAT model calibration was conducted by comparing simulated monthly discharge data (flow out from the rch file) with observed monthly discharge data from 2011–2020 at the Jurug post-Surakarta City. Calibration involved adjusting parameters related to river flow (*.bsn), base flow (*.gw), main channel (*.rte), land management (*.mgt), and soil properties (*.sol). Table 7 presents the parameters used during calibration along with their final calibrated values. The parameter codes 'R' and 'V' for each parameter indicate their respective characteristics. 'R' refers to calibration factors adjustable in the SWAT model, while 'V' denotes parameters associated with groundwater flow (baseflow). The performance the SWAT model in simulating monthly discharge for the Bengawan Solo River is quantitatively evaluated using statistical indicators summarized in Table 8.

The calibration results show an R² value greater than 0.50 (Table 9), indicating that the simulated discharge is acceptable (Santhi et al. 2001), while the KGE value above 0.70 falls into the good performance category (Pechlivanidis and Arheimer, 2015). The parameters that most strongly influence the simulated discharge to match the observed values are SOL_K, CH_N1, CH_K2, and CN2. These parameters significantly impact the simulated discharge based on land use and soil type, particularly SOL_K, which represents saturated hydraulic conductivity. An increase in SOL_K leads to greater percolation, thereby increasing base flow values (Widyastuti et al. 2018). Fig. 10 illustrates the comparison of monthly discharge data before and after calibration for the period 2011-2020. The calibration model demonstrates a good agreement between observed and simulated discharge. This can be seen in the discharge time series from 2011 to 2020, where the simulated discharge fluctuations closely follow the observed data pattern. Fig. 11 presents a comparison between observed discharge and the calibrated model results over the 2011–2020 period.

The calibration results demonstrate that the SWAT model can effectively simulate the overall maximum discharge events. The model accurately capture extreme discharge values, such as those occurring in the 3rd and 67_th months. For instance, the simulated discharge reaches a value of 174.7 m³/sec, closely matching the observed discharge of 185.7 m³/sec during these peak events.

3.5 Impact of LULC and climate change on groundwater recharge

This study demonstrates that both climate change and LULC change significantly contribute to the decline of groundwater recharge in the study area. Hydrological model analysis reveals that climate-driven shifts in rainfall patterns have a greater impact on the variability of annual water infiltration. Rising temperature lead to increased evaporation, while more intense and concentrated rainfall events reduce infiltration efficiency during the rainy season. In addition, the conversion of agricultural land to built-up areas exacerbates the decrease in infiltration due to the expansion of impervious surfaces that directly limit water percolation. Sensitivity analysis indicates that although both factors influence groundwater infiltration, LULC changes exert a more substantial long term impact because they are permanent, whereas climate change effects fluctuate annually depending on rainfall variability.

The quantity of groundwater infiltration was calculated based on the results of LULC projections and climate projections using two climate scenarios: SSP2-45 and SSP5-85. As shown in Figs.12 and 13, the average monthly infiltration in Surakarta City has decreased across all periods: historical, current and projected, with values of 86.04 mm/month, 85.91 mm/month, and 84.98 mm/month, respectively, over a 10-year span. The projected infiltration results indicate a significant decline under both climate scenarios, with infiltration under SSP5-85 consistently lower than under SSP2-45. This difference is attributed to the lower annual average precipitation in SSP5-85 (2,076.51 mm/year) compared to SSP2-45 (2,149.74 mm/year). Although the built-up area shows a 1% decrease in 2030, the overall trend from 2000 to 2040 is an increase in built-up areas, which reduces the surface available for groundwater infiltration (Suharyanto, 2018). Additionally, the decline in infiltration is further aggravated by reductions in soil physical quality, such as texture, structure, and porosity, caused by development activities, which diminish the soil's infiltration capacity even in areas still covered by green vegetation (Lei et al. 2019). Tables 10 and 11 present detailed annual groundwater recharge projections considering both LULC and climate change; the numbers 245 and 585 in the legend refer to the SSP2-45 and SSP5-85 rainfall projection scenarios, respectively.

Fig. 14 shows areas that consistently experienced declining infiltration values from 2000 to 2020, particularly in sub-basin 16. This area is predominantly covered by forest and shrub landcover types and features relatively steep slope at the foot of Mount Lawu. In flat areas such as Surakarta City, low groundwater infiltration estimates can be attributed to clay soils, which have poor infiltration capacity (Saputra et al. 2020). This soil type tends to hinder water from seeping into the ground. In addition, reduced infiltration is also linked to high evapotranspiration rates (Luetkemeier et al. 2022), where water evaporates from soil and plant surfaces, decreasing the amount available for groundwater recharge. Increasing infiltration and managing runoff effectively are essential strategies for improving water resource sustainability, tailored to the specific characteristics of each watershed. To mitigate the decline in groundwater recharge, appropriate water management strategies, such as water use regulation, reforestation, and land conservation practices, must be implemented to enhance water availability in the region.

Land use change has a direct impact on groundwater recharge and dynamics (Scanlon et al. 2005). Additionally, the lack of effective land use policies contributes significantly to groundwater degradation in urban areas (Sajjad et al. 2022). Fig. 15 illustrates the spatiotemporal pattern of groundwater recharge in Surakarta City. As previously discussed, infiltration is highly dependent on both rainfall and LULC conditions. Surakarta City is located within sub-basins 1–6, highlighted in the red box, while the other sub-basins are excluded from the analysis, as this study focuses solely on Surakarta City. As shown in Fig. 13, the projected conditions for both the 2021–2030 and 2031–2040 periods, under both climate scenarios, appear visually similar on the maps and indicate a continuing decline in groundwater infiltration conditions compared to previous years.

This finding aligns with research by Alfandhani et al. (2021), which attributes critical infiltration conditions to the density of built-up areas and the scarcity of vegetation. The type of land cover significantly influences infiltration rates, dense vegetation promotes greater infiltration, while sparse cover results in lower rates (Abbas et al. 2021). Furthermore, soil type plays a crucial role in determining the effectiveness of the water catchment areas (Pandiangan et al. 2021). In Surakarta, the area is predominantly composed of grumusol soil, which is known for its low infiltration capacity. In contrast, soils with a higher proportion of sand generally exhibit better infiltration and water absorption potential (Sejati, 2020).

Geologically, Surakarta City is located above a complex aquifer system consisting of Quaternary alluvial deposits in the upper layer, underlain by volcanic materials and Tertiary sedimentary rocks in the deeper layer. The expansion of built-up area has a direct impact on the uppermost Quaternary alluvial deposits. These deposits, consisting of clay, silt, sand and gravel, have high porosity and permeability, making them highly conducive to groundwater infiltration (Putranto and Kusuma, 2009). However, urban development has led to the extensive coverage of these permeable layers with impermeable surfaces such as asphalt and concrete. As a result, the natural capacity of the alluvial layer to absorb and transmit rainwater into the underlying aquifer system is significantly diminished.

Furthermore, the reduction in surface infiltration adversely affects the recharge of the volcanic layer beneath. This volcanic layer, composed of tuff, breccia and lahar sediments, is crucial for conducting water to deeper aquifers. With less water percolating from the surface, the effectiveness of these volcanic layers in transmitting water to deeper aquifers is also reduced.

4 Conclusions

This study evaluated the impacts of climate change and Land Use/Land Cover (LULC) changes on groundwater recharge in Surakarta City using an integrated modeling framework comprising the SWAT hydrological model, the CA-ANN based LULC projection model, and an ensemble of five RCM models under two SSP scenarios.

The results indicate that future rainfall is projected to increase under both scenarios, with the SSP2-45 scenario showing higher rainfall values compared to SSP5-85.

The CA-ANN model, implemented through the MOLUSCE plugin in QGIS, was employed to project Land Use and Land Cover (LULC) changes using various driving factors, including physical, accessibility, and socioeconomic data. The model demonstrated good simulation performance, with a Kappa accuracy value exceeding 0.61, indicating reliable predictive capability. According to the CA-ANN projections, Surakarta City will experience a 4.15% and 4.9% increase in built-up area and a 7.32% and 8.58% reduction in green space by 2030–2040, relative to the year 2000.

The SWAT model was used to assess groundwater recharge under historical, current and projected conditions based on LULC and climate change scenarios. The projection indicate a consistent decline in groundwater recharge across both timeframes (2021–2030 and 2031–2040), with more pronounced decreases during the rainy season, while slight increases may occur during the dry season.

Overall, the findings confirm that both climate change and LULC transformation significantly affect groundwater recharge in Surakarta City. The average monthly groundwater recharge has gradually declined from 86.04 mm/month in 2000, to 85.91 mm/month in 2010, and 84.98 mm/month in 2020. Under future scenarios, recharge is expected to further decrease to 83.85 mm/month (SSP2-45) and 78.25 mm/month (SSP5-85) by 2030, and to 82.10 mm/month (SSP2-45) and 77.44 mm/month (SSP5-85) by 2040. This downward trend is closely linked to urban expansion and the loss of vegetated areas, which reduce infiltration and increase surface runoff.

Given Surakarta City's limited spatial capacity and growing population, it is particularly vulnerable to these impacts. Therefore, the implementation of effective and stringent urban planning policies is essential to mitigate groundwater recharge decline and ensure long-term water sustainability.

References

[1] Abbas Z, Yang G, Zhong Y, et al. 2021. Spatiotemporal change analysis and future scenario of lulc using the CA-ANN approach: A case study of the greater bay area, China. Land, 10(6): 584. DOI: 10.3390/land10060584

[2] Adhikari RK, Mohanasundaram S, Shrestha, S. 2020. Impacts of land-use changes on the groundwater recharge in the Ho Chi Minh city, Vietnam. Environmental Research, 185: 109440. DOI: 10.1016/j.envres.2020.109440

[3] Alfandhani RS, Hizbaron DR, Widyastuti M. 2021. Kajian Pengaruh Kondisi Resapan Air pada Pola Pemanfaatan Ruang di Sub DAS Jlantah-Walikun pada Wilayah DAS Bengawan Solo Hulu. Jurnal Pendidikan Tambusai, 5(3): 11236–11244. (in Indonesian

[4] Arnold JG, Kiniry R, Srinivasan R, et al. 2012. Input/output documentation soil & water assessment tool.Texas Water Resources Institute: Thrall, TX, USA, 1-650.

[5] As-syakur AR, Suarna IW, Adnyana IWS, et al. 2008. Studi Perubahan Penggunaan Lahan di DAS Badung. Bumi Lestari, 10(2): 200–207. (in Indonesian

[6] Boretti A, Rosa L. 2019. Reassessing the projections of the world water development report. NPJ Clean Water, 2(1): 15. DOI: 10.1038/s41545-019-0039-9.

[7] Bucton BGB, Shrestha S, Mohanasundaram S, et al. 2022. Impacts of climate and land use change on groundwater recharge under shared socioeconomic pathways: A case of Siem Reap, Cambodia. Environmental Research, 211: 113070. DOI: 10.1016/j.envres.2022.113070.

[8] Desky AF. 2022. Sosiologi pedesaan dan perkotaan. (In Indonesian

[9] DLH Kota Surakarta. 2021. Indeks Kualitas Lingkungan Hidup Kota Surakarta Tahun 2021. (In Indonesian

[10] DLH Kota Surakarta. 2022. Indeks Kualitas Lingkungan Hidup Kota Surakarta Tahun 2022. (In Indonesian

[11] Ferreira RS, Uagoda RE. 2017. Revista Brasileira de Geografia Física Análise da predição do balanço hídrico da bacia do ribeirão do Gama-DF através do modelo SWAT. In Revista Brasileira de Geografia Física, 10(03): 880-893. (in Portuguese

[12] Fitriana AL, Subiyanto S, Firdaus HS. 2017. Model Cellular Automata Markov untuk Prediksi Perkembangan Fisik Wilayah Permukiman Kota Surakarta menggunakan Sistem Informasi Geografis. Jurnal Geodesi Undip Oktober, 6(4): 246-253. (In Indonesian

[13] Guo H, Hu Q, Jiang T. 2008. Annual and seasonal streamflow responses to climate and land-cover changes in the Poyang Lake basin, China. Journal of Hydrology, 355(1-4): 106−122. DOI: 10.1016/j.jhydrol.2008.03.020.

[14] Guo XD, Liu Q, Li WP, et al. 2024. Groundwater imbalance and its mutual feedback relationship with land use in West Liaohe Plain. Hydrogeology & Engineering Geology, 51(4): 77−87. DOI: 10.16030/j.cnki.issn.1000-3665.202310001.

[15] Halefom A, Sisay E, Khare D, et al. 2017. Hydrological modeling of urban catchment using semi-distributed model. Modeling Earth Systems and Environment, 3(2): 683−692. DOI: 10.1007/s40808-017-0327-7.

[16] Indhanu N, Chalermyanont T, Chub-Uppakarn T. 2025. Spatial assessment of land use and land cover change impacts on groundwater recharge and groundwater level: A case study of the Hat Yai basin. Journal of Hydrology: Regional Studies, 57: 102097. DOI: 10.1016/j.ejrh.2024.102097.

[17] Landis JR, Koch GG. 1977. The measurement of observer agreement for categorical data. Biometrics, 159-174.

[18] Lei Z, Yang J, Yang Y, et al. 2019. Impact of urban sprawl on dry island effects and its countermeasures. Sustainable Cities and Society, 51.

[19] Luetkemeier R, Söller L, Frick-Trzebitzky F. 2022. Anthropogenic pressures on groundwater. In Encyclopedia of Inland Waters, Second Edition (3): 548-559. DOI: 10.1016/B978-0-12-819166-8.00183-3

[20] Meaurio M, Zabaleta A, Uriarte JA, et al. 2015. Evaluation of SWAT models performance to simulate streamflow spatial origin: The case of a small forested watershed. Journal of Hydrology, 525: 326−334. DOI: 10.1016/j.jhydrol.2015.03.050.

[21] Mengistu TD, Chung IM, Kim MG, et al. 2022. Impacts and implications of land use land cover dynamics on groundwater recharge and surface runoff in east african watershed. Water (Switzerland), 14(13): 2068. DOI: 10.3390/w14132068.

[22] Moriasi DN, Arnold JG, Van Liew MW, et al. 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Transactions of the ASABE, 50(3): 885−900. DOI: 10.13031/2013.23153.

[23] Muhammad R, Zhang W, Abbas Z, et al. 2022. Spatiotemporal change analysis and prediction of future land use and land cover changes using QGIS MOLUSCE Plugin and remote sensing big data: A case study of Linyi, China. Land, 11(3): 419. DOI: 10.3390/land11030419.

[24] Mustafa A, Van Rompaey A, Cools M, et al. 2018. Addressing the determinants of built-up expansion and densification processes at the regional scale. Urban Studies, 55(15): 3279−3298. DOI: 10.1177/0042098017749176.

[25] Nabila DA. 2023. Pemodelan prediksi dan kesesuaian perubahan penggunaan lahan menggunakan Cellular Automata-Artificial Neural Network (CA-ANN). Tunas Agraria, 6(1): 41−55. DOI: 10.31292/jta.v6i1.203 (In Indonesian

[26] Neitsch SL, Arnold JG, Kiniry JR, et al. 2011. Soil and water assessment tool theoretical documentation version 2009.

[27] Pandiangan NL, Diara IW, Kusmiyarti TB. 2021. Analisis Kondisi Daerah Resapan Air Kecamatan Sukasada Kabupaten Buleleng menggunakan Sistem Informasi Geografis. Jurnal Agroekoteknologi Tropika, 10(3): 324–336. (In Indonesian

[28] Pechlivanidis LG, Arheimer B. 2015. Large scale hudrological modelling by using modified PUB recommendations: The India-HYPE case. Hydrology and Earth System Sciences, 19: 4559−4579. DOI: 10.5194/hess-19-4559-2015.

[29] Piani C, Haerter JO, Coppola E. 2010. Statistical bias correction for daily precipitation in regional climate models over Europe. Theoretical and Applied Climatology, 99(1-2): 187−192. DOI: 10.1007/s00704-009-0134-9.

[30] Piyumi MMM, Abenayake C, Jayasinghe A, et al. 2021. Urban flood modeling application: Assess the effectiveness of building regulation in coping with urban flooding under precipitation uncertainty. Sustainable Cities and Society, 75: 103294. DOI: 10.1016/j.scs.2021.103294.

[31] Putranto TT, Krisna W, Revina R, et al. 2016. Groundwater vulnerability to contamination assessment using GOD method in Surakarta city, central Java. Proceeding PIT-PAAI.

[32] Putranto TT, Kusuma KI. 2009. Permasalahan airtanah pada daerah urban. Teknik, 30(1): 48−56.

[33] Putranto TT, Widiarso DA, Udin MI. 2017. Zonasi Potensi Air Tanah Kota Surakarta, JawaTengah. Proceeding, Seminar Nasional Kebumian. (In Indonesian

[34] Ridho FT, Nurhidayati E, Damayanti HN, et al. 2020. Komparasi Model Cellular Automata dalam Memprediksi Perubahan Lahan Sawah di Kabupaten Purworejo. Jukung Jurnal Teknik Lingkungan, 6(2): 157−167. DOI: 10.20527/jukung.v6i2.9259 (In Indonesian

[35] Sajjad MM, Wang J, Abbas H, et al. 2022. Impact of climate and land-use change on groundwater resources, study of Faisalabad District, Pakistan. Atmosphere, 13(7): 1097. DOI: 10.3390/atmos13071097.

[36] Santhi C, Allen PM, Muttiah RS, et al. 2008. Regional estimation of base flow for the conterminous United States by hydrologic landscape regions. Journal of Hydrology, 351(1-2): 139−153. DOI: 10.1016/j.jhydrol.2007.12.018.

[37] Santhi C, Arnold JG, Williams JR, et al. 2001. Validation of the SWAT Model on a Large Rwer Basin with Poin and Nonpoint Sources. Journal of the American Water Resources Association, 37: 1169−1188. DOI: 10.1111/j.1752-1688.2001.tb03630.x.

[38] Saputra DR, Yudono ARA, Partoyo. 2020. Assessment of the groundwater recharge potential areas using GIS in Kajor Kulon Hamlet, Selopamioro, Imogiri, Bantul, Yogyakarta. Jurnal Geografi Lingkungan Tropik, 4(2): 4. DOI: 10.7454/jglitrop.v4i2.89

[39] Scanlon BR, Reedy RC, Stonestrom DA, et al. 2005. Impact of land use and land cover change on groundwater recharge and quality in the southwestern US. Global Change Biology, 11(10): 1577−1593. DOI: 10.1111/j.1365-2486.2005.01026.x.

[40] Sejati AE. 2020. Analisis Perubahan Tutupan Lahan dan Pengaruhnya terhadap Kemampuan Resapan Air di Daerah Aliran Sungai. Jurnal Planologi, 17(1): 1−12. DOI: 10.30659/jpsa.v17i1.7555. (In Indonesian

[41] Siddik MS, Tulip SS, Rahman A, et al. 2022. The impact of land use and land cover change on groundwater recharge in northwestern Bangladesh. Journal of Environmental Management, 315: 115130. DOI: 10.1016/j.jenvman.2022.115130.

[42] Sisay E, Halefom A, Khare D, et al. 2017. Hydrological modelling of ungauged urban watershed using SWAT model. Modeling Earth Systems and Environment, 3(2): 693−702. DOI: 10.1007/s40808-017-0328-6.

[43] Stull R. 2017. Practical meteorology an algebra-based survey of atmospheric science. Available on http://www.eos.ubc.ca/~geol/softeng/practical.html

[44] Suharyanto. 2018. Analisis Dampak Perubahan Tata Guna Lahan terhadap Daerah Resapan Air di Wilayah Perkotaan. Jurnal Teknologi Lingkungan, 19(2): 149−158. (In Indonesian

[45] Tamm O, Kokkonen T, Warsta L, et al. 2023. Modelling urban stormwater management changes using SWMM and convection-permitting climate simulations in cold areas. Journal of Hydrology, 622: 129656. DOI: 10.1016/j.jhydrol.2023.129656.

[46] Tegegne G, Melesse AM, Worqlul AW. 2020. Development of multi-model ensemble approach for enhanced assessment of impacts of climate change on climate extremes. Science of the Total Environment, 704: 135357. DOI: 10.1016/j.scitotenv.2019.135357.

[47] Wibisono P, Miladan N, Pamardhi-Utomo R. 2023. Hubungan Perubahan Kerapatan Vegetasi dan Bangunan terhadap Suhu Permukaan Lahan: Studi Kasus di Aglomerasi Perkotaan Surakarta, 5: 148–162.

[48] Widyastuti MT, Taufik M, Santikayasa IP. 2018. Prediksi Debit Jangka Panjang untuk Sungai Bengawan Solo. Jurnal Geografi, 15(2): 71−82. DOI:10.15294/jg.v15i2.15387 (In Indonesian

[49] Willmot CJ, Matsuura K. 2005. Advantages of the Mean Absolute Error (MAE) over the Root Mean Square Error (RMSE) in assessing average model performance. Climate Research, 30(1): 79−82. DOI: 10.3354/cr030079.

[50] Xu T, Gao J, Coco G. 2019. Simulation of urban expansion via integrating artificial neural network with Markov chain–cellular automata. International Journal of Geographical Information Science, 33(10): 1960−1983. DOI: 10.1080/13658816.2019.1600701.

[51] Yunus HS. 2005. Manajemen kota: Perspektif spasial. Yogyakarta: Pustaka pelajar, 2005. (In Indonesian

[52] Zhang JK, Wang JP, Shi JS. 2024. Attribution analysis of water-sediment variation under the influence of climate change and human activities in the Kuye River Basin. Hydrogeology & Engineering Geology, 51(6): 47−59. DOI: 10.16030/j.cnki.issn.1000-3665.202312053.

[53] Zhang Q, Liu J, Singh VP, et al. 2016. Evaluation of impacts of climate change and human activities on streamflow in the Poyang Lake basin, China. Hydrological Processes, 30(14): 2562−2576. DOI: 10.1002/hyp.10814.

RIGHTS & PERMISSIONS

Journal of Groundwater Science and Engineering Editorial Office