Using Random Forest Regression to Model the Spatial Distribution of Concentrations of Selected Metals in Groundwater in Forested Areas of the Wielkopolska National Park, Poland

Abstract

Monitoring groundwater pollution is an important issue in terms of analyzing threats to protected, environmentally valuable areas. The topographical and environmental characteristics of a given area are often mentioned among the factors affecting the dynamics and chemistry of groundwater. In this study, the random forest regression (RFR) model was used to determine the spatial distribution of selected metals, such as aluminum, calcium, iron, potassium, magnesium, manganese, sodium, and zinc. In the role of indicators describing terrain variability, derivatives of the digital elevation model (DEM) were employed, with a spatial resolution of 5 m, describing the topography of the terrain on a local scale, such as, among others, slopes, the aspect and curvatures of slopes, the topographic position index, and the SAGA wetness index, as well as generalized values determined for each sampling point of the areas contributing their runoff. In addition, environmental parameters were taken into consideration: forest habitat types, the structure of soil cover, and the seasons when samples were collected. This study used samples collected from 15 wells located in forested areas of the Wielkopolska National Park on seven dates. The results obtained show that random forest can be used with very good results to model the spatial variability of the concentrations of aluminum, potassium, magnesium, manganese, and sodium in groundwater. However, in the case of calcium and zinc, no correlations were found between the adopted indicators describing the spatial variability of the area and their concentrations in groundwater. In addition, the degree of importance of each predictor was determined in order to rank their importance in modeling the concentration of each of the metals in groundwater. The summary ranking of predictors indicates that the strongest influence on the predicted concentration of metals in groundwater is exhibited by profile curvatures, planar curvatures, multiscale TPI, and then the habitat type of the forest. On the other hand, curvature classifications, soil composition, and seasonality exhibit the smallest generalized impact on the results of modeling.

1. Introduction

Groundwater is one of the most important resources of the natural environment, constituting the largest source of fresh water [1]. The quality of groundwater is influenced by a number of factors, both environmental and anthropogenic [2,3]. They can significantly limit the availability of this water, and consequently affect the surrounding environment. Nowadays, an ever-growing number of soluble chemicals resulting from urban and industrial activities and modern agricultural practices pose a threat [4]. At the same time, groundwater, being a component of the critical zone, affects the development of plants [5].

The chemical composition of groundwater in a given place depends on the soil and rock substrates, the time of contact with the soil or rock substrates, biological processes, and the mixing of water coming from the surface and from under the surface [6,7,8]. This affects the spatial variability of groundwater chemistry. Depending on the location in the relief, a greater role may be played by biogeochemical processes [9] or processes related to the accumulation of organic matter [10]. Research conducted in the United States has shown that geochemical processes occurring in the unsaturated zone can significantly affect the composition of groundwater [11]. Knowledge of the spatial variability of the occurrence of metals in groundwater, combined with their location in the landscape, is necessary to identify possible sources of pollution in order to develop appropriate management strategies aimed at remediating or regulating these pollutants [12].

Monitoring groundwater quality provides important information that can be used for activities aimed at maintaining ecological systems, including those in environmentally valuable areas. Assessing temporal and spatial changes in water quality parameters is of fundamental importance to controlling and preventing water pollution [13,14,15]. Traditional geostatistical and numerical models are typically very complex, expensive, time-consuming, and require extensive data [16]; they are typically unable to reflect the complex biological, chemical, and physical characteristics used to describe water quality [17,18]. The use of remote sensing data describing the terrain can help to more easily and quickly assess the spatial variability of groundwater contamination by metals. This is especially true for areas where access for systematic sampling is difficult or even impossible [19].

Machine learning (ML) algorithms are increasingly used to analyze the variability of water quality parameters. ML models use the inductive hypothesis to analyze and “learn the rules” based on data without relying on a specific set of equations. Haggerty et al. [20] identified four main groups of ML models used in modeling groundwater quality. The first two groups are partially supervised models and ensemble models. The group of models with unsupervised learning included self-organizing map models, grouping, and multiple frameworks. However, the most numerous group of models consists of supervised models, which include decision tree models and random forests, artificial neural networks, support vector machine models, adaptive neuro-fuzzy inference systems, deep learning, regression models, comparative studies, and optimization techniques.

