New Interpretable Framework for Clustering Spatial Hydrogeochemical Data and Assessing Groundwater Quality and Chemical Evolution Factors: A Topological Synthesis Approach

1. Introduction

The groundwater analyzed in this study is the result of a complex interaction of natural and anthropogenic factors and is widely used by the local population for domestic, agricultural, and industrial use. The study area is 150 km² and was declared an environmental emergency zone in August 2000 owing to increasing groundwater contamination and deteriorating health conditions in some communities. Therefore, our research prioritizes the examination of the factors controlling the chemical composition of groundwater in the Southern Bug and Sinyukha interfluves area (Ukraine). The factors governing the groundwater chemistry of natural water has been focus of many studies [1,2,3,4], with their number having increased considerably in last two–three decades [5,6,7,8,9,10,11,12,13,14]. Many of them involve graphical methods, such as the Piper diagram, Gibbs diagram, ion ratio analysis, and rock weathering plots, as well as statistical techniques, such as principal component analysis (PCA) and cluster analysis. When combined, these methods help evaluate the quality of groundwater and determine the geological and anthropogenic impacts on its chemical composition, including weathering and dissolution of minerals, evaporation, atmospheric precipitation, ion exchange and pollution [12]. In recent years, the utilization of machine learning and artificial intelligence methods for groundwater quality analysis has gained popularity [7,9,13,14,15,16,17,18,19]. Such methods may involve cluster analysis to identify representative clusters (also known as water types or hydrogeochemical facies) reflecting the processes generating natural variability in hydrogeochemical parameters. These representative clusters, which help identify key chemical trends, can provide insight into the heterogeneity and connectivity of aquifers, as well as the physical and chemical processes that control water chemistry. Specifically, Nguyen et al. [7] applied the self-organizing map (SOM) framework for the first time to decipher complex seasonal and spatial hydrogeochemical characteristics within confined aquifers. By coupling SOM clustering with classical Gibbs diagrams, they successfully identified how recharging factors shift water types from dry to rainy seasons, and distinguished between rock weathering, evaporation, and marine intrusion processes. Using the SOM framework, Nakagawa et al. [9] classified groundwater hydrogeochemical data into eight clusters, identifying three main water types based on salinity: high salinity, low salinity, and freshwater. Liu et al. [16] demonstrated the effectiveness of integrating SOMs and multivariate statistics for elucidating the complex interaction between natural hydrogeological processes and human activities affecting groundwater geochemistry in an urban environment. However, most studies utilize raw ion concentrations, pH, and total dissolved solid (TDS) as isolated parameters within the SOM framework, which does not allow for the explicit identification of dominant geochemical processes. To address this limitation, the input parameters in our study are ion combinations rather than individual concentrations, allowing the neural network to operate with process indicators rather than simple concentrations. Moreover, the proposed method, based on the projection of stoichiometric equality lines and the principal component vectors across SOM component planes, allows for the identification of dominant (competing) processes in a single format, while simultaneously revealing the underlying interrelationships. Furthermore, unlike other methods, our method, which combines log-ratio analysis, PCA, SOM, and graphical cross-validation (GCV), focuses on identifying physical stoichiometric equality isolines as well as main component isolines, converting the algorithmic “black box” into a transparent model for expert analysis. For example, the application of Support Vector Machines (SVM; [20]) via kernel functions (the “Kernel Trick”), according to [21], maps data into infinite-dimensional feature spaces to identify an optimal separating hyperplane. However, as noted in works on Interpretable Machine Learning (XAI) (e.g., [22]), such transformations render the visual tracking of physical thresholds impossible, effectively turning the process into a “black box”. The primary advantage lies in the ability to delineate cluster boundaries—uniquely facilitated by the intersection points of stoichiometric equilibrium lines and principal component vectors projected directly onto the component planes—while simultaneously distinguishing between anthropogenic and natural factors driving groundwater chemistry. Furthermore, for an independent verification of the topological conclusions and a robust assessment of hydrogeochemical process intensities, we incorporated GCV plots into our model and constructed distribution profiles of normalized balances (f_i). Available hydrogeochemical research on the groundwater in the studied Pleistocene aquifer has primarily examined the chemistry of major ions, the steady-state flow, and the ecological–hydrogeological state of the groundwater [23,24]. In our previous study [25,26] we revealed the processes and factors controlling the chemical evolution of groundwater in modern and Pleistocene deposits in the Southern Bug Sinyukha interfluve area. This study investigates the role, dynamics, and interaction of these factors as well as their impact on TDS and pH. The studied unconfined aquifer, owing to its proximity to the surface and location in low-lying, topographically depressed areas, is one of the causes of groundwater pollution. Furthermore, groundwater functions as a protective barrier and source of pollution for the fissure waters of the underlying Precambrian crystalline aquifer, which is used for supplying drinking water. In this context, comprehensively examining the chemical composition of groundwater and the factors controlling its formation is a priority for the study area and wider region. The objective of this study was to cluster groundwater chemistry controlling factors and investigate their impact on groundwater quality in the unconfined aquifer in the Southern Bug and Sinyukha interfluve area using log-ratio analysis, PCA, SOM, and GCV. This study's findings will yield significant insights into the spatial hydrogeochemical properties of groundwater within the Pleistocene unconfined aquifer of the study area. Furthermore, the proposed method outperforms standard “black box” solutions in interpretability and provides a robust framework for verifying the boundaries of hydrogeochemical data clusters, which will be essential for accurately evaluating the quality and spatiotemporal evolution of groundwater chemical composition across the region.

