The Impact of Catchment Characteristics and Weather Conditions on Heavy Metal Concentrations in Stormwater—Data Mining Approach

: The dynamics of processes a ﬀ ecting the quality of stormwater removed through drainage systems are highly complicated. Relatively little information is available on predicting the impact of catchment characteristics and weather conditions on stormwater heavy metal (HM). This paper reports research results concerning the concentrations of selected HM (Ni, Cu, Cr, Zn, Pb and Cd) in stormwater removed through drainage system from three catchments located in the city of Kielce, Poland. Statistical models for predicting concentrations of HM in stormwater were developed based on measurement results, with the use of artiﬁcial neural network (ANN) method (multi-layer perceptron). Analyses conducted for the study demonstrated that it is possible to use simple variables to characterise catchment and weather conditions. Simulation results showed that for Ni, Cu, Cr, Zn and Pb, the selected independent variables ensure satisfactory predictive capacities of the models (R 2 > 0.78). The models o ﬀ er considerable application potential in the area of development plans, and they also account for environmental aspects as stormwater and snowmelt water quality a ﬀ ects receiving waters.


Introduction
Urbanisation and increasingly frequent extreme weather events result in higher volumes of stormwater flowing from catchments. That, together with the deterioration of air quality in urban areas, leads to growing amounts of pollutants that are washed off surfaces in rainfall events. Consequently, they negatively affect the aquatic environment of the receiving waters. Various mathematical models are currently developed in order to predict negative stormwater impacts on the natural environment [1][2][3]. They are useful when making decisions, in adherence to sustainable development principles, which concern land use and urban space management. To predict stormwater quality, these models primarily take into account deposition of pollutants on surfaces [1], their runoff, and transportation by the drainage system. The Stormwater Management Model (SWMM) software is commonly used; however, simplifications [2] adopted in the model can make it difficult to fit simulation results to measurements. Literature overview [4] shows that, in many cases, unsatisfactory prediction results can be caused by the simplified model of pollutant deposition. The complex process of pollutant accumulation on the catchment surface is influenced by numerous factors, including air quality, weather patterns, and local conditions. Air quality and pollutant deposition conditions vary in time [5,6]; therefore, their proper modelling is not easy and requires advanced tools. Consequently, statistical models are more frequently used to account for the stochastic character of determinants that affect stormwater

Characteristics of the Study Area
The investigations into stormwater quality and quantity were conducted for three urban catchments in the city of Kielce, Poland ( Figure 1). In the catchments, which differ in land use, five characteristic surface types with respect to runoff were identified (Table 1): asphalt and gravel road surfaces, roofs, car parks, and green spaces. The first catchment, having a total area of 62 ha, is located in the central-eastern part of the city. It holds main transportation routes, service sector areas, and high-rise residential buildings. Overall, impervious areas with a high runoff coefficient (>0.8) constitute 51.5% of the total area of the catchment, which indicates its typical urban character. The highest point of the catchment is at 271.2 m above sea level and the lowest at 260 m a.s.l. The average slope of the land is 0.71% [15]. The stormwater from the catchment area is collected by the existing stormwater drainage system and conveyed by a main sewer to the stormwater treatment plant (SWTP) in IX Wieków Kielc Street (Figure 1a).
The Witosa SWTP catchment, which has an area of 83 ha, is located at the outskirts of the city. On the east, and partially north side, the catchment is surrounded by an open ditch collecting stormwater flowing from a dense forest area. The ditch turns into a ø 800 mm closed sewer, which is connected to a sewer conveying effluent from the Witosa treatment plant (Figure 1b). All the stormwater is piped (ø 1400 mm) to the receiving water of the river Silnica. Land development includes primarily one-and multi-family housing, and the share of impervious surfaces is 25.9% of total catchment area Appl. Sci. 2019, 9, Table 1). The highest point in the catchment is 365.5 m a.s.l., and the lowest is 291.25 m a.s.l. The average slope of the land is 8.2% [16]. Jesionowa 400 11.3 8.4 11.5 11.2 57.6 Jesionowa SWTP catchment occupies the largest area (400 ha). It is located in north-western part of Kielce and includes highly urbanised areas. The catchment land use is dominated by industrial zones with large commercial buildings, and low residential buildings. As a result, the share of impervious surfaces is 42.4% (169.6 ha). The highest point in the catchment is 315 m a.s.l., and the lowest is 265 m a.s.l. The average slope of the land is 2.65%.
In all three catchments, a separate sewage system is used, in which drains collect stormwater and sanitary wastewater is carried by sewers.