The random forest (RF) model used in this study is a classification algorithm consisting of many decision trees. RF randomly selects observations and traits to build each individual tree and, consequently, create an uncorrelated forest of trees [21,22]. Due to its non-parametric nature, in the regression of random trees in random forest regression (RFR), the data do not have to meet the conditions for distribution; the data may also be collinear. The RFR method also proves its worth if a large number of predictors are used [23]. A random forest can be treated as a form of a “grey box” [24,25]. Although individual trees cannot be checked separately, we obtained some measures for the interpretation of results. One of these is the importance of the variable describing how much worse the prediction would be if the data for this predictor were randomly shifted [25]. Random forests have the additional advantage of dealing with “small n, large p” problems, a common scenario in studies on groundwater quality where observational data are scarce compared to the number of potentially influential variables [26].

Numerous studies indicate that RF exhibits better quality when modeling pollutants in groundwater compared to other methods. In the case of modeling nitrates, they allowed better results to be achieved compared to, e.g., MRL, kriging, SVM, naïve Bayes, a generalized additive model (GAM), and a boosted regression tree (BRT) on the scale of both a single well [27] and a continent [28]. RF was also used to determine the spatial distribution of nitrate concentrations on a large scale using only spatial predictors describing the environment [29].

The data used to teach the models have a specific spatial resolution, which may affect the quality of the results obtained from the model [30]. Depending on the size of the area subjected to analysis, the type of data used to describe the terrain and climate, and its availability, the spatial resolution of the data used for modeling may range from several dozen meters to a resolution expressed in dozens of kilometers [28,31,32,33].

The aim of this study was to determine the possibility of describing the spatial variability of selected metal concentrations in the groundwater within the forested areas of the Wielkopolska National Park using parameters describing the terrain topography determined on the basis of remote sensing data. The analysis took into account both parameters determined locally and the impact of areas neighboring the sampling sites. The final aim of this study was also to determine the possibility of using spatial data with a fairly high spatial resolution of 5 m in modeling the spatial variability of metal concentrations in groundwater using RFR. This is particularly important in valuable, protected areas where access for collecting water samples is difficult.

2. Materials and Methods

2.1. Study Site

This study used data describing concentrations of metals in groundwater samples collected in the Wielkopolska National Park, located in the western part of Poland (52.27° N, 16.79° E) (Figure 1). The area of the Park is 76 km². The Wielkopolska National Park is located in the immediate vicinity of the city of Poznań, which has over 500 thousand inhabitants. In total, the agglomeration of the city of Poznań has over 1 million inhabitants and it is one of the few areas in Poland with a continuous increase in the number of people living there. According to the Köppen–Geiger classification, the climate in this area exhibits characteristics of an oceanic climate (Cfb) characterized by warm summers and mild winters [34]. The average annual air temperature is 8.5 °C. The coldest month is January (−1.0 °C) and the warmest is July (18.1 °C). The average annual precipitation is 507 mm, which is one of the lowest in Poland. The highest precipitation is observed in July (76.0 mm) and the lowest in February (22.9 mm) [35]. The average annual groundwater table depths measured in observation wells for the period from 2015 to 2022 ranged from 0.88 m to 17.93 m [36].

2.2. Flowchart

Figure 2 shows a workflow diagram for an example of one metal. Part A, aimed at determining the topographic and environmental parameters of the catchment and then determining the area for which modeling of metal concentrations dissolved in groundwater is possible, is common to all metals and was performed only once. The remaining part, covering the construction of the RFR model and the production of maps of the spatial distribution of metal concentrations in groundwater, was performed for each of the metals.

2.3. Groundwater Samples

The results of the analyses of Al, Ca, Fe, K, Mg, Mn, Na, and Zn concentrations were made available by the Wielkopolska National Park. The data used in the analysis are the only data regarding metal concentrations in groundwater in the Wielkopolska National Park. Water samples were taken from shallow groundwater, the depth of which ranged from 0.5 to 4.6 m. Groundwater samples were collected at 15 sites located in the forested areas of this national park on 7 dates: March, July, and November 2016, and February, May, August, and November 2017. This gives a total of 105 samples for each of the metals analyzed. These sites are marked with numbers 1, 2, 4, 5, and 6–17 (Figure 1). Water samples were collected after pumping water out of the well at least three times. Concentrations of elements dissolved in water were determined using atomic absorption spectrometry in certified laboratories in accordance with the W-METAXFL1 methodology.

