## 2 Materials and Methods
### Statistical Modeling  
The effect of the water source type and year of measurement on the gross alpha and gross beta activity and radon-222 content considering local variations was evaluated. Because the differences between deep (Deep) and subsurface (SSurf) water were not significant, we considered the water source as a factor with two levels, groundwater (Ground) and surface (Surf). The year of measurement was considered as a numeric variable starting with 2017 (the beginning of the study), because we proposed to test the potential linear trend of the radiation level. We used the GLMM with gamma distribution (the data were overdispersed and residuals were not normally distributed) and logarithmic link function, including locality as a random factor, because of repeated sampling. When the effect of the random factor was not significant, we reported the results of the corresponding generalized linear model (GLM). The explained variation in GLMs was expressed as the explained deviation. In GLMMs the partition of the variation explained by the fixed and random effects was carried out using the marginal and conditional lognormal R squared. The correlation between water physical-chemical parameters was evaluated based on the repeated measures (within localities) correlation coefficient calculated with package mrcorr [56] in R version 4.1.0 [57], using bootstrapping with 100 resamples. We ran the correlation analyses both with and without the outliers, and reported all results when significant (or marginally significant). To evaluate the overall response of physical-chemical parameters to water source and sampling year, we used multivariate linear redundancy analysis (RDA), performed using the Canoco 5.12 software [58]. In the RDA, multiple numerical response variables (in our case the physical-chemical parameters and, separately, the radioactivity level) are regressed against one or, more often, several predictors (explanatory variables) that may be of different type (numerical or categorical--in our case, water source and sampling locality, year and month), which are combined in independent (orthogonal) constrained ordination axes (usually the first two are represented in the ordination diagrams) which explain most variation in the response variables, their explanatory power decreasing gradually. In the ordination diagrams the response variables are represented by arrows. The longer the arrow, the better the response variable is explained by the predictors. The projection of the arrowhead on the axes gives the degree of dependence of the response variable to the constrained axis. The angle between arrows indicates the correlation between the response variables (positive for acute angles and negative for obtuse ones). Numerical predictors are also represented by arrows and their projection on the axes illustrates their contribution to those axes. However, in our case the numerical predictor (the year) was not significant in the multivariate models. Levels of factors are represented by the centroids of the data points corresponding to each level. Their effect on the response variables is given by the projection on the arrows. When the differences between deep and subsurface water were significant, we considered the water source as a factor with all three levels (Deep, Ssurf, and Surf). The significance of ordination axes was tested by the Monte Carlo permutation test with 999 unrestricted permutations per test [59]. The significance of physical-chemical parameters' responses (either positive or negative) to individual predictors was evaluated visually, constructing the t-value biplots (with van Dobben circles). To evaluate the