Measurement Apparatus
In rainfall events, the flow was measured, and at the same time, samples were collected for analytical tests. The measurement apparatus was installed in the separation chambers (monitoring  Jesionowa SWTP catchment occupies the largest area (400 ha). It is located in north-western part of Kielce and includes highly urbanised areas. The catchment land use is dominated by industrial zones with large commercial buildings, and low residential buildings. As a result, the share of impervious surfaces is 42.4% (169.6 ha). The highest point in the catchment is 315 m a.s.l., and the lowest is 265 m a.s.l. The average slope of the land is 2.65%.
Appl. Sci. 2019, 9, 2210 4 of 16 In all three catchments, a separate sewage system is used, in which drains collect stormwater and sanitary wastewater is carried by sewers.

Measurement Apparatus
In rainfall events, the flow was measured, and at the same time, samples were collected for analytical tests. The measurement apparatus was installed in the separation chambers (monitoring point- Figure 1) of all treatment plants. The ISCO AV 2150 module flow meters (Teledyne ISCO, Lincoln, NE, USA) were employed to measure stormwater quantity. Stormwater samples were collected using ISCO 6712 portable samplers (Teledyne ISCO, Lincoln, NE, USA). Unstabilised samples were immediately transported to a chemical laboratory in order to determine the following quality indicators, i.e., HM: Ni, Cu, Cr, Zn, Pb, and Cd. The pH value was measured in accordance with the PN-EN ISO 10523:2012 method using SevenMulti™ meter (Mettler Toledo, Columbus, OH, USA) [17]. Samples were digested in nitric acid using microwave oven (Multiwave 3000, Anton Paar, Graz, Austria) and filtered with membrane filters (0.45 µm). Concentrations of HM were determined by atomic emission spectrometry with inductively coupled plasma ICP Optima 8000 (Perkin Elmer, Waltham, MA, USA) with certified multi element standards (Perkin Elmer) [18]. The tipping bucket type rain gauge RG50 from SEBA Hydrometrie GmbH (Kaufbeuren, Germany) was used to measure rainfall depth. Measurement frequency ranged from 2 to 5 min, with a resolution of 0.1 mm.

Rainfall Events
The analyses were based on the results of rainfall depth measurements recorded since 2009. The measurements were taken using a rain gauge located in Kielce (Figure 1), within the Jesionowa SWTP catchment. The rain gauge is situated 1.2 km from the boundaries of the Witosa SWTP catchment and 2.5 km from the IX Wieków SWTP catchment. For calculation purposes, individual precipitation events were extracted from rainfall data. Based on the analysis of stormwater drainage system operation, the following parameters were adopted in the study-minimum inter-event (rainless) time (t bd ): 4h (DWA A-118 [19]) and minimum rainfall depth (P tot ): 2.0 mm. Although average weighted catchment retention is higher and amounts to 3.81 mm [20], rainfall with smaller P tot value may induce runoff, which is determined by specific rainfall distribution within catchment area. The number of precipitation events ranged yearly from 52 to 75. Precipitation events were parameterised with respect to total rainfall depth and rainfall duration (t r ). Values t r for the observed precipitation events varied in the range 10-2366 min., t bd was 0.16-60 days, and P tot 2.0-45.2 mm.

Meteorological Conditions and Parameters of Hydrographs
In the city of Kielce, monitoring of air quality and weather conditions was conducted in the years 2009-2018. Daily measurements included visibility, temperature and wind speed, and the occurrence of rainy, snowy and foggy days was noted. In this study, visibility is the measure of air pollution level [5]. The measurements demonstrated that in the years 2008-2018, selected parameters ranged as follows: visibility 8.7 to 33.0 km, average air temperatures from −1.3 • C to 22.0 • C, average wind speed was from 0.6 to 11.4 m/s, annual number of rainy days ranged from 189 to 201 and snowy days ranged from 36 to 78. In addition, it was found that in the periods between individual precipitation events, the number of days with rainfall ranged from 60 to 75, those with snowfall in the range of 1-70 and with fog 25-71.
The hydrographs analysed for the study showed different patterns, depending on the maximum flow values and runoff duration in the sewer. The highest flow (Q) values were observed in the Jesionowa SWTP (maximum 4.53 m 3 ·s −1 ), which results not only from precipitation characteristics, but also from the catchment size (Table 1). For Witosa SWTP and IX Wieków Kielc SWTP, Q values amounted to slightly over 0.5 m 3 ·s −1 . Hydrographs duration values ranged from 60 to 850 min.

Statistical Analysis
Statistical analyses were employed to determine mean, median and range values. In order to identify the relationships between variables under consideration (HM, percentage of each land use type in a given catchment, and weather characteristics), Pearson correlation coefficient and regression analysis were used. Additionally, it was checked beforehand whether empirical distributions of the independent variables of concern do not deviate from normal distribution. The Kolmogorov-Smirnov test was used for that purpose. When the determined value for individual independent variables was p ≤ 0.05, no grounds were found to reject the assumption that the empirical distribution of variables was normal. Otherwise (p > 0.05), the data were subjected to Box-Cox transformation [21]. Initially, to develop the ANN model for predicting HM concentrations, the independent variables describing catchment characteristics, precipitation and weather conditions [1,3,10,22] were adopted. Prior to further analyses, those data were normalised and standardised [23]. Then, the independent variables that have negligible effect on the simulation results were eliminated from the set of potential variables using the Fisher-Snedecor test. The next stage involved the development of ANNs models-the multi-layer perceptron (MLP), and structure optimisation (conjugate gradient, gradient descent, number of neurons in the hidden layer and activation functions in the hidden and output layers) [23]. Figure 2 shows the diagram of the adopted computational procedure [24].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 5 of 13 variables was p ≤ 0.05, no grounds were found to reject the assumption that the empirical distribution of variables was normal. Otherwise (p > 0.05), the data were subjected to Box-Cox transformation [21].

The Model for Stormwater Quality Prediction
Initially, to develop the ANN model for predicting HM concentrations, the independent variables describing catchment characteristics, precipitation and weather conditions [1,3,10,22] were adopted. Prior to further analyses, those data were normalised and standardised [23]. Then, the independent variables that have negligible effect on the simulation results were eliminated from the set of potential variables using the Fisher-Snedecor test. The next stage involved the development of ANNs models-the multi-layer perceptron (MLP), and structure optimisation (conjugate gradient, gradient descent, number of neurons in the hidden layer and activation functions in the hidden and output layers) [23]. Figure 2 shows the diagram of the adopted computational procedure [24].
where HMC (t)-heavy metal concentrations in time t, G-event type (rainfall, snowmelt), Acatchment area; Zi-land use; Q-stormwater flow rate in the SWTP gauging section of the analysed catchment; Pi, tr, tbd-variables describing rainfall characteristics, including Pi (total depth-Ptot, maximum ten minutes rainfall depth-P(10)); tr-rainfall duration, tbd-inter-event time; Vi, T, w, Tsn, Tf-variables characterising meteorological conditions and air quality based on Vi-visibility, T-air temperature, w-wind speed, Tsn,f-number of days with snowfall (or fog) preceding the event and ts -time measured from rainfall onset to the instant, for which stormwater quality is to be determined (precipitation events).

Selection of Variables for the Model
To predict HM concentrations in stormwater, the general formula was used: where HMC (t)-heavy metal concentrations in time t, G-event type (rainfall, snowmelt), A-catchment area; Z i -land use; Q-stormwater flow rate in the SWTP gauging section of the analysed catchment; P i , t r , t bd -variables describing rainfall characteristics, including P i (total depth-P tot , maximum ten minutes rainfall depth-P (10) ); t r -rainfall duration, t bd -inter-event time; Vi, T, w, T sn , T f -variables characterising meteorological conditions and air quality based on Vi-visibility, T-air temperature, w-wind speed, T sn,f -number of days with snowfall (or fog) preceding the event and t s -time measured from rainfall onset to the instant, for which stormwater quality is to be determined (precipitation events).
In order to thoroughly identify the physics of the analysed phenomena, the following parameters were studied for the inter-event (rainless) period (t bd ): average values, variance in the values (x i -independent variables found in Equation (1) Based on the measurements of the variables above x i (t 0 ) (where t 0 = 0 − t bd ), trends in variable changes were shown in the period of concern by determining the values of coefficients in the linear model (estimated with the least squares method) x i = a i t + b i (Figure 3).  The analyses demonstrated that the values of coefficient a = tg α range from −3.48 to 5.55, from −8.5 to 1.6, and from −2.85 to 1.4. These values indicate considerable variability of trends in changes in weather conditions during the analysed inter-event periods, which may have significant impact on pollutant deposition on the catchment surface.
For snowmelt events, t is measured from the minimum flow value Qmin (t = 0), which is the basis for identifying a specific event of given character. HMC(t) values are calculated for successive values of Q(t+i) > Qmin obtained from the measurements in snowmelt events.
Two hypotheses are formulated in the Fisher-Snedecor test: (i) the null hypothesis (H0), in which structural parameters are not significantly different from zero (α1 = α2 =…αk = 0) and (ii) the alternative hypothesis (HA), according to which at least one parameter is found that is significantly different from zero (α1 ≠ 0, α2 ≠ 0,…, αk ≠ 0).

Artificial Neural Networks
Artificial neural networks are frequently used for modelling different phenomena. Although many more complex and less complex ANN have been developed, one of the most popular is the unidirectional multi-layer network-the multi-layer perceptron (MLP) (Figure 4). The analyses demonstrated that the values of coefficient a = tg α range from −3.48 to 5.55, from −8.5 to 1.6, and from −2.85 to 1.4. These values indicate considerable variability of trends in changes in weather conditions during the analysed inter-event periods, which may have significant impact on pollutant deposition on the catchment surface.
For snowmelt events, t is measured from the minimum flow value Q min (t = 0), which is the basis for identifying a specific event of given character. HMC(t) values are calculated for successive values of Q(t+i) > Q min obtained from the measurements in snowmelt events.

Artificial Neural Networks
Artificial neural networks are frequently used for modelling different phenomena. Although many more complex and less complex ANN have been developed, one of the most popular is the unidirectional multi-layer network-the multi-layer perceptron (MLP) (Figure 4). different from zero (α1 ≠ 0, α2 ≠ 0,…, αk ≠ 0).

Artificial Neural Networks
Artificial neural networks are frequently used for modelling different phenomena. Although many more complex and less complex ANN have been developed, one of the most popular is the unidirectional multi-layer network-the multi-layer perceptron (MLP) (Figure 4). In ANN model, individual input signals (xi) that reach the input layer are multiplied by the values of weights (wij), then by products xi•wij. Next, they are transmitted to the hidden layer neurons, where they are summed up (Figure 4). The sums obtained undergo transformation by the activation In ANN model, individual input signals (x i ) that reach the input layer are multiplied by the values of weights (w ij ), then by products x i ·w ij . Next, they are transmitted to the hidden layer neurons, where they are summed up (Figure 4). The sums obtained undergo transformation by the activation function f (linear or nonlinear), and then they are transmitted to the output layer neurons. In the MLP model, the values of weights (w ij ) are estimated at the training stage. The estimation is intended to minimise the following function [23]: where E(RMSE)-target function, in the case under consideration equivalent to RMSE, n-number of data in a set, y s -value y obtained through iteration and y m -value y obtained through measurements. That is achieved with appropriate numerical methods, for instance, that of the Broyden-Goldfarb-Shanno.
The values of the dependent variable outputs (y) are determined using the following formula [23]: where I-number of model inputs, J-number of neurons in the hidden layer, w ij -the values of weights between inputs and neurons in the hidden layer, b j -loads of neurons in the hidden layer, w j -the values of weights between neurons in the hidden layer and the neuron in the output layer and f-the activation function. For the set of independent variables (x i ), the Fisher-Snedecor test was used to find the optimal structure of the MLP model for the number of the hidden layer neurons ranging (i−2)·(i+1) [24]. In the simulation tests, the following activation functions were taken into account: linear, exponential, tangent-hyperbolic and sine. The Broyden-Goldfarb-Shanno algorithm [23] was employed to determine the values of weights in the model. Due to substantial impact of the initial values of weights on the simulation results, and also problems with their optimisation, each ANN (having a specified number of neurons and activation function) was generated many times (5000 times) for different boundary conditions (initial values of weights for individual neurons). The above-mentioned algorithm was computed using the MATLAB software. The program automatically generates initial values [25], and the user can only choose the decay values for the hidden and output layers, which were assumed to be 0.001 and 0.0001 [8], respectively. The assessment of the predictive capacities of the models was based on coefficient determination (R 2 ), mean absolute error (MAE) and root-mean-square error (RMSE), as discussed in [26].
For the assumed values of neuron numbers and activation functions, the optimum values of weights were sought. The computations were performed using the method of successive approximations (for various numbers of iterations). In succession, for a number of neurons of the hidden layer, assumed a priori, activation functions were selected (hidden and output layers). With adopted assumptions, the computations of weight values were made using MATLAB software. Assuming successive values of neuron numbers to range from i to 2i+1 (1 neuron difference between successive models of concern), ANN models were determined for different combinations of the activation function. At the analysis stage, measures of fit of the computational results to measurements were obtained for the models produced. Structure and model considered optimal were those for which the lowest values of RMSE and MAE were found. The approach, described in details in [8], is frequently used when developing ANN models.

Heavy Metal Concentrations in Runoff
Heavy metals found primarily in the airborne particulate matter are particularly hazardous substances. Deposited on the ground together with particulate matter, HM are then washed off the catchment and transported, with the stormwater, into the wastewater receiving body. High concentrations of Zn, Pb and Ni in the stormwater result from the widespread use of these elements in the automotive and fuel industries. The presence of copper and chromium in the atmosphere is mainly caused by coal combustion and industrial activity [27][28][29][30][31].
Zinc had the highest concentration of all HM analysed ( Table 2). In stormwater samples collected from the Jesionowa SWTP catchment, Zn concentrations ranged 319-5731 µg·dm −3 , and were higher than the values observed for the IX Wieków Kielc SWTP (0-3873 µg·dm −3 ) and for Witosa SWTP (45-3160 µg·dm −3 ). These values are similar to the range specified for Bialystok, Poland (200-6000 µg·dm −3 ) [27], and for Belgrade (284-6200 µg·dm −3 ) [32]. The high concentrations of Zn in the studied catchments result from significant percentage (9.4%-14.3%) of the area of roof covers in land development [33], especially in the industrial part of the SWTP Jesionowa. The highest concentration of Pb (1405 µg·dm −3 ) was observed in the IX Wieków Kielc SWTP catchment, and it fits within the limits for runoff from motorways and major roads (3-2410 µg·dm −3 ) [34]. The high Pb values can be caused by a large volume of vehicular traffic in the city centre-exhaust gas, among others, is a source of Pb [28]. Conversely, the lowest Pb concentrations (1-343 µg·dm −3 ) were observed in the Witosa SWTP catchment, which may result from the prevalent single-family housing and local character of roads. However, Pb concentrations are higher than limits for Paris [35,36] and Genoa [37], although they fit within the range for runoff from suburban and residential roads (10- Cr concentration ranges in IX Wieków Kielc SWTP (0-350 µg·dm −3 ) and Witosa SWTP (3-319 µg·dm −3 ) catchments were similar, whereas the highest concentration was found for the Jesionowa SWTP catchment (2236 µg·dm −3 ). It was greater than that found in Serbia [32] (1350 mg dm −3 ) for asphalt surface samples collected from the car park of the University of Belgrade. So high concentrations of Cr can be associated with the industrial character of the catchment and carbon incineration in household boilers [28]. Additional source of Cr is a municipal heat and power plant located less than 1.0 km from the studied catchment.

Correlation Significance Analysis
The existence of statistically significant relations between stormwater HM concentrations and catchment characteristics and weather conditions may deliver information on HM sources (e.g., pollution from transport, atmospheric deposition, etc.) [43,44]. Additionally, strong correlation of individual HMs (Zn-Cu, Cr-Cu, Pb-Ni, Cd-Ni, Pb-Cd and Cr-Zn) may indicate common sources [45]. Analysis of data in Table 3 suggests that HM present in stormwater from the three catchments originate in either traffic, pollutant wash-off from impervious surfaces, or atmospheric deposition.
Pb, Ni and Cd concentrations in stormwater are positively correlated with the area of roofs (RF), roads (RD) and car parks (CP), which suggests that atmospheric deposition is not the main source of those elements. In the catchments of concern, Pb, Ni and Cd originate primarily from traffic and roof covers [46]. Concentrations of HM in storm water are not determined by rainfall event duration (t r ) as shown by a low correlation coefficient (Table 3). They rather depend on precipitation depth (P tot , P (10) ). This observation complies with the results reported by Rocher et al. in [47]. Based on the data obtained from Paris, the researchers observed that the process of pollutant removal from the atmosphere with rainfall depends mainly on its depth rather than other characteristics. Permeable surfaces (green spaces-GS and gravel roads-GRD) within the catchment, characterised by high water retention capacity, to some extent reduce the content of trace elements in stormwater intercepted by watertight drainage systems. This observation is confirmed by the occurrence of negative correlation relationships between the analysed pollution indicators and the share of permeable surfaces in total area of the catchments.
Stormwater quality also depends on wind velocity, which has substantial impacts on transport and deposition of atmospheric pollutants (Table 3) and may provide a reasonable explanation for the negative correlation between the catchment area and the studied HM. With respect to the analysed catchments, large areas are relatively less urbanised and have a considerable share of green spaces. In addition, the negative correlations (from −0.34 to −0.71) between wind velocity and HM concentrations in stormwater may indicate that wind has a substantial impact on the volume of dry pollutant deposition.

Prediction of HM Concentrations in Stormwater
The results also indicate that it is possible to develop a model for assessing stormwater HM concentrations on the basis of precipitation characteristics, catchment parameters, and weather conditions, including visibility. The model proposed in the study can have practical applications to urban catchments. Additionally, the model makes it possible to examine the impact and interactions between variables analysed, both for rainfall and snowmelt events. These issues were investigated in many research studies; however, they focused on homogeneous catchments [46,48]. Yet, with respect to urban catchments, a large number of problems have not been tackled. The model proposed allows analysis of the impact of many different phenomena (e.g., deposition of airborne pollutants and their washing off the ground surface) and their variability, which depends on land use type. Comparison of simulation results and literature data shows that the ANN model proposed in the study accounts for nonlinear variables, which contributes to its predictive capacity (Table 4, Figure 5).
The parameters that characterise the structure of MLP models, namely number of neurons, activation functions in the hidden and output layers, are shown in Table 5. The values of measures of the model fit (R, MAE, RMSE) to experimental data are also given in the table. Based on the data from Table 5, it can be seen that for individual networks, a number of neurons in the hidden layer ranges from 28 to 35. Additionally, the analyses demonstrate that the number of neurons obtained for individual MLP models is not greater than that recommended [25]. The results show that for Ni and Zn, selected independent variables ensure satisfactory predictive capacities of HMC mathematical models-R 2 > 0.90 ( Figure 5). With regard to Cr and Pb, the fit of simulated values to the measured ones is also satisfactory (0.83 < R 2 < 0.86). In addition, Figure 5 shows that the smallest scatter is found for concentrations lower than 0.08 mg·dm −3 (Cr) and 0.4 mg·dm −3 (Pb). With respect to Cu, good predictive abilities of the model are observed only up to the value of approx. 0.2 mg·dm −3 . The value of R 2 = 0.256 for Cd shows poor fit between ANN predictions and measurements, especially for concentrations higher than 0.03 mg·dm −3 ( Figure 5).
The analyses indicate ( Table 4) that the pattern of variation in HM concentrations for catchments with diversified land use is much more complex than in homogeneous catchments [46]. Additionally, in heterogeneous catchments, mechanisms affecting the quality and quantity of stormwater are far more complex. Consequently, taking into account simulation results, it is difficult to specify variables that decidedly determine HM concentrations. The concentrations of the analysed HM are affected to the greatest extent by the following values: G, P (10) , T me , w me and A (Table 3). It is necessary to conduct further investigations and select adequate computational models to identify key independent variables.

Conclusions
The computations performed for this study demonstrate it is possible to develop an ANN model to simulate variation in HM concentrations in stormwater and meltwater. The model is innovative because it attempts to simulate different mechanisms that affect HM concentrations. The use of a single mathematical model for that purpose has not been reported in the literature so far.
The computations show it is possible to simulate HM (Ni, Cd, Cu, Zn, Cr and Pb) concentrations in stormwater flowing from catchments with different types of land use. The ANN-based model employed to that end ensures satisfactory reliability. The results reported in the study demonstrate that HM concentrations in stormwater are affected by the physical and geographical characteristics of the catchment, and also by parameters describing atmospheric conditions and air quality. It is also shown that for heterogeneous catchments (as those investigated in the study), it is difficult to clearly specify key variables that decide the selected HM concentrations. To validate the ANN models developed for the study, it is necessary to conduct further field investigations.
The methodology proposed in the study is universal in character. The analysis of spatial data making up the catchment characteristics, which are retrieved from available databases (orthoimagery, base maps, digital elevation models-http://www.gis.kielce.eu/), measurements of atmospheric conditions and air quality, (alternatively available monitoring network data-https://en.tutiempo.net/ climate/poland, http://smjp.kielce.pios.gov.pl/) and measurements of HM concentrations in stormwater can be combined. That provides a tool which can be applied to model HM concentrations in different scenarios. The ANN method employed in the study is commonly used in computational packages. As a result, the development of the model and the optimisation of its structure are not complex tasks (the parameters are automatically estimated by the programs, e.g., MATLAB), and such an approach will gain in popularity in engineering practice. The model can be applied to spatial development and planning. It can be used to seek optimal solutions in stormwater management intended to reduce the HM content flowing from the catchment. Funding: This research was funded from the programme of the Minister of Science and Higher Education entitled "Regional Initiative of Excellence" in 2019-2022 project number 025/RID/2018/19 financing amount PLN 12,000,000".

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.