2.4. Terrain Topography

In order to analyze the topography of the terrain, the author’s original digital elevation model (DEM) was employed, with a resolution of 5 m, covering the entire National Park area (Figure 1). The results of the LiDAR flight pass, made available by the Wielkopolska National Park in the form of “las” files, were used as source data. When constructing the DEM and determining the parameters of the terrain, the SAGA GIS 9.3.1 program was used [37] (SAGA User Group Association, Hamburg Germany).

The topography of the terrain influences groundwater recharge from water infiltrating from the surface but also determines groundwater flow patterns [38,39] and the movement of substances dissolved in it [40,41]. In the case of using topographic parameters obtained from remote sensing sources to describe a phenomenon, many parameters can be used [42,43]. It should be remembered that the determination of topographic parameters using DEM is subject to uncertainties and errors contained in the data and in the algorithms used to determine these parameters [44,45,46]. Therefore, the use of multiparametric models can reduce the impact of a single parameter error on the final result.

The topographic parameters used in this work allow for the assessment of water movement conditions near the points where water samples are collected for analysis. The description of the terrain topography for each of the raster cells included terrain slopes, aspect, general curvature, profile and planar curvature, curvature classification, the topographic position index (TPI), TPI-based landform classification, multiscale TPI, and the SAGA wetness index.

In order to describe the water sample collection sites, locally defined indicators were used, such as profile curvature, planar curvature, general curvature, curvature classification, the topographic position index (TPI), TPI-based landform classification, multiscale TPI, and the SAGA wetness index. The slopes and the aspect of the slope are among the basic factors describing the terrain [47]. Consideration of slopes is important for assessing the rate of infiltration. Areas with low slopes are characterized by high infiltration and low surface runoff. Areas with high slopes show the opposite relationship [48]. In turn, the slope’s exposure can be used to assess its insolation, which can affect, among other things, evapotranspiration. On southern slopes, stronger insolation leads to stronger water uptake by vegetation, as well as from the more heated soil in capillary rise.

Profile curvature measures the topographic curvature along a flow line, i.e., the steepest descent path. Planar curvature measures the curvature of contour lines on topographic maps. It is directly related to the convergence and divergence of flow lines [49]. Terrain slopes, the general profile, and planar curvatures were determined using the 9-parameter model presented by Zevenberg and Thorn [50]. These indicators allowed for determining the local humidity conditions.

The TPI index and the derived multiscale TPI and TPI landform classification indicators allowed for determining the location of cells in the relief of the terrain at a different level of generalization [51,52]. TPI compares the elevation of each cell with the average elevation of cells located in a specific neighborhood. Positive TPI values indicate cells that are positioned higher than the average for the neighborhood, which may indicate hill ridges. Negative TPI values correspond to locations located below the neighborhood (valleys). TPI values close to zero indicate flat areas or areas characterized by a constant slope of the terrain. TPI landform classification, on the other hand, additionally uses TPI standard deviation values to delineate discrete location classes on the slope [51]. A description of terrain forms was also carried out using a fuzzy classification based on terrain slopes and a number of different descriptions of terrain curvature [53].

The SAGA wetness index (SWI), a modification of the Topographic wetness index, was used to describe water flow conditions. This index employs a multi-directional water flow algorithm and iteratively modifies the contributing area based on the local slope of the terrain and the accumulation of flow in cells adjacent to the cell being analyzed, in accordance with the formula [54,55,56]:

$$SWI=ln\left({\displaystyle \frac{{SCA}_M}{\tan \beta }}\right)$$

$${SCA}_M={SCA}_{max}{\left({\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$t$}\right.}\right)}^{\beta exp\left(t^{\beta }\right)} \mathrm{f}\mathrm{o}\mathrm{r} SCA<{SCA}_{max}{\left({\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$t$}\right.}\right)}^{\beta exp\left(t^{\beta }\right)}$$