2. Materials and Methods

2.1. Study Area

The study area is in the southern part of the Ukrainian reservoir of strata-fissure waters on the territory of the Southern Bug and Sinyukha interfluves (Figure 1). This area is a typical agricultural and industrial area, located in the steppe zone of southern Ukraine, and belongs to the arid zone of insufficient moisture. Here, the dry season sets in during the second half of June and the first half of November, with precipitation decreasing, and evaporation reaching peak values.

The study area is characterized by gently rolling plains dissected by river valleys and gullies, with the slopes being eroded by gullies and ravines. Absolute relief marks range from 130 to 247 m. The lithologic–stratigraphic complex of the Pleistocene and Holocene includes the following [24]: an aquifer in the alluvial, proluvial–deluvial sediments of the river floodplains and the gully bottoms; and a locally weak aquifer in the aeolian–deluvial deposits of the watershed plateau (Figure 1). The lithology of the Pleistocene aquifer comprises mainly sand, sandy loam, loam sometimes with carbonate, and clay with gypsum crystals and silt sometimes with crystalline rocks. These sediments comprise mainly quartz, calcite, dolomite, gypsum, and feldspar. Weathering and dissolution of these minerals is a possible source of groundwater enrichment with different ions. The depth of the groundwater level varies from 0.0–5.0 to 10.0–15.0 m. The thickness of the water-filled layer is 0.5–3.0 m—sometimes up to 5.0 m. Groundwater recharge is facilitated by atmospheric precipitation, surface water, and inflow of pressure water from the crystalline base through decompaction zones. During low-water periods, discharge occurs into rivers and underlying aquifers—provided there is no low-permeability rock present. The groundwater is widely used by the local population in numerous boreholes for domestic and agricultural purposes. In this area there are about 35 anthropogenic sources of groundwater pollution such as industrial and domestic landfills, pesticide and fertilizer storage, various livestock farm complexes, settling tanks, wastewater and solid waste discharges, and so on. These anthropogenic items have an impact on the geological environment—especially concerning its hydrosphere. Incisions, such as river valleys and ravines, which host populated regions and their agricultural complexes, are particularly prone to technogenic influences. In addition, groundwater serves as a protective screen and source of contamination for the underlying Precambrian Crystalline Aquifer, which is used for supplying drinking water. These aquifers are hydraulically connected mainly in gullies or river valleys.

2.2. Sample Collection and Analysis

Groundwater samples were processed for a short chemical analysis, estimating the content of nitrogen group compounds and spectral analysis of the dry residue [24]. Sampling was performed at eight monitoring wells and 25 dug wells, which are located in the territory of Boleslavchik Village. Monitoring wells were drilled along the slope of the central part of the village (north–south profile; profile line P-1–P-6 in Figure 1) and along the bottom of the gully (east–west profile; profile line P-7–P-8 in Figure 1). Sampling was conducted during the late-summer low-flow stage (specifically from August 30 to September 21, 2000) in accordance with the guidance for sampling, preservation and processing of groundwater samples recommended by the State Agency for Water Resources of the Ministry of Environmental Protection and Natural Resources of Ukraine. Chemical analyses were performed by the accredited Research Central Laboratory of the Black Sea Coastal Regional Geological Enterprise “Prichomomor DRGP” (Odessa, Ukraine), following the standard examination methods for drinking water quality assessment (aligned with DSanPiN 2.2.4-171, [27]). Temperature and pH were measured at the sampling site to account for possible physicochemical changes. All samples were filtered through a 0.45-μm membrane filter and collected in clean, dry plastic bottles made of polyethylene or high-density polytetrafluoroethylene.

2.3. Data Used and Pre-Processing

Primary data used in our study comprised 91 groundwater samples with 10 parameters including pH, TDS and major ions ($C{a}^{2+}$, $M{g}^{2+}$, $N{a}^{+}$, $K^{+}$, $HC{O}_3^{-}$, $C{l}^{-}$, $S{O}_4^{2-}$, and $N{O}_3^{-}$). In terms of the charge balance error (CBE), groundwater quality was satisfactory (CBE < ±2.5%). The average chemical composition of groundwater is listed in Table 1. According to Piper diagrams in our previous study, the dominant hydrochemical types of the studied shallow groundwater are sodium bicarbonate, sodium sulfate and sulfate–magnesium–sodium, sulfate–calcium–sodium types, and sometimes nitrate–hydrocarbonate and nitrate–chloride–sodium types [26,27].

Building upon our previous work defining hydrogeochemical facies [26,27], this study employs a modified approach combining log-ratio analysis, PCA, SOM, and GCV to analyze the evolution of the dynamic system. For the application of SOM, we adopted six dimensionless parameters—i.e., f₁ – f₅ (Table 2, column 2) and ln(pH). Parameters f₁ – f₅ can be interpreted as factor or criteria of hydrogeochemical dynamic similarity. Previously, we identified these factors using graphical methods [1,4,12,28], including Gibbs diagrams, as well as binary plots of the major ions [25]. In this study, these factors were used as input parameters for the neural network. The concentrations of all ions for the factors are expressed in milliequivalents per litter (meq/L, Table 2). To satisfy multivariate normality and reduce outlier effects [29,30,31], a log-transformation was applied to both the parameters f₁ – f₅ and pH. Column 5 in Table 2 is the result of a generalization of numerous scientific results [1,4,12,26,28].

