## 2 Materials and Methods
### Statistical Analysis and Mapping  
The study of the spatial distribution of trace elements in lagoon waters was carried out using geostatistical analysis. The geostatistical method consists of describing spatial patterns (semi-variograms) and attribute value prediction (kriging). It uses the available sampling point values to estimate those of the unsampled locations to produce maps showing the trends of spatial distribution of the studied parameters in the area. Indeed, distribution prediction maps were performed using the Ordinary Kriging Interpolation (OKI) technique of the ArcGIS 10.2.2 software. The ordinary kriging is estimated by a linear combination of the observed values with weights using the following equation [41; 42; 43; 44; 45].  
\[Z^{*}(x_{0})=\sum_{i=1}^{n}\lambda_{i}Z(x_{i})\]  
where Z* (x\({}_{0}\)) is the estimated value of Z at point x\({}_{0}\), Z(x\({}_{i}\)) is the sampled value at the point xi, and \(\lambda_{i}\) is the weight placed on Z(x\({}_{i}\)).  
The particularity of the OKI technique is to produce more accurate data prediction and provide checking elements of the accuracy of predictions [29; 43]. For that, it is the most used method in geostatistical analysis. The most accurate spatial autocorrelation (variogram) model of each trace element is selected using cross-validation. The mean errors (ME) and root mean square error (RMSE) values were used to check the accuracy of the results of the models and test their robustness. The Nugget (\(C_{0}\)) and Sill (\(C_{0}+C\)) ratio was calculated in order to determine the spatial dependence (SD) between sampling points or the degree of spatial variability between them. The SD is calculated using the following equation [29; 43; 42]:  
\[SD\ (\%)=\frac{C_{0}}{C_{0}+C}\]where the nugget (\(C_{0}\)) is the variability at zero distance; sill (\(C_{0}+C\)) is the maximum variance between pairs of data. The spatial dependence (SD) expressed in percentage (%) is interpreted according to [PERSON] et al. (2018) as follows: 0-25% (high), 25-75% (medium), and 75-100% (low) spatial dependence between sampling points. The maximum separation distance between the pairs of data is assessed using the range values ([PERSON] et al., 2019; [PERSON] et al., 2020; [PERSON] et al., 2020).  
The interrelationships and similarities between the different variables and sampling sites studied were determined by performing multivariate analyses such as principal component analysis (PCA), Pearson correlation, and Cluster analysis (CA) using [PERSON]'s method ([PERSON], 1974; [PERSON], 1975; [PERSON], 1976). Before carrying out the principal component analysis (PCA), the adequacy of the data was assessed using the Kaiser-Meyer-Olkin (KMO) test and the sphericity test of Bartlett. The PCA was used to assess the typologies of water quality and highlight the interrelationships between physicochemical parameters and trace elements in order to deduce their possible common sources and the determinism of their distributions in the different samples. The grouping of sampling sites using PCA was carried out manually by observing the trends of raw data and the projection of sampling sites in the factorial plan. The CA method allows one to group observations or variables according to the similarities that exist between them or not. Pearson correlation analysis allowed us to highlight the links between variables. Since these multivariate analyses are highly sensitive to extreme values, all the data were normalized by calculating z-score values before carrying out the analyses according to the following equation ([PERSON] et al., 2019; [PERSON] et al., 2020):  
\[K_{ij}=1+\frac{\left(X_{ij}-X\right)}{S_{i}}\]  
where \(K_{ij}\) is the normal value of \(K_{ij}\) for the \(i^{th}\) variable of the \(j^{th}\) individual, X is the mean of the \(i^{th}\) variable, and \(S_{i}\) the standard deviation of the \(i^{th}\) variable. The multivariate analyses were performed using the STATISTICA 6.1 software.