where: SCA_M is a modified specific catchment area of each DEM grid cell (m), SCA_max is the maximum stable value for SCA_M obtained through the iterations, β is the slope angle (arcs), and t is described as a suction factor. A high SWI value indicates that the region has a greater potential for soil water saturation [44].

The chemical pattern of groundwater is related to groundwater flow, and hydrochemical changes are controlled primarily by geochemical processes. The chemistry of shallow local groundwater may change systematically along the flow path [57]. In order to determine the possible impact of the areas surrounding sample collection sites on the concentrations of metals in the water, areas contributing their surface runoff were delineated for points 1, 3, and 8–17. These areas were determined using the Upslope Area tool with the multiple flow direction method. For points 4, 5, and 7, located in alluvial areas, a buffer zone with a radius of 50 m surrounding each point was adopted as the area that may influence the concentrations of metals in groundwater. The descriptions of the areas defined in this way included determining each contributing area’s maximum and average slopes, the standard deviation of slopes for the contributing area, and the average aspect. The standard deviation of the slope and mean gradient can be used to describe the conditions for runoff formation and infiltration in the contributing area.

2.5. Environmental Parameters

The distributions of forest habitat types and topsoil cover were made available by the Wielkopolska National Park in the form of shp files [58]. The data were rasterized with a resolution of 5 m, in a grid compatible with the DEM raster grid. Then, the dominant type of FHT and soil cover was determined for the previously described contributing areas of each of the points. In Poland, FHT is the basic unit in the forest habitat classification system. It includes forest areas with similar site conditions resulting from soil moisture and fertility, terrain shape, and geological structure [59,60]. In the case of shallow groundwater levels, forests can affect the water–salt balance of soils by increasing salt concentrations [61]. Different tree species affect the salinity of groundwater in different ways [62].

In addition, the analysis took into account the dates of measurements of metal concentrations in groundwater, classifying the months of sampling as spring: March, April, and May; summer: June, July, and August; autumn: September, October, and November, and winter: December, January, and February. This may allow the assessment of the seasonal variability of metal concentrations in groundwater.

2.6. Random Forest Regression

All statistical analyses were performed using the program R, version 4.4 [63] (R Foundation for Statistical Computing, Vienna, Austria).

The random forest regression (RFR) method was used to determine the relationships between the concentrations of metals in groundwater and the geomorphometric and environmental parameters described above. RFR is a machine learning method that allows for classification and regression, consisting of the construction of many decision trees [22]. The philosophy behind the learning techniques used in RFR assumes that its accuracy is higher than that of other machine learning algorithms because the combination of predictions works more accurately than any single component model [28]. RFR can be used to model the dependencies of non-linear variables. It allows for determining complex relationships between variables and is not affected by the collinearity of variables [64]. RFR requires two defined parameters—the number of factors to be used in each tree-building process (mtry) and the number of trees created in the forest (ntree).

The analysis employed the “randomForest” package in version 4.7-1.1 [22], which is part of the R program.

The measurement dataset for each metal was 105 measurements. For each of the metals, a subset of 75% of the data (i.e., 79 items of measurement data) were randomly separated from each dataset and was then used to teach the models. The remaining 25% of the data (i.e., 26 items of measurement data) constituted a test dataset that was used to assess the quality of the models. The accuracy of the model was then assessed by comparing the predicted value with the actual values of the test data using the Pearson correlation coefficient.

RF modeling allows for testing the significance of all predictive variables in relation to the response variable. The increase in the purity of nodes from the division of decision trees based on the analyzed variable (IncNodePurity) is used to measure the importance of variables [65]. This metric of the importance of variables was used to determine the overall validity of environmental variables in the model.

3. Results and Discussion

3.1. Topography of the Area

Figure 3, Figure 4 and Figure 5 show the topographical and environmental parameters of the WNP area. These values were used to determine the conditions describing the sites of water sampling, as presented in

Table 1. Topography of sampling points.

