The study area is the Upper Bengawan Solo Watershed (BSH), the upstream part of the Bengawan Solo system, the largest river basin on Java Island Figure 2. Geographically, this area is located in Central Java and covers the upstream region, supplied by the Merapi–Merbabu volcanic complex to the west and the Lawu to the east. Generally, the area's terrain is dominated by relatively flat plains; however, the northeastern and northwestern sections feature undulating to rolling terrain near the mountain foothills. The hydrological segmentation classifies BSH as one of the three main sub-watersheds within the Bengawan Solo system. The study boundaries are defined by hydrological (watershed) boundaries based on the digital elevation model (DEM) and align with administrative borders ((Absori et al., 2023); (Anna et al., 2023); (Anna & Priyana, 2015); (Santhyami et al., 2025)).

The topographic conditions and river networks primarily influence flood dynamics in BSH. Past studies have shown varying vulnerability to flooding: areas around Surakarta (the urban core) tend to be more susceptible, while certain regions in Wonogiri, Karanganyar, and Boyolali are comparatively less at risk. This pattern aligns with the morphological gradient (flat versus undulating or steep terrain), drainage density, and land use/cover structure along the main river corridor. Variations in vulnerability among districts highlight the need for locally tailored management, especially during the rainy season when increased rainfall strains urban drainage systems (Jumadi et al., 2025). Administratively, BSH includes at least six districts/cities: Boyolali, Surakarta, Karanganyar, Sukoharjo, Klaten, and Wonogiri. The Surakarta urban area acts as a hub for population and economic activity, leading to higher exposure of residents and property in this region (Santhyami et al., 2025).

1. Population Centers (GHS-POP R2023A)

The population center dataset is derived from the latest Global Human Settlement Layer (GHSL) release (GHS global population grid multitemporal, GHS-POP R2023A) (Pesaresi et al., 2024)(Schiavina et al., 2023). The population pixel values from GHSL, typically at 100 m resolution, were converted into point centers using Raster to Point in ArcGIS Pro. Each point retains a center_id, pixel center coordinates, and population count (the number of residents in that pixel). These points are used as references to generate synthetic individuals within a 100 m radius.

2. Attribute Frequency Tables (Age, Gender, Education, Occupation)

Four tables display the proportions of demographic and social categories at the watershed level, which combines four districts or cities, for age, gender, education, and occupation. For age, the labels are given as age ranges; during synthesis, these are converted to specific ages by sampling integers within each range. For gender, education, and occupation, category selection is weighted by frequency. Substantive rules are applied to maintain social realism, such as automatically assigning individuals under 15 years old to the occupation 0 (no occupation).

3. Flood Hazard Data Based on RS–GIS

The hazard index Figure 3 is compiled from physical-environmental parameters derived from open-access spatial data ((Jumadi et al., 2024); (Jumadi et al., 2024)). The core framework uses eleven parameters: elevation (El), slope (Sl), flow accumulation (FA), distance to river (DR), drainage density (DD), topographic wetness index (TWI), curvature (Cu), land cover/use (LULC), NDVI, rainfall (Rf), and soil moisture (SM). The primary data sources include SRTM 30 m (El, Sl, FA, DR, DD, TWI, Cu), Sentinel-2 10 m (LULC, NDVI), GPM v6 (Rf), and SMAP (SM). All parameters are projected to the working CRS, standardized in spatial resolution, and then reclassified to an ordinal scale of 1–5 (very low – very high) with local thematic thresholds. Weights among parameters are explicitly assigned (e.g., El/Sl/FA/DR as major controlling factors), and the weighted overlay results in a hazard surface H. Physically, low or flat locations and those near rivers tend to experience flow accumulation and low flow velocity; high NDVI values indicate vegetative cover that may delay runoff; TWI and SM capture soil saturation; flat or concave curvature leads to pooling (see further explanations in (Jumadi et al., 2024); (Jumadi et al., 2024)).

3. Supporting Data