Importantly, to ensure conceptual clarity, the specific values on the SOM feature maps (e.g., 0.31 for f₁ and 0.14 for f₃ in Table 2) are defined as normalized stoichiometric 1:1 thresholds (Th_norm). These values do not represent classic hydrogeochemical characteristics for mineral equilibrium derived from laboratory experiments. Instead, each Th_norm represents dataset-specific topological projection of the physical 1:1 ionic equivalence line (ln(Y/X) = ln(1) = 0)). The spatial representation of this 1:1 boundary on the SOM is defined as the stoichiometric equality line. Since each log-ratio parameter has its own unique empirical minimum and maximum range, the physical zero maps onto a different normalized value for each individual feature. Following log-ratio pre-processing, each f_i was treated by PCA and SOM as an independent feature. Min–Max scaling was used to preserve the original data distribution shape and correlation structures within the feature space F. Many researchers apply isometric (ILR) or centered (CLR) log-ratio transformations to classify water types based on percentage compositions [32,33]; however, in this study, we did not adopt these approaches. Since our data represent arbitrary geochemical intensity indicators rather than parts-of-a-whole (proportions), CLR or ILR would be counterproductive, confounding specific thresholds (e.g., f_i =0 for stoichiometric equality) with the variance of other features. In our approach, the logarithmic space ensures a unique global equilibrium at f_k = ln(1) = 0, preserving the physical interpretability of SOM isolines to distinguish competing processes. Furthermore, the model is robust against translation or scaling, while shifting the normalized stoichiometric 1:1 threshold does not distort the data. Consequently, multiplicative relationships can be transformed into additive ones.

2.4. Efficiency of SOMs for Hydrogeochemical Data Analysis: Configuration and Training

The SOM framework, developed by Kohonen [34,35], is a type of artificial neural network trained through unsupervised learning. The SOM algorithm projects high-dimensional data into a low-dimensional latent space, creating a regularized map based on data similarity and competitive neural principles. Unlike traditional deterministic models, the SOM framework preserves the inherent universality and nonlinear complexity of relationships between hydrogeochemical characteristics. Geological similarity is predominantly approximate, as most geological processes are multifactorial and operate across vast spatiotemporal scales. The dual nature of geological formations—i.e., “locally heterogeneous, yet globally homogeneous”—requires a tool capable of capturing generalized features (the “homogeneous in total”), while accommodating local variances. In this study, the SOM provides a synergistic effect, integrating input parameters into a unified topological structure that mirrors the evolutionary complexity of groundwater composition. Furthermore, the topological consistency of the SOM allows for the overlay of theoretical stoichiometric equality lines and PC projections, transforming the visual map into a robust tool for identifying phase transitions, hydrogeochemical bifurcation points, and so on. The SOM is a robust tool for clustering and data mining facilitating the extraction of essential information in hydrogeochemistry [7,9,10,13,14,15,16,18,30,36,37,38,39,40]; it maps high-dimensional input space onto reference (weight) vectors through three stages: node competition, winning node selection (best-matching unit), and weight updating via Euclidean distance [41]. Designing an appropriate SOM structure—specifically the number of nodes (m) and the grid’s aspect ratio—is crucial [42]. In this study, the total number of nodes was determined using the heuristic rule—i.e., $m=\sqrt{n},$where n is the number of samples. The grid’s side lengths were set according to the square root of the ratio between the two largest eigenvalues ($\sqrt{{\lambda }_1/{\lambda }_2}$) derived from PCA [36,43]. All log-ratio data were normalized using Min–Max scaling to a range of [0, 1]. Given the limited dataset (n=91), linear initialization was preferred [30,41]. Training was performed in batch mode for 6000 iterations using a Gaussian neighborhood function. The initial learning rate and neighborhood radius were 0.37 and 2.40, respectively. Map accuracy was evaluated through the quantization error (QE) and topographic error (TE) [35]. For 91 groundwater samples, the optimal map size was calculated as 8×6 (48 nodes). The resulting QE and TE values were 0.13 and 0.00, respectively. Adjacent node distances were visualized using a U-matrix.

2.5. Multivariate Statistical and Graphical Methods. Topological Invariance Analysis