Sampling Point	General Curvature	Profile Curvature	Planar Curvature	TPI	Multiscale TPI	TPI Landform Class	Profile Curv. Class	Tangential Curv. Class
1	0.0139	0.0093	−0.0155	−3.92	−2.24	2	1	2
2	−0.0036	−0.0021	0.0042	−1.09	−0.74	4	0	1
4	0.0031	0.0009	0.0363	−0.33	0.29	5	2	2
6	−0.003	0.0016	−0.0406	−0.15	−0.15	4	0	1
7	0.0014	0.0008	−0.0155	−0.32	0.1	3	1	2
8	0.0043	0.0004	0.0567	0.41	0.49	5	2	2
9	−0.0076	−0.0058	0.0324	−1.37	−0.95	1	0	1
10	0.0074	0.0024	0.0197	−2.76	0.77	1	2	2
11	−0.0014	−0.0006	−0.0002	1.64	1.79	8	1	2
12	0.0076	0.0007	0.0284	−3.81	0.6	0	2	2
13	−0.0089	0.0038	−0.073	−0.53	−2.38	3	2	0
14	−0.0035	−0.0028	0.0126	0.22	−0.18	4	0	2
15	0.0015	−0.0006	0.0422	−1.24	−0.05	1	2	2
16	−0.0039	−0.0015	−0.0209	−1.83	−0.84	1	0	2
17	−0.0096	−0.0048	−0.0062	−2.07	−3.19	1	0	1

and

Table 2. Topography of point catchments.

Sampling Point	Area (m2)	Max Slope (%)	Mean Slope (%)	Slope SD (%)	Average Aspect (°)	Mean SAGA WI	Forest Habitat Type	Soil Type
1	5350	37.5	16.4	7.7	121	5.25	2	1
2	14,125	39.5	15.3	7.8	196	5.85	2	3
4	8850	18.7	4.9	3.7	178	8.5	1	1
6	8850	11.3	4.1	2.3	172	8.43	5	5
7	8850	19.6	3.7	4.2	223	9.39	5	2
8	57,275	5	1.6	0.7	73	11.19	3	2
9	6550	18.5	6.8	4.6	188	7.37	4	3
10	6850	41.5	13.1	10.8	143	7.24	2	2
11	5400	18.6	6.4	4.6	179	7.65	1	1
12	47,400	51	10.2	10.3	130	8.17	3	2
13	2875	19.9	8.5	4.6	208	6.39	3	6
14	12,650	14.4	4.2	2.8	265	8.29	1	1
15	47,025	23.8	7.2	5.2	243	8.13	2	2
16	8375	48.5	25.9	11.6	136	5.32	1	1
17	30,950	90	20.7	19.8	149	5.99	3	2

SD—standard deviation; FHT: 1—fresh mixed coniferous forest, 2—fresh mixed broadleaved forest, 3—fresh broadleaved forest, 4—moist broadleaved forest, 5—riparian forest; soil: 1—loose sand, 2—light clay sand, 3—loamy sand, 5—sandy silt, 6—sandy loam.

. The smallest topographic catchment area, covering 2875 m², belongs to point 13. However, the largest area, amounting to 57,275 m², is characteristic of point 8. The average slope for the catchment area ranges from 1.6% for catchment area 8 to 25.9% for catchment area 16. The highest variability of slopes is shown by the catchment area in point 17, where the standard deviation is 19.8%, and the smallest is in point 8 with 0.7%.

The top layers of the soil cover of the analyzed areas are dominated by formations characterized by excellent water permeability. Most of the soil cover of the catchment area in the measurement points consisted of slightly loamy sands (six catchments) and loose sands (five catchments). The rest predominantly consisted of loamy sands, silty sands, dust, and sandy loams (Table 2).

Due to the habitat types of the forest, the area of four catchments is dominated by fresh broadleaf forests, fresh mixed forests, and fresh coniferous forests. In the case of three catchments, riparian forests dominate, and in one catchment, moist broadleaf forests are the most prevalent.

Regarding the points where the sites of water sampling were located, they were situated in various landforms. The wells were located both on concave slopes (negative curvature values) and convex slopes (positive curvature values), as well as on flat land (≈0 values).

3.2. Description of Metal Concentrations in Groundwater

Table 3. The statistical description of metal concentrations in the groundwater (mg dm⁻³).

	Mean	Median	Min	Max	SD	Skewness	CV