The research uses administrative boundaries (districts), river networks, watershed boundaries, and ESRI Sentinel 2 Land Cover. Administrative boundaries are essential for reporting, especially when comparing results with different regions. The river network helps interpret the results. ESRI Sentinel-2 land cover is needed to validate the distribution of SSP.

The characteristics and sources of all datasets are summarised in Table 1.

1. Spatial Harmonization

Preprocessing starts with aligning the coordinate system and analysis grid so that all layers function properly without distortion. All data is projected to UTM 49S (EPSG:32749) as a standard CRS to ensure consistency in distance and area measurements. The extent is defined by the mask of the Upper Solo River Basin (BSH), which is then used as a clipping boundary for all raster and vector data.

2. Extraction of GHSL Population Data

The GHSL population raster is converted into points using ArcGIS Pro (Raster to Point). These points are clipped by the watershed mask and assigned a unique ID (center_id). The population attribute (number of people per pixel) is cleaned by removing negative or empty values. Points with a population of zero are discarded to improve computational efficiency. If there is an off by one issue at the watershed boundary where GHSL cells cross it, a small adjustment is applied to the analysis grid. All points are stored in the working CRS and exported to CSV for use in the Python pipeline (SSP Algorithm).

1. Principles and Assumptions

The development of SSP aims to represent the population at the individual level in a synthetic way while maintaining overall consistency with official statistics and spatial accuracy regarding settlement areas. The main principles are: (i) aggregate consistency—matching the number of individuals generated at each center point to the population count from GHS-POP; (ii) micro-spatial realism—dispersing individuals around the pixel center within a fixed radius (100 m) with a denser radial distribution at the center; (iii) attribute consistency—sampling age, gender, education, and occupation based on weighted standardized frequency tables; and (iv) substantive rules—to avoid unrealistic attribute combinations (e.g., children working), such as age < 15 → education = 0, occupation = 0. Other key assumptions include: (a) category frequencies reflect the population makeup at the watershed level; (b) micro-spatial dispersion within a 100 m radius sufficiently represents the intrapixel variability observed in GHS-POP; (c) ages within intervals are assumed uniform unless additional information is available.

2. Spatial Algorithm (Individual Placement)

For each center point (cx, cy) with a population n, generate n individual positions using radial jitter, where the angle theta is uniformly distributed between 0 and 2π, and the radius r is drawn from a Beta(2, 2,5) distribution scaled by RADIUS_M. This creates a denser clustering around the center. Offsets are calculated as dx = r* cos (theta) and dy = r* sin (theta).

3. Demographic-Social Attribute Sampling

Individual attributes are assigned through weighted random sampling based on frequency tables. For age, the ID label can indicate a range (e.g., “a–b”), a single value, or “a+.” The range “a–b” is converted into specific ages by selecting a random integer within [a, b]; labels “a+” are limited to [a, a+4] to avoid an infinite tail; single values are used directly. For gender, education, and occupation, categories are selected based on their relative probability (frequency divided by total). After age is assigned, the rule age < 15 is applied to disable education and occupation (set to “0”).

4. Computational Implementation

The implementation uses Python on Google Colab Pro (High-RAM). It reads population centers from DBF files with dbfread and converts data types to numeric. The source code for SSP generation is available at https://doi.org/10.17605/OSF.IO/2KEJQ.

5. Validation of Aggregate Consistency and Spatial Distribution

Validation involves three layers. (i) Aggregate count: the number of SSP per point must match the GHSL count; a difference greater than 0 indicates anomalies (e.g., non integer population values or casting issues). (ii) Attribute distribution: calculate the proportion of results for age, gender, education, and occupation, then compare these with frequency tables using MAE and RMSE (Equations 1 and 2), where y is the observed proportion and ŷ is SSP proportion in category i, and n is the number of categories. (iii) Compatibility of the spatial distribution of SSP with the built-up land positions according to ESRI Sentinel-2 Land Cover.

(1)

(2)

1. Determination of Scores per Attribute Class