Hierarchical Cluster Analysis (HCA), using Euclidean distance and Ward’s linkage, grouped six hydrogeochemical parameters (Table 2). The optimal cluster count, determined by a dendrogram phenon line, was validated using the Silhouette coefficient [44,45]. PCA revealed the dataset’s latent structure via orthogonal variables [43]. In this study, a key methodological innovation is the integration of the Th_norm line and PC1 and PC2 projections into the SOM structure. Th_norm acts as a topological boundary distinguishing process dominance (f_i > Th_norm or f_i < Th_norm). Owing to log-ratio transformation and Min–Max scaling, the SOM topology remains invariant to scalar multipliers of initial concentrations. This allows dynamic rethinking of process dominance zones—by shifting the stoichiometric threshold Th_norm—without recalculating the model. To verify these conclusions, GCV analysis and normalized Δf_i balance profiles were constructed against TDS and pH, with the line ln(Y/X) = 0 becoming the universal axis of stoichiometric equality (Th_norm), dividing the zones of unconditional dominance (> 0 for Y_i; <0 for X_i), and enabling the evaluation of parameter dynamics within SOM planes. HCA and PCA were performed using a modified SciPy package (https://scipy.org/, accessed on September 2024) in Python environment.

2.6. Topological Synthesis of PCA Projections and Hydrogeochemical Stoichiometric Equality Lines on the SOM Surface. Nonlinear Projection of Balance Relations

We define the feature space$F\subseteq R^N$ as a set of N logarithmic parameters (in this study, N = 6): $F={\left\{f_k\right\}}_{k=1}^N$, where each feature f_k is calculated as the logarithm of the ratio of a pair of initial concentrations of hydrogeochemical elements Y_i, X_j (measured in meq/L) or as the logarithm of the hydrogen exponent (pH):

f_k = ln(Y_i /X_j)  or  f_k = ln(pH),

(1)

where i, j ∈ {1,…, n}, and n is the number of initial components in the system. Despite f_k being derivative functions, they are treated by the mathematical models (i.e., SOM and PCA) as statistically independent coordinates of the input vector. In additive balance space, the value f_k = 0 strictly corresponds to the point of physical equilibrium (Y_i = X_j), justifying the application of the Euclidean metric in the SOM and PCA algorithms for the objective assessment of process dominance. Each feature f_k was transformed into a dimensionless value ${\tilde{f}}_k$ within the range of [0, 1] according to the following expression:

$${\tilde{f}}_k={\displaystyle \frac{f_k-f_{k,\min }}{f_{k,\max }-f_{k,\min }}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\tilde{f}}_k\in \left[0,1\right].$$

(2)

where f_k,min and f_k,max are the minimum and maximum values of the feature f_k across the entire dataset, respectively. The SOM is defined as a nonlinear mapping Φ from the high-dimensional space of normalized features $\tilde{F}$ onto a discrete two-dimensional lattice of neurons M (latent space) [35,46]:

$$\Phi :\tilde{F}\to M,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}M=\{{\mathbf{w}}_1,{\mathbf{w}}_2,\dots ,{\mathbf{w}}_m\},$$

(3)

where w_p is the weight vector of neuron p (p = 1,…, m), and m is the total number of neurons. The Th_norm on the component plane for feature f_k is defined as the set of neurons where the weights satisfy the condition of physical equality between components (Y_i =X_j):

$${\mathit{Th}}_{\mathit{norm}}=\left\{p\in M\left|\right.w_{p,k}=s_k\right\}$$

(4)

The threshold value s_k in the normalized space is calculated as follows:

$$s_k={\displaystyle \frac{0-f_{k,\min }}{f_{k,\max }-f_{k,\min }}}.$$

(5)

where 0 corresponds to the value of ln(Y_i /X_j) at equilibrium, and f_k,min and f_k,max are the extrema of the initial log-ratio feature f_k. The principal component projections onto the SOM surface are calculated as the dot product of the neuron weight vectors and the PCA loading vectors:

$${\mathit{V}}_{\mathit{P}\mathit{C}\mathbf{1}}\mathbf{(}\mathit{p}\mathbf{)}=〈{\mathbf{w}}_{\mathit{p}},{\xi }_{\mathbf{1}}〉={\displaystyle \sum _{i=1}^N{\mathrm{w}}_{\mathit{p},\mathit{i}}}\cdot {\xi }_{1,\mathit{i}},$$

(6)

$${\mathit{V}}_{\mathit{P}\mathit{C}\mathbf{2}}\mathbf{(}\mathit{p}\mathbf{)}=〈{\mathbf{w}}_{\mathit{p}},{\xi }_{\mathbf{2}}〉={\displaystyle \sum _{i=1}^N{\mathrm{w}}_{\mathit{p},\mathit{i}}}\cdot {\xi }_{2,\mathit{i}},$$

(7)

where ξ₁ and ξ₂ is the first and second eigenvectors (loadings) of the covariance matrix of the normalized data$\tilde{\mathbf{Z}}$, and N is the number of features. For GCV analysis, the transition to logarithmic relations in normalized space is conducted as follows:

$$\Delta f_{ik}={\tilde{f}}_{ik}-s_k.$$

(8)