Al	0.034	0.021	0.005	0.208	0.043	2.77	126
Ca	133	129	52.6	260	47.7	0.69	36
Fe	0.543	0.023	0.001	7.2	1.14	3.34	210
K	3.8	2.52	0.664	23.2	3.09	1.84	81
Mg	17.5	17.3	7.95	38.1	6.41	0.69	37
Mn	0.31	0.19	0.002	1.61	0.36	1.85	116
Na	29	23.2	6.9	97.4	19.7	1.3	68
Zn	0.0063	0.0033	0.001	0.093	0.012	5.4	190

SD—standard deviation; CV—coefficient of variation.

presents descriptive statistics of the concentrations of eight metals in groundwater. The average concentrations of Al, Ca, Fe, K, Mg, Mn, Na, and Zn in groundwater were 0.034, 133, 0.543, 3.80, 17.5, 0.31, 29.0, and 0.0063 mg dm⁻³, respectively. These values arrange the concentrations of metals in the following series: Ca > Na > Mg > K > Fe > Mn > Al > Zn. Analyzing the values of the CV coefficient of variation, one can notice that Fe and Ca concentrations exhibit the greatest variability, for which the CV is, respectively, 210% and 190%. Ca and Mg concentrations, however, are characterized by the lowest variability, for which the CV values are 36% and 37%, respectively.

3.3. RFR Modeling Results

The models that were created based on the learning data were assessed by comparing the concentrations of metals in groundwater, which were calculated for the test data with the measured concentrations (Figure 6). The best compliance of the results calculated by using the model with the measured results was obtained for manganese (Mn), for which the Pearson correlation coefficient was r = 0.96. Slightly lower levels of correlation occurred for potassium K (r = 0.91), aluminum Al (r = 0.86), magnesium Mg (r = 0.84), and sodium Na (r = 0.82). For these compounds, the level of significance is p < 0.001. Lower compliance of the concentrations obtained from the model with the measured concentrations was obtained for Fe/iron. The correlation coefficient of r = 0.44, with a significance level of p = 0.03, indicates that the results obtained are fairly low-quality. On the other hand, the quality of RFR models generated for calcium and zinc does not allow their use in modeling spatial variability. It should be emphasized that in the analysis, all topographic and environmental parameters were determined for a resolution of 5 m. The research conducted by Wu et al. [31] indicates that the optimal resolution of determining topographic parameters for determining the concentrations of heavy metals in the soil may vary not only for individual parameters but also across different metals. For example, the optimal spatial resolution determined by them in the case of a terrain aspect ranged from 90 m for Zn to 3000 m in the case of Cu and Pb. On the other hand, in the case of zinc, the optimal resolution varied from 90 m for aspect to 250 and 1500 m for the depth of the valley, 1500 m in the case of TWI, and 3000 m for the slopes of the terrain.

In addition, the quality of the results obtained in the modeling process may be influenced by both the number of variables and the nature of the terrain features they describe. The study focuses on 15 parameters describing the topography of the terrain, and the remaining features include the top layers of the soil cover, the date of sampling, and the habitat types of the forest. These features provide a very good description of the spatial distributions of K, Al, Mg, and Na. The improvement of results for Fe, or the possibility of modeling Ca and Zn, may be influenced by supplementing the range of predictors with variables describing geology, soil properties, or groundwater table depths [66]. However, increasing the range of predictors may be challenging due to the lack of detailed information on their spatial variability. The effect of groundwater depth on dissolved metal concentrations was not considered in this study because modeling the spatial distribution of concentrations, which is the final output of this work, would require knowledge of the groundwater depth for each raster cell.

3.4. Importance of Predictors

The random forests algorithm allows for determining the contribution of individual parameters describing the terrain to the results obtained from the model. As shown in Figure 7, for Ca and Zn, it was not possible to identify topographic parameters that would allow the spatial variability of these metals to be modeled. On the other hand, for other metals, the variability of concentrations is dependent to varying degrees on various parameters. In the case of Al, the greatest impact on the concentration is caused by the profile curvature and the maximum slope in the watershed. The profile curvature influences the erosion and sedimentation processes taking place on the slope [67] and also determines the moisture conditions of the topsoil. Studies conducted in Brazil have shown lower aluminum contents on concave slopes than on convex ones [68].