The method for calculating scores per class follows the schema you established in the table. To ensure consistency and facilitate multicriteria aggregation, all scores are expressed within the range [0, 1]. If an initial score column falls outside [0, 1], linear normalization can be applied without altering the actual ordering: sᵢ(norm) = (sᵢ − min sᵢ) / (max sᵢ − min sᵢ). Table 2 illustrates this basic schema.

Regarding age groups, the young (0–10 years) and elderly (>60 years) are assigned the highest scores due to mobility limitations, dependency, and increased support needs during emergency responses. The middle working age group (31–40 years) receives a lower score because of its relatively high responsiveness. This mapping also assumes that higher levels of education are linked to greater risk literacy, better access to information, and improved ability to use services—thereby reducing vulnerability scores. Jobs with low income stability or high physical exposure are assigned higher scores (e.g., category 3) because of increased economic and operational risks during flood disruptions. From a gender perspective, women are considered more vulnerable than men.

2. Determination of Vulnerability Parameter Weights

The established core preferences are: Age (very important), Occupation (important), Education (somewhat important), Gender (less important) Table 3. These preferences are linked to the Saaty scale as anchor vectors s = [7, 5, 3, 1]. The paired comparison matrix A = [aᵢⱼ] is created as anchor ratios, where aᵢⱼ = sᵢ/sⱼ, ensuring the matrix's transitivity and mathematical consistency Table 4.

The AHP weights w are derived from the normalized principal eigenvector of matrix A Table 5. When constructing anchor ratios, the maximum eigenvalue becomes λmax = n; as a result, the consistency index and ratio are: CI = (λmax − n)/(n − 1) = 0 and CR = CI/RI = 0 (with RI for n = 4 being 0.90).

3. Determination of Vulnerability Index

Vulnerability (V) is the tendency of a population to experience adverse impacts when exposed to hazards, influenced by demographic characteristics and socio-economic conditions. In this context, V does not measure the likelihood of floods (hazard domain), but rather the internal capacity of individuals to anticipate, absorb, and recover. The scope of indicators is limited to attributes directly available from the SSP-namely, age, gender, education level, and employment status -where Si is the score for the i-th attribute category and Wi is the weight for the i-th attribute (Equation 3).

(3)

Risk analysis combines hazards (H) and vulnerabilities (V) to evaluate the potential impact of flooding on populations. We use a complementary approach: (1) Relative Risk Index: R = H × V, with H and V scaled from 0 to 1.

The SSP modeling result is shown in Figure 4. The resulting SSP includes 4,459,936 individual entities geolocated on a 100 m grid in the Upper Bengawan Solo watershed, derived from GHS-POP R2023. Each individual is assigned directional jitter within a radius of ≤100 m of the pixel center and is assigned attributes such as age, sex, education, and occupation, which are calibrated against official district or city level statistics via marginal alignment (IPF) until they are consistent.

The results of attribute validation are shown in Table 5 and Figure 5. Based on these results, the errors for each attribute are as follows: education (MAE 4.89%; RMSE 5.19%), occupation (MAE 6.52%; RMSE 6.93%), gender (MAE 0.0103%; RMSE 0.0103%), and age (MAE 0.0137%; RMSE 0.0163%). Notably, gender and age are nearly identical to the observed values (a difference of about 0.01), highlighting the success of mapping age bins and maintaining consistent category labels. This alignment is important because these attributes are often key factors in flood vulnerability indices (such as the proportion of children and the elderly as vulnerable groups); therefore, the stability of their values helps prevent errors from affecting the overall risk results.

Regarding education and occupation, the MAE of 4.9–6.5 reflects different dynamics: the SSP distribution closely matches the observed distribution (the pattern of comparisons between categories is preserved), but there is an offset in some high frequency categories. In practice, this often happens when category ontologies are not perfectly aligned, such as combining High School/Vocational School into a single source and distinguishing High School from Vocational School in another. In these cases, a few percentage points difference in one category is usually offset by a comparable difference in a nearby category.