Eventually, the SOM map is segmented into physically distinct operational regimes: domains Ω₁$\{{\tilde{f}}_k>s_k\}$ and Ω₂$\{{\tilde{f}}_k<s_k\}$, representing the dominance of competing hydrogeochemical processes. The modified software MiniSom (https://pypi.org/project/MiniSom/, accessed on September 2024) was utilized for constructing the SOM model.

3. Results and Discussion

3.1. SOM and Clustering Results: SOM-Based Visualization of Multivariate Data

The component planes for the 48 reference vectors (nodes) of the six log-ratio parameters (f₁, f₂, f₃, f₄, f₅, and pH), normalized to a range of [0, 1], are presented in Figure 2. HCA was applied to the SOM prototypes using Ward’s linkage algorithm and the Euclidean distance measure. The first principal component accounts for 42.6% of the total variance, largely driven by feature f₃ (loading value of 0.82). The second principal component accounts for 34.4% of the total variance, with f₅ being the parameter most closely associated (loading value of 0.80). The optimal number of clusters (three) was determined by maximizing the silhouette coefficient (Figure 3a), which reached its peak score of 0.40 with three clusters. This optimal partitioning is visually confirmed by the phenon line in the dendrogram (Figure 3c), drawn at a linkage distance of 17. The resulting low QE of 0.14 and perfect TE of 0.00 further indicate optimal model training and high mapping quality. The orange isoline on the component planes indicates the projection of the hydrogeochemical stoichiometric equality line (Th_norm). The threshold values (s_k) for Th_norm for each SOM plane were calculated according to Equation (5). The green solid isolines represent multiples of the s_k threshold. The red and blue dotted lines represent the projections of PC1 and PC2 onto the SOM plane, calculated using Equations (6) and (7), respectively.

The U-matrix (Figure 2g) visually confirms high distances between the three identified clusters, which delineate physically distinct hydrogeochemical facies. The coincidence of Th_norm isolines with zones of high inter-neuron distance on the U-Matrix proves that the equilibrium line acts as a natural topological barrier. A preliminary visual qualitative analysis of the component SOM maps (Figure 2) and clustering results reveals the specific characteristics of each cluster. For example, Cluster 1 (left of Figure 2h) is associated with high values of f₁, f₂, f₃ and f₅, primarily observed in the upper-left regions of the respective component planes (Figure 2a–c, 2e). Cluster 2 is characterized by high values of f₁, f₃, and f₄ (the zones below Th_norm), relatively high values of f₅ (lower part of Cluster 2), and high values of pH (hexagon almost in the center of the pH SOM), with low values being concentrated in the right corner. Cluster 3 is characterized by high values of f₂, f₄ (the zones below Th_norm), and f₅ and relatively high values of f₁ (to the right of the Th_norm). Also highlighted here are the low pH values for a small number of samples in the lower right corner of the SOM map (Figure 2f). Since SOM clusters data based on the similarity measure of multidimensional feature vectors, the clustering of factors (features f₁– f₅) in this study can also be interpreted as a clustering of hydrogeochemical dynamic similarity criteria.

3.2. SOM and Clustering Results: Topological Discontinuity and Singularity Analysis

This section synthesizes system-scale analysis (PC1/PC2 projections) and local hydrogeochemical constraints (Th_norm) to implement phase boundary mapping. The Th_norm isolines on the SOM plane visualize objective hydrogeochemical boundaries dictated by shallow groundwater dynamics, where comparing these lines with PC vectors differentiates between equilibrium processes and external or anthropogenic drivers. Since SOM weights are normalized within the range of [0, 1], Th_norm acts as a fundamental topological boundary dividing the map into two phase zones (Table 2; Figure 2). The dominance of the process is characterized by the zone occupying a larger map area where weight vectors w_p,k exceed the specific threshold s_k of the specific feature. Equal zone areas signify process equivalence, and the shift of the Th_norm line on the SOM plane is mathematically equivalent to rotating the equilibrium ray in the original Cartesian coordinate system. The topological discontinuity in the PC1 projection represents a first-order mathematical singularity, corresponding to a sharp or inverse gradient in weight vectors where smooth evolutionary trajectories are interrupted by a hydrogeochemical barrier. Its physical nature as a genuine hydrogeochemical barrier—rather than a multidimensional projection artifact—is confirmed by contour crowding (Figure 2a and Figure 2c) and topological loops or intersections (Figure 2b, Figure 2d, Figure 2e). Contour crowding before the discontinuity and abrupt value changes indicate a sharp barrier, typically representing water mixing, phase transitions (e.g., rapid pH shifts), or intense anthropogenic pressure (point-source pollution or localized over-pumping). This is confirmed by Figure 2e and Figure 2f, where neurons containing samples 57, 72 and 75 show a sharp jump in pH, making linear interpolation of isolines impossible. The area between PC1 discontinuities represents a maximum entropy zone where dissolution forces and anthropogenic pressures counteract, creating inverted U-shape patterns (Figure 2a). According to [47,48], this inverse U-shaped pattern of PC1 and isolines that are proportional to Th_norm indicates a topographic break in a nonlinear data manifold, which is not a dimensionality artifact [49]. In a hydrogeochemical context, it confirms a sharp phase transition where the primary vector (PC1) fails to continuously approximate the process, visualized as a first-order discontinuity. For instance, limit analysis (one-sided limits) in Figure 2a reveals a functional jump typical of physical–chemical barriers: Values increase on the left side of the break (0.50, 0.60, 0.70, 0.75, and 0.80), whereas the right side begins at a different order (0.70, 0.75, and 0.80). Furthermore, the PC1 isoline runs parallel to Th_norm, defining the first cluster boundary and proving it is a determinant for PC1. Thus, PC1 isoline segmentation demonstrates system nonlinearity and stepwise hydrogeochemical evolution through critical barriers. Plotting the blue dashed PC2 isoline across all planes is essential to identify secondary processes masked by the dominant PC1 trend. The intersection of systemic axes (PC1, PC2) with Th_norm isolines enables a comprehensive topological model of phase transitions. In PC1 discontinuity zones, the variance of secondary factors (PC2) increases and overcomes the global trend resistance threshold, suggesting that the local physical–chemical potential of a specific feature (e.g., f₅) becomes dominant at critical points. In the middle upper layer of SOM maps (Figure 2), where the left branch of PC1 and the PC2 line approach orthogonally to the boundary separating Clusters 1 and 3, PC2 intensity resists the global PC1 trend and reaches a critical level. At this junction, the global trend weakens, causing an isoline break where the secondary process (feature f₅) becomes sufficient to pierce the cluster boundary, transitioning the system into a new state. The alkaline core region, where PC1 and PC2 approach the boundary of Clusters 2 and 3 (Figure 2f), locally represents a bifurcation point where the system determines its evolutionary direction. Inter-cluster boundaries (Figure 2h) are determined by the approximation of principal components and Th_norm isolines. For instance, the boundary between Clusters 1 and 2 coincides with the PC1 trajectory and the multiple stoichiometric equality isoline Th_norm = 0.5. The boundary between Clusters 2 and 3 is defined by the PC1 vector up to the third row of the neural grid (prior to the PC2 and f₅ stoichiometric isoline intersection); after the intersection, it is dictated by the PC2 trajectory. Topologically, boundaries are fixed where PC1 and PC2 become collinear with stoichiometric equality isolines. Furthermore, the intersection points of isolines s_k, Th_norm, and PC2 with the PC1 global trend consistently pass along the cluster boundaries, thereby playing a constructive role in the formation of cluster boundaries. Such critical threshold points are clearly visible in Figure 2b, Figure 2d and Figure 2e and on the distance matrix (U-matrix) in Figure 2g. The system includes two functional relations based on PC1 loadings. The first group includes features with significant positive loadings (e.g., f₁, f₃, and the more complex f₅), where PC1 isolines do not intersect Th_norm as these parameters themselves dictate the PC1 trajectory. The second group consists of features with negative loadings (f₂ and, f₄) that oppose the global trend; their Th_norm isolines intersect the PC1 vector, creating bifurcation points and complex cluster boundaries. This interaction assesses the resistance of local geochemical balances to systemic forcing, marking structural transformation zones. Thus, this topological synthesis transforms the passive SOM visualization into a phase diagram.

3.3. Phase Boundary Mapping via ln(Y_i/X_i) Topological Synthesis of Feature Balances (f₁–f₅)

On the f₁ component plane (Figure 2a), most samples lie above the Th_norm line (s_k = 0.31; Table 2), indicating $(S{O}_4^{2-}+HC{O}_3^{-})$ dominance over $(C{a}^{2+}+M{g}^{2+})$ driven by sulfate dissolution or $N{a}^{+}/(C{a}^{2+}/M{g}^{2+})$ cation exchange. The area below Th_norm reflects carbonate dissolution, while proximity to the line indicates mixed mineral weathering. The convex s_k = 0.5 isoline (inverted U-shape) between discontinuous PC1 segments marks a phase transition from natural weathering to an anthropogenic regime. This shift features a superposition effect where sulfate dissolution increases ionic strength, triggering ion exchange, while localized pollution shifts the system from its primary PC1 trajectory, causing sharp pH fluctuations (Figure 2f). Alignment between the left PC1 branch and the s_k = 0.5 isoline shows that the latter determines PC1, confirming that PC1 segmentation reflects nonlinearity and stepwise evolution across discontinuities. This topological gap represents a hydrogeochemical barrier created by mineral dissolution, ion exchange, and wastewater input, which the system overcomes in a stepwise, nonlinear manner. The biplot (Figure 3b) confirms f₁ as a primary driver with the second-highest positive PC1 loading. The PC1 discontinuity signals critical instability in the $(S{O}_4^{2-}+HC{O}_3^{-})/(C{a}^{2+}+M{g}^{2+})$ growth process, where the system deviates from its trajectory, crosses an equilibrium zone, and restabilizes in Cluster 3. At the C1/C3 boundary (Figure 2f), converging PC1 and PC2 trajectories suggest a regime transition where dominance shifts to feature f₅ (maximum PC2 loading). The f₂ plane indicates dolomite dissolution and cation exchange beyond Th_norm lines, while the segment between s_k=0.59 and s_k =0.72 reflects calcite dissolution. Decreasing f₂($C{a}^{2+}/M{g}^{2+}$) alongside increasing f₃$(S{O}_4^{2-}/C{l}^{-})$ signals external sulfate inputs driving magnesium dominance via cation exchange. The negative PC1 loading of f₂, intersecting multiple stoichiometric isolines, confirms a shift toward a magnesium-dominated regime via dedolomitization and evaporative concentration. The f₃ map (Figure 2c) shows sulfate dominating over calcium, linking ion sources to sulfate mineral dissolution. High topological similarity of f₁, f₂, and f₃ in upper-left quadrants confirms synergistic processes in Cluster 1. This phase exhibits $S{O}_4^{2-}$ and $M{g}^{2+}$ enrichment with calcium suppression, linked to Na⁺/K⁺ ion exchange (Figure 2d). Absent halite and marine impact involving elevated sulfate with trace elements (Cr, Li, Ni, Ti, and Zr) suggests sulfide oxidation and anthropogenic influences. Feature f₃ has the maximum positive loading on PC1, driving the system’s main evolutionary vector. The alignment of PC1 parallel to isolines confirms that the system variability depends on the f₃ distribution. This isoline fusion with the PC1 projection marks the cluster boundary as a phase threshold defining sulfate regime stability. Conversely, PC2 projection intersecting isolines perpendicularly in Cluster 2 indicates that the secondary factor acts as a disruptor to the equilibrium. The f₄ component plane (Figure 2d) reflects sulfate and chloride dominance (likely from anthropogenic runoff and fertilizers) over silicate weathering, with $(S{O}_4^{2-}+C{l}^{-})>(N{a}^{+}+K^{+})$. A 180° inversion between f₄ and f₂, relative to f₁, indicates a facies shift from bicarbonate to sulfate–chloride–sodium types. Holding the second-highest negative PC2 loading, f₄ regulates local salinization independently of the primary PC1 trend. The Th_norm and PC1 projection intersection in the bottom-right quadrant shows this feature’s physical threshold resists the main evolutionary vector. Isoline refraction and intersections confirm high nonlinearity within Cluster 2, matching extreme weight gradients on the U-matrix (e.g., neurons with samples 39, 73, and 89; Figs. 2g, 2h). Intersections of the PC1 and PC2 projections with the Th_norm lines on f₄ and f₂ maps pinpoint zones where the system balances between $(N{a}^{+}+K^{+})$ and $(C{a}^{2+}+M{g}^{2+})$ and define the cluster boundaries. Finally, the Th_norm intersection with the right branch of the PC1 projection (Figure 2d), correlated with pH fluctuations (Figure 2f), marks the termination of the stable regime and a transition to a localized, different geochemical phase. The f₅ component plane (Figure 2e) illustrates the synergy between anthropogenic impact and decarbonization, holding the maximum positive PC2 loading. Technogenic anions $(N{O}_3^{-}+C{l}^{-})$ dominate over bicarbonates in 56% of samples, whereas the natural bicarbonate component prevails in the remaining 44%. An orthogonal intersection of the PC1 trajectory with Th_norm line at the boundary between Clusters 1 and 3 marks a pronounced geochemical barrier. The topological discontinuity of this line indicates a zone of intensive mixing and metamorphism where the system cannot evolve smoothly without surpassing a critical nitrate–chloride saturation threshold. Comparative analysis of the f₃ and f₅ planes reveals that within the phase transition zone, control is intercepted by the secondary process (PC2), making the nitrate–chloride complex the dominant factor in local variability. Consequently, feature f₅ serves as a primary indicator of technogenic metamorphism, driving localized shifts in groundwater composition. The f₆ (pH) component plane (Figure 2f) reveals an anomalous alkaline zone within Cluster 2, situated between discontinuous PC1 isolines. The coincidence of this node with the pH extremum identifies a critical point where anthropogenic factors intercept natural evolution. The convergence of the PC1 and PC2 trajectories above this alkaline core forms a bifurcation point, marking maximum instability where the alkaline trigger and the nitrate–chloride complex (f₅) drive the topological discontinuity of the primary vector. Conversely, the lower right corner shows a localized pH decrease, signaling an acid-oxidative impulse. This acidity promotes intensive ion exchange and calcium leaching, leading to magnesium and $N{a}^{+}+K^{+}$ dominance, reflected in high weights of alkali metals, magnesium, sulfate, and bicarbonate components. The localized anomaly and limited sample size in Cluster 3 confirm the technogenic origin, interpreted as point-source pollution.

3.4. Cross-Validation and Geochemical Characteristics of the Clusters

Figure 4a and Figure 4b illustrate the spatial distribution plots of hydrogeochemical features (f_i), TDS, and pH relative to a unified equivalence line. To enhance the visualization of the dynamics and variability of these features, the data were categorized into two groups: f₁, f₂, f₄ (Figure 4a) and f₃, f₅ (Figure 4b). Above the equivalence line, TDS dynamics are determined by technogenic factors, whereas below it, they are driven by natural processes. Features f₁, f₂, f₃ and f₅ show the highest fluctuation amplitudes, whereas f₅, f₂, TDS, and pH show a stable trend within Сluster 3. Positively correlated features with TDS show collinearity of isolines with respect to the PC1 vector (e.g., f₁ and f₃ in the border zone of Cluster 1).

Conversely, negatively correlated features show orthogonal or intersecting isolines relative to the PC1 vector, as seen for f₂, f₄, and f₅ in Cluster 3 (neurons with samples 22 and 26), where PC1 intersects Th_norm isolines corresponding to points 23 and 27 on the GCV (Figure 4a and Figure 4b). High isoline density marks zones of intense feature growth. Parallel isolines to the PC1 vector indicate a contribution to the upward trend of TDS and pH, corresponding to progressive mineralization reaching extreme values in the spatial graphs (Figure 4a and Figure 4b, points 61, 64, 65). Conversely, orthogonal intersections of dense isolines by the PC1 vector signal a correlation inversion, exemplified by the f₅ topology (Figure 2e), where the PC1 trajectory intersects Th_norm almost orthogonally, corresponding to an anticorrelation at point 51 (Figure 4b).

According to the monitoring data [23], in the study area surface water and deep groundwater exhibit only low or localized (point-source) contamination, in contrast to the shallow groundwater, which is heavily polluted. Since the study area is located within the southern part of the Ukrainian basin of stratum-fracture waters and represents a zone of intense water exchange, the pollution predominantly affected the areas with poor water exchange. This is evidenced by the elevated values of TDS and features f₁, f₃, and f₅ observed from the north to the east and southeast in some dug wells and monitoring wells located in low water exchange conditions. The monitoring results of the study area [23] show that the soil-grounds within the vadose zone remain largely unpolluted by nitrates and are non-saline. However, the areal soil pollution by heavy metals (Ti, Cr, and Li) and point-source pollution by (Ni and, Zr) are clearly traced. This represents fingerprint of technogenic pollution, where the active phase has migrated from the ground surface into the unconfined groundwater aquifer and has been partially immobilized within the zones of poor water exchange. A joint analysis of SOM planes and GCV, while considering also Figure 1, allowed the identification of the patterns of the spatial distribution of the three identified clusters. On the GCVs, samples numbered 1–60 correspond to groundwater from dug wells, whereas samples 61–91 correspond to monitoring wells located under conditions of delayed water exchange in the Pleistocene aquifer of the Southern Bug and Sinyukha interfluve area within the territory of Boleslavchik Village. Cluster 1 includes groundwater samples from both dug and monitoring wells. The analysis of SOM planes for features f₁–f₃ and f₅–f₆ (Figure 2 and Figure 5, respectively) shows that groundwater in the northern part of the study area (Figure 1) is contaminated in both dug and monitoring wells.

The highest contamination levels occur in monitoring wells P-1–P-3, indicating a primary contamination zone. This is reflected in neurons containing samples numbered 60–70 on SOM planes 2a, 2b, 2c, and 2e (and 61–71 on the GCV plot, Figure 5), which are characterized by high TDS value. Since f₃ and f₁ have high loadings on PC1, the PC1 component within Cluster 1 reflects the vector of intense sulfate influx, and Cluster 1 is identified as the formation zone of the sulfate–chloride–magnesium water type (or salinity water type). Cluster 2 characterizes groundwater contamination in the eastern and south-eastern parts of the study area for both well samples (Figure 1) and dug wells. High pH and TDS values are observed in dug wells 24 and, 31, located in the southern part of the study area (neurons with samples 51–59 in Figure 2h; samples 52–60 in Figure 6), and extremely high values are observed in monitoring wells P-4 – P-6 (neurons with samples 72, 73, and 75 in Figure 2h; samples 73, 74, and 76 in Figure 6), located in the eastern part of the study area, which is characterized by low water exchange. Here, the delayed hydrodynamic regime leads to the retention and accumulation of nitrate–chloride components, effectively depleting the natural carbonate buffering capacity. The bifurcation point above the alkaline core identified on the SOM plane for pH, along with the topological discontinuity of the primary vector PC1 on the SOM planes for f₁–f₅, mark the physical boundary where dominant hydrodynamic process shifts from regional northern infiltration to local stagnant accumulation, driven by PC2. Since f₅ has the maximum positive loading on PC2, the nitrate–chloride complex serves as a key factor of local dispersion, and f₅ itself serves as the main indicator of anthropogenic metamorphism of groundwater in the eastern and south-eastern part of the territory. Thus, Cluster 2 is identified as a dynamic zone of technogenic facies formation (or of high salinity water type), where the PC1 projection discontinuity is driven by the system's loss of self-regulation capacity. Cluster 3 contains groundwater samples from monitoring well P-7 (intensive water exchange zone) and dug wells in the central and south-western areas (Figure 1). Spatial analysis (Figure 7) records TDS and pH stabilization at moderate values, with natural processes dominating. The intersection of the PC2 projection with Th_norm visually captures a critical hydrogeochemical threshold. The left segment of the PC2 projection in Cluster 2 characterizes low water exchange and contaminant stagnation in the eastern area. Conversely, the shift of the other PC2 projection segment into Cluster 3 indicates a changed migration vector, marking intensive water exchange in the central, western, and south-western parts. Here, active drainage and hydrodynamic dilution (e.g., monitoring well P-7) successfully reduce the contaminant load, restoring equilibrium via the influx of natural waters enriched with bicarbonates and$N{a}^{+}+K^{+}$.

In Cluster 3, the specific TDS dynamics for well P-7 are driven by extreme technogenic pollution (samples 84, 85, and 87 in Figure 7; neurons with samples 83, 84, and 86 in Figure 2h). Despite the system's tendency to restore natural equilibrium, high environmental aggressiveness cannot be rapidly neutralized owing to kinetic limitations of the system. Thus, Cluster 3 is identified as a zone of stable bicarbonate–sodium–calcium water formation (or relatively fresh water), where parameter f₅ serves as the vector of systemic evolution.

4. Conclusions

Based on the six log-ratio relative variables comprising hydrogeochemical data from 91 groundwater samples, the SOM was combined with PCA, and GCV to perform clustering of the hydrogeochemical data and assess the groundwater quality and factors governing groundwater geochemistry in the Pleistocene unconfined aquifer of the Southern Bug and Sinyukha interfluve area, Ukraine. Hydrogeochemical processes were analyzed by superimposing theoretical balance isolines (ln(Y/X) = 0) and principal component (PC1, PC2) projections onto the nonlinear SOM plane. This approach yielded the following key findings:

1. The synthesis of topological mapping and spatial relative series analysis allowed the identification of three clusters. Cluster 1 revealed salinity water types forming under the synergistic influence of sulfate mineral dissolution, ion exchange, and anthropogenic sulfate inputs. The latter exerts a dominant influence on groundwater in the northern part of the study area, which is characterized by intense vertical infiltration and constitutes the primary pollution zone. Cluster 2 identified high-salinity water types and a shift in dominant influence from sulfate pollution to intense nitrate and chloride pollution, indicating technogenic groundwater metamorphism in the eastern and southeastern parts of the study area. Maximum pollution is observed in the eastern part of the region, which physically corresponds to the zone of delayed water exchange. Cluster 3 revealed relatively fresh water types governed by silicate weathering and natural carbonate dissolution. The latter plays a leading role in stabilizing and restoring the natural equilibrium of the hydrogeochemical system in the central, western, and southwestern parts of the study area.

2. The proposed approach allowed us to establish cluster boundaries with high accuracy, track the dynamics of groundwater formation factors, and identify bifurcation points with phase transition zones on the SOM plane.

3. Cross-validation analysis revealed the intensity and trends of the hydrogeochemical processes and their impact on the TDS and pH variations, confirming the robustness of the boundaries for Clusters 1–3.

4. Although SOMs are traditionally classified as unsupervised machine learning techniques, this study successfully transformed the SOM framework into an XAI model. The proposed approach provided a robust methodology to verify cluster boundaries, offering a transparent identification mechanism that exceeds the interpretability of standard “black box” methods.