In the case of Fe, the influencing factors are seasonality and the average slope of the contributing area. Strong seasonal variability of Fe concentrations in groundwater, unlike other metals, was also observed in earlier studies [69]. This variability can be explained by the influence of the high variability of redox dynamics in soils related to oxygen availability and changes in groundwater depth, as well as by seasonal changes in anaerobic microbial activity [70,71].

The second metal whose concentration in groundwater is directly related to redox changes is manganese. Depending on the oxidation state, its solubility changes, and consequently its concentration in groundwater is also altered [72]. However, analyses of data available for the WPN area indicate that for Mn, the most important parameters are the multiscale TPI and TPI.

For potassium, these factors are the general curvature, the SAGA WI, and the TPI landform. In the case of aluminum, the studies conducted by dos Santos [68] showed a relationship between K content and terrain curvature. However, in this case, they observed the opposite relationship, i.e., a higher potassium content occurred on concave slopes.

In the case of Mg, the strongest impact is caused by the planar curvature and the multiscale TPI. In the case of sodium, the average slope of the catchment and the fuzzy landform are key (Figure 7). These parameters indicate the influence of parameters describing water infiltration not only in the immediate vicinity of the sampling site but also include the influence of the contribution area. It is indicated that sodium concentrations are strongly influenced by the processes of percolation and dissolution of salt contained in the soil [73].

As presented in Figure 8, the ranking of predictors indicates that the strongest influence on the predicted concentrations of metals in groundwater is exhibited by profile curvatures, planar curvature, the multiscale TPI, and then the forest habitat type (FHT). It should be noted that the FHT is most commonly associated with local humidity conditions. Di Stefano et al. [67] indicate that the relationship between different forms of terrain curvature and soil chemical properties, and consequently groundwater chemistry, can be linked to the dynamics of erosion and sedimentation processes and the distribution of vegetation in a given area. Conversely, curvature classifications, soil composition, and the season of sampling have the least generalized impact on the modeling results. The low impact of soils may result from the low diversity of the surface layers of the soil cover of the analyzed areas; in 11 out of the 15 places where groundwater sampling occurred, highly permeable materials (loose sands and loamy sands) predominated. In addition, studies conducted by Kendall et al. [74] and Walker et al. [75] confirm a much greater spatial variability of groundwater chemistry compared to time variability.

3.5. Result Maps

The models obtained during the analysis, which depict the relationships between metal concentrations and indicators describing the terrain, were used to create maps of the spatial variability of Al, Fe, K, Mg, Mn, and Na concentrations for the forested areas of the Wielkopolska National Park (Figure A1 and Figure A2). The maps were created at a resolution of 5 m, for areas for which the variability range of the predictors was consistent with the ranges used during random forest model development, as shown on Figure 2. This is due to the limitations of RF when data outside the range used for training are used for modeling. The results obtained in the extrapolation process may be completely wrong in such cases [76]. The maps allow for the assessment of the spatial distribution of metal concentrations in groundwater in areas where very detailed field studies are very difficult. This applies, for example, to very valuable natural protected areas in a national park, access to which may be difficult or impossible due to legal restrictions.

4. Conclusions

This study assesses the possibility of using random forest regression (RFR) techniques to model the spatial variability of concentrations of aluminum, calcium, iron, potassium, magnesium, manganese, sodium, and zinc in groundwater. The results of the analysis of data from samples collected in forested areas of the Wielkopolska National Park indicate that for the adopted data resolution of 5 m, RFR allows for modeling the variability of Al, Fe, K, Mg, Mn, and Na concentrations in the groundwater. However, the model did not detect a relationship between the adopted predictors and Ca and Zn concentrations. Among the 18 indicators describing the variability of the terrain morphology and environmental variables, the most important factors in the creation of models are local parameters: profile and planar curvature (expressed in absolute values) and the multiscale TPI, as well as the FHT. In contrast, curvature classifications, soil composition, and sampling season have the least impact.

It should be noted that all sampling points are located in forested, protected areas of the National Park, which minimizes the risk of the occurrence of point sources of metal contamination. The results obtained indicate that modeling the spatial variability of metal concentrations in groundwater provides information about the state of pollution of these waters in the studied areas. The maps of the spatial distribution of metal concentrations in groundwater prepared as the final result of this work can be a valuable source of information, allowing protective measures to be planned in the Wielkopolska National Park.