The SSP validation results against GHS-POP for 30 sample pixels show a plausible pattern Figure 6. The relationship between the two is strong (r = 0.935), yet there is a systematic overcount with a bias of +10.4 people per pixel (MAE 12.93; RMSE 22.22), and a 95% agreement limit that remains wide (−28.74 to 49.54), indicating significant residual uncertainty in some pixels. Structurally, the regression slope is less than 1, and increasing error at higher density bin points to heteroskedasticity and the potential for spatial “bleeding” caused by jitter points extending beyond pixel limits and non mass preserving allocations per pixel. The SSP has successfully captured the density gradient but needs scale calibration, and jitter restrictions within pixel polygons and allocations hinder consistency at 100 m. With these adjustments, especially for aggregations ≥300 m, the accuracy is expected to improve, but expanding the sample will be necessary to stabilize the metrics and reduce the Level of Acceptance range.

Additionally, the distribution of SSP locations compared to ESRI Sentinel-2 Land Cover (10 m) shows high consistency: out of a total of 4,235,296 points, 91.95% (3,894,372) are on built-up land and 8.05% (340,924) on non-built-up areas. Treating this as an estimate of binomial proportions, its statistical uncertainty is minimal (95% confidence interval approximately 91.92–91.98%), indicating a strong “location accuracy” signal for the SSP. Substantively, the roughly 8% off-built-up figure remains reasonable given (i) thematic errors in land cover maps (such as commission/omission of built classes, especially in rural residential areas mixed with vegetation), (ii) the impact of the 10 m boundary—misregistration and mixed pixels at settlement edges, and (iii) jitter or bleeding of SSP points near class boundaries. These results strongly support the conclusion that the placement of SSP aligns with built-up spaces; the remaining 8% likely results from class boundary issues and thematic uncertainties, which could be minimized through small buffering, datum/registration alignment, or the use of building footprints in the test zone.

The demographic vulnerability analysis indicates a very “low” profile (Figure 7 and Figure 8). The average weighted population is 0.2160 (scale 0–1), with 99.63% of the population categorized as Very Low (44.43%) and Low (55.20%), while only 0.37% falls into the Moderate category. The AHP weights are consistent (CR=0), confirming that age is the primary factor (w=0.5333), followed by occupation (0.2667), education (0.1333), and gender (0.0667).

Consequently, pockets of vulnerability mainly develop in areas with higher proportions of children and the elderly, or in areas dominated by low adaptive capacity jobs (informal/daily workers) and low education levels, even though their overall contribution is smaller than their share of the population. Spatially, the impact is that vulnerability hotspots are limited but operationally significant: although they account for only 0.37% of the population, the intersection of moderate pockets with high hazard zones (H) could prioritize interventions such as increasing household preparedness, improving access to emergency services, and strengthening social networks. These results highlight that most regions possess reasonable demographic resilience reserves; however, risk management must continue to focus on micro locations with vulnerable age groups and job/education vulnerabilities to deploy mitigation resources effectively.

The distribution of the population's flood risk based on SSP (approximately 4.23 million entities) is heavily concentrated in the Low Risk category (3,391,760; 80.21%), with a substantial portion in the Moderate Risk category (830,469; 19.64%). Meanwhile, the tail distribution is nearly zero for the Very Low (5,026; 0.12%) and High (1,602; 0.04%) categories, indicating that overall exposure is primarily in the low risk class. However, about 832,071 (±19.68%) of the population falls into the moderate high risk category, which should be the focus of priority interventions, such as improvements to early warning systems, micro infrastructure protection, and spatial planning in hotspots (see Figure 9).

The Getis–Ord Gi analysis reveals that risk hotspots form significant clusters along river network corridors and in urban core areas Figure 10. The directional pattern along water pathways indicates the dominance of river flooding in floodplains, while the clustering in urban zones reflects pluvial flooding caused by high impervious surfaces and limited drainage capacity; both are worsened by high exposure (population/asset density) and local vulnerabilities (informal settlements, fragile buildings, service inequality). Therefore, the interaction of H×V explains the consistently positive z-score of Gi observed in riverbanks and urban pockets. Operationally, this involves addressing setbacks and restoring floodplain connectivity, improving drainage systems and green-blue infrastructure in cities, and prioritizing early warning systems and protecting critical facilities in hotspots. To ensure the accuracy of these conclusions, sensitivity testing should be performed on spatial neighbor weighting schemes (including river network based methods), along with corrections for multiple testing and cross validation using historical flood records.

The metrics findings show that the synthetic population distribution (SSP) has been accurately reconstructed regarding attribute spread and spatial distribution, although some residual errors remain in several categories. In high resolution population mapping, harmonizing category definitions and maintaining consistent classification schemes are essential steps to reduce artifact related residual bias (Rubinyi et al., 2022).

Substantively, age and gender attributes exhibit practical errors of approximately 0.01% (≈0.01%). This confirms the success of age bin mapping, the consistency of category labels, and the normalization of proportions. High accuracy in these two attributes is essential because the proportions of children and the elderly often influence response capacity and evacuation needs, while demographic structure also impacts risk projection by incorporating population dynamics and climate change (Shu et al., 2023). Therefore, the propagation of errors to the overall vulnerability index should be minimal. Conversely, in education and employment, a mean absolute error (MAE) of around 5-7 percentage points was observed. While this is generally considered a good fit for categorical shape, some dominant classes still show offset levels. From a policy standpoint, such deviations tend to affect thematic evaluations-such as post flood job training programs or school level disaster literacy initiatives-more than they impact spatial risk rankings between zones, provided these two attribute weights do not dominate the composite index.

From a utilization standpoint, current performance is sufficient for operational risk mapping at a resolution of ≥100 m (or village level). The research on composite based flood risk assessment models indicates that aggregate outcomes are fairly resilient to minor variations in less sensitive components, as long as the most sensitive factors-in this case, age and gender-are accurately measured. However, it is also important to note that fluctuations in population datasets alone have been shown to alter flood exposure estimates across different regions (Li et al., 2024)(Zhang et al., 2025)(Zhang et al., 2025). Consequently, documenting and monitoring residual patterns remain essential for quality control, especially in urban and peri urban areas, which often display greater social and economic diversity.

Regarding sources of bias, three main factors should be noted. First, misalignment between category lexicons (observed vs. SSP) could cause discrepancies of several percentage points in large categories. Second, imperfections in conditioning employment attributes on demographic predictors and regional contexts (e.g., age, education, and urban rural typologies) may lead to small but systematic distortions. Third, rounding in reference data can mathematically amplify differences in high frequency classes without changing the overall shape of the distribution. Two low cost steps are recommended that do not alter the SSP architecture: (i) post synthesis ranking or Iterative Proportional Fitting (IPF) in each territorial unit to "attach" margins to the targets; this method is common and effectively reduces the mean absolute error (MAE) of categories to <3 percentage points (Nejad et al., 2021), and (ii) one to one category ontology harmonization (consistently merging or crossing "similar" classes), as recommended in high resolution synthetic population mapping (Rubinyi et al., 2022). Both steps improve accuracy without sacrificing spatial consistency or computational efficiency.

As an additional measure, two more tests are recommended. First, a sensitivity analysis of the weight assigned to the vulnerability index: assessing the stability of spatial rankings across reasonable variations in educational and employment weights will confirm the robustness of the results. Second, constructing a residual map (observed – SSP, in %) for each main category can help identify spatial bias patterns—such as residual concentration in urban cores versus peri urban areas—that may influence the next ranking stage. This can also help distinguish between methodological bias and actual demographic heterogeneity (Zhang et al., 2025). In relation to medium to long term planning, combining demographic structure accuracy with these light balancing mechanisms is also important for a risk projection framework that considers future population dynamics and climate related hazard scenarios (Shu et al., 2023). Overall, these results support the conclusion that the SSP has adequately captured key demographic structures for area prioritization and operational level mitigation planning, with room for improvement through IPF/ranking and category harmonization. This approach aligns with industry best practices and offers a way to increase precision without adding computational complexity or model architecture overhead ((Nejad et al., 2021); (Rubinyi et al., 2022); (Shu et al., 2023); (Zhang et al., 2025); (Zhang et al., 2025)).