Application of Environmetrics tools for geochemistry, water quality assessment and apportionment of pollution sources in Deepor Beel, Assam, India

BOD₅

5-day Biochemical Oxygen Demand

Ammonia

Cd

Cadmium

Ca²⁺

Calcium

COD

Chemical Oxygen Demand

Cl⁻

Chloride

Cr

Chromium

Cu

Copper

DO

Dissolved oxygen

EC

Electrical conductivity

F⁻

Fluoride

Fe

Iron

Pb

Lead

Mg

Magnesium

Mn

Manganese

Nitrate

OrgN

Organic nitrogen

Phosphate

K⁺

Potassium

Na⁺

Sodium

Sulphate

TA

Total alkalinity

TDS

Total dissolved solids

TH

Total hardness

TKN

Total Kjeldahl nitrogen

TSS

Total suspended solids

WQ

Water quality

Rapid and significant changes in the land-use and land-cover (LULC) patterns across the globe have been observed to have profound effects on the natural water bodies. Numerous studies citing the alteration of the natural state of our water bodies, both quantitative as well as qualitative, have raised serious concerns among researchers regarding global water sustainability (Shi et al. 2017; McCartney et al. 2018). The quantitative revival of the water bodies can be achieved through change in our land-use patterns and judicious utilization of our water resources as well. The qualitative restoration process, however, is a very tedious task and requires a more comprehensive and detailed study on the pollution sources and contaminating factors responsible for the continuous degradation of the water quality. This is owed to the dynamics of the natural systems. Thus, a continuous monitoring and assessment process, both spatial and temporal, for quality check purposes becomes inevitable. Continuous monitoring provides a reliable yet sophisticated water quality (WQ) dataset, which becomes quite ineffective when it comes to interpretation, owing to its complexity (Vega et al. 1998; Simeonov et al. 2003; Iscen et al. 2008; Chow et al. 2016; Hajigholizadeh & Melesse 2017; Singh et al. 2019). However, in recent years an Environmetrics approach (that is, application of various statistical techniques), has eased how the water quality datasets are understood, including the classification of spatially distributed monitoring sites and the primary factors or contaminants responsible for the deterioration of the water quality (Jha et al. 2014; Machiwal & Jha 2015; Bodrud-Doza et al. 2016; Chow et al. 2016; Kumar et al. 2017).

Hierarchical clustering of sampling sites through cluster analysis (CA) and the identification of probable pollution sources through principal component analysis (PCA) have been widely used and accepted. The use of discriminant analysis (DA) as a supervised pattern recognition tool for the recognition of the most significant water quality variables accountable for spatial and temporal variability has also been used more recently compared to the other two methods (Hajigholizadeh & Melesse 2017). However, these statistical tools cannot solely quantify the contribution of potential pollution sources. For this purpose, various receptor models, such as the positive matrix factorization (PMF) model, can be made use of. The PMF models were initially applied to the dataset pertaining to the atmospheric pollution to determine how much is the contribution of various pollution sources. Only in recent times, have they been applied to the water quality (WQ) datasets along with PCA and CA for quantifying the contributions of pollution sources (Zhao et al. 2013; Mustaffa et al. 2014; Chen et al. 2015; Gholizadeh et al. 2016).

Deepor Beel is one of the largest freshwater wetlands of the Brahmaputra Valley, located in the heart of Guwahati, India. It can be considered as the lungs of the bustling city, assimilating a large amount of municipal and industrial wastewater along with agricultural runoff daily (MoEF 2008). The dilution and dispersion of pollutants in the wetland reduces the pollutants’ concentration in the waters, whose ultimate fate lies in the Brahmaputra river. However, in recent times, Deepor Beel has been in a sorrowful state, primarily because of various anthropogenic reasons. Once a 40.14 km² spread wetland is now just a 4.1 km² plot of wetland with human dwellings and industries surrounding it from almost all directions (MoEF 2008). Thus, there is a need for an intensive scientific research with continuous monitoring of the wetland water quality to ensure that the conditions do not deteriorate any further and at the same time help the authorities in the decision-making process to implement better restoration programs.

Given the above considerations, the vast WQ dataset matrix generated as a result of continuous monitoring from 23 sampling locations for 28 parameters from October 2017 to February 2019 (about 24,000 observations) were subjected to (i) hierarchical clustering through CA for identification of the sampling locations having similarities and classifying them accordingly; (ii) confirmation of the grouping of sites from CA through discriminating functions obtained from DA; (iii) identification of presumable contributing pollution sources based on the results obtained through PCA and (iv) quantification of the contribution from various sources through source apportionment from a PMF model. The results obtained from this study are expected to be of substantial help for efficient water management in Deepor Beel.

Deepor Beel (Figure 1) is situated in the heart of the Indian city of Guwahati, Assam. It holds an essential place in the international space owing to its inclusion as a Ramsar site under the Ramsar Convention of Wetlands (Ref. No. 1207) (Bhattacharyya & Kapil 2010). Two primary outlets act as sources of water for Deepor Beel, one being the Basistha River, joining the wetland from the south-eastern part, and the other being Khanajan, connecting Deepor Beel to the Brahmaputra River in the north. The entire watershed of the wetland has been delineated using ArcGIS (Version 10.2) and the land-use and land-cover (LULC) map prepared using ENVI (Version 4.7) (Figure 1(a)). It is clear from the LULC map that a significant portion of the watershed is either vegetation (Western part) or forest cover (the Rani and Garbhanga Forest Reserve in the South); however, the wetland is also surrounded by a dense urban built-up area, which has been a primary reason for its continuous deterioration, both in expanse as well as quality of water. The presence of a landfill site (Boragaon landfill in the East) and the growth of small and large scale industries (in the West) in recent times have further added to the degradation of water quality. Given these problems, there is a need for an immediate water quality monitoring program to identify the plausible pollution sources and determining which parameters are primarily responsible for the deterioration of the water quality and also to quantify their contribution.

Hence, based on the reconnaissance survey of the entire wetland, 23 sampling sites were chosen such that they are entirely accessible and are also proximate to various pollution sources (Figure 1(b)). Sampling was conducted continuously from October 2017 to February 2019, twice every month. The water samples were collected from a depth of 30 cm below the top water level in decontaminated and thoroughly scrubbed and washed (in laboratory-grade detergents) reagent bottles. Two bottles were utilized for carrying samples from each location with proper labelling; one for all the physiochemical analyses, and the other for metal analyses. The containers used for metal analyses of water samples were acidified with HNO₃ after the collection to avoid precipitation. To minimize the changes in the samples, chemically or biologically, they were preserved in an icebox during the transfer from the sampling sites to the laboratory and later in the refrigerator when analyses were not carried out. The entire processes of sample collection, preservation, transfer, and investigations complied with the standard protocol (APHA 2012). The WQ parameters investigated in this study, together with their units of measurements and the methodologies adopted for their determination, are listed in Table 1.

Hierarchical agglomerative clustering technique is a simple algorithm for partitioning the datasets into groups (Singh et al. 2004; Shrestha & Kazama 2007). Each dataset is first assumed to be a unique cluster. Each of these clusters is then joined step by step depending on their proximity to each other, till a stage is arrived at when all the datasets are placed under a single large cluster. The fundamental algorithm for hierarchical agglomerative clustering depends on large inter-cluster and low intra-cluster variance. The basic steps involved in the clustering process have been explained below:

The matrix obtained is called the distance matrix, which provides the unit of distance between each dataset (Equation (1)). It is always a lower triangular matrix (m × m), with its diagonal elements being zero, which is quite apparent.

(1)

From the distance matrix, the clusters closest (i.e. the minimum value in the matrix) to each other are categorized under a single cluster (in this case, let's suppose d_(m−1)3) (Equation (2)). A distance matrix is again created for (m − 1) observations {(m − 1) × (m − 1)}.

(2)

The clusters closest to each other are further categorized to a single group, and the process is repeated ‘m’ times till a single value is obtained; that is, the entire dataset is grouped into a single cluster.

The results of hierarchical clustering is usually represented by a dendrogram which provides essential visual insight into the clustering process.

DA was carried out on the raw dataset for both standard, as well as stepwise modes. The standard method involves a step-by-step build-up of the model of discrimination in which all the variables are assessed and evaluated to determine which one will contribute most to the discrimination between groups individually, at each step. That variable is then incorporated in the model, and the procedure starts again. On the other hand, in stepwise mode, all variables are included in the model initially, after which the variable contributing the minimum to the estimation of group membership is eliminated at each step. Thus, only the ‘important’ variables in the model which contribute the most to the discrimination between groups are kept as the result of a successful discriminant function analysis (Astel et al. 2009).

The principal component analysis (PCA) forms a powerful tool for minimizing the problem of data overfitting to a model, which has led to several confusions due to the inclusion of too many attributes and variables into it. The primary objective of the PCA is to convert a dataset having high dimensionality to a lower dimensionality dataset, without actually losing vital information, thus effectively reducing the problem of overfitting (Vega et al. 1998; Wunderlin et al. 2001; Chen et al. 2007). The reduced set of variables formed are termed as principal components (PCs), which depend on the number of attributes the dataset possesses, as well as its dimension. The number of PCs formed for a particular dataset is always less than or at a maximum scale, equal to the number of attributes given for model preparation. Also, the PCs formulated should be independent and orthogonal to each other. The necessary steps involved in the formulation of PCs is described below.

Suppose we have ‘m’ attributes (A₁, A₂, …, A_m), the values of those for ‘n’ sampling locations (variables) (S_{11, …,} S_nm) can be represented as shown in Equation (7)., the first operation performed is the computation of the mean values for all the attributes (⁠⁠).

(7)

The next step includes the preparation of the covariance matrix for all the attributes (Equation (8)) based on the formula given in Equation (9). The matrix formulated will be an (m × m) matrix (C).

(8)

(11)

PCs depend on the eigenvalues of the corresponding attributes. The eigenvalues are first considered in the descending order of their numerical values. The eigenvectors of the corresponding largest eigenvalue become PC1 and subsequently, other PCs are also determined (the corresponding eigenvectors) based on their descending order of eigenvalues.

Varimax rotation, Keiser-Meyer-Olkin (KMO) (Sahu et al. 2013) and Bartlett's sphericity (Dalal et al. 2010) were employed for the orthogonal rotation, evaluating the adequacy of sampling and estimating the applicability of PCA to the crude dataset, respectively. Eigenvalues greater than 1 were only considered for accountability of the factor loadings.

Considering a matrix A_a×b; m and n denoting the dependent and independent variables respectively (Equation (15)),

(15)

The following inferences can be established:

• A_a×b covers a multi-dimensional space (A ɛ ℜ^a×b).

• A_a×b may be sparse.

• This may lead to noisy features in additional data processing, ultimately providing poor classification.

A multiplicative update algorithm is usually employed for solving this minimization problem, which is carried out using the following steps:

Iterations are carried out using step 2 till W and H converge with each other and become stable.

All the statistical techniques were carried out using various software; Hierarchical clustering, DA (both standard and stepwise) and PCA using IBM – SPSS Statistics (Version 25) and the PMF model was formulated using EPA – PMF (Version 5.0).

A summary of the overall descriptive statistics of various WQ parameters is shown in Table 1. It was observed that most of the elements exceeded the acceptable limit as set by IS:10500 (2012), which is an indicator for water contamination of the wetland. Only Cd exhibited negative skewness, while all other parameters exhibited positive skewness, with a majority lying in the highly positively skewed (skewness > 1) category. This indicates that they have a long right-hand tail curve, when normally distributed. As far as the kurtosis of the observed dataset is concerned, about 57.2% of the parameters were observed to be of the Plytokurtic category (Kurtosis < 3), thus indicating a flatter, normally distributed curve. On the other hand, the rest were found to be Leptokurtic (kurtosis > 3), indicating a higher peak representation of the normally distributed data.

The hierarchical clustering of the 23 sampling locations of Deepor Beel resulted in classifying them into three major groups (clusters 1, 2, and 3), represented through a dendrogram (Figure 2(a)). The clustering was then validated using the agglomeration schedule coefficients, which generated a plot like a scree plot, only backward (Figure 2(b)). The number of agglomerative schedule coefficients represents the amount of heterogeneity in a cluster solution (known as the Stopping Rule in hierarchical clustering). Hence, our main objective was to maximize the heterogeneity (distinction) for obtaining better clustering results. A single factor solution (single cluster) is the most distinct, with the next best being the split-up of the cluster into two parts; that is, the formation of two clusters. However, often, two numbers of clusters are not considered to be ideal. Hence, the check is carried out for three clusters, then four, and so on until no significant jump (flatter slope) arises, wherein the distinction is considered unsuitable. In the current study, it was observed that a steep slope was produced until the third cluster process, after which the slope of the curve became relatively flatter, thus indicating a lack of difference in the homogeneity among the clusters. This can also be visualized in Figure 2(a), where distinct clusters were observed until the third clustering, after which the distance between the clustering schedule decreased significantly with the formulation of new clusters. This provided a relatively unclear picture about the further clustering process. Thus, three clusters (Figure 2(a)) were considered for the analysis. The sites categorized into three clusters were then plotted on a GIS platform for better visualization of the sampling locations (Figure 2(c)). It was observed that the three locations belonging to cluster 3 were in the closest proximity to Boragaon landfill site, thus providing evidence for maximum possible contamination among all the locations. Furthermore, cluster 2 categorizes all the sites located in the middle portion of the wetland, which has remained devoid of any significant anthropogenic contamination, thus indicating the sites having the least pollution. Lastly, the locations categorized under cluster 1 were observed to be primarily belonging to the areas proximate to the industrial complex (Western side of the wetland), or between the landfill site and the middle area. This is evidence of a moderate level of pollution, the sections between the landfill and the central portion of the wetland acting as a transition zone between high and medium contamination, respectively.

The discriminant analysis (DA) technique was employed on the raw WQ dataset to determine the spatio-temporal variation of the WQ parameters. Standard and stepwise modes were used for constructing the classification functions (CFs) (Table 2). It was observed that only 9 (EC, TDS, TSS, ⁠, Na⁺, Mg, Cd, Pb and OrgN) out of 28 parameters (based on Fisher's Linear Discriminating Functions) were responsible for the discrimination, as well as the spatial variability among the three clusters. The generation of the CFs was validated using the classification matrix (Table 3), which provided correct classifications among the discriminating functions using the stepwise mode of DA for 73.9% of the cross-validated grouped cases, thereby providing an excellent result for the spatial variations. The scores of the two functions were plotted (Figure 3(a)) along with the Wilk's lambda (λ) values of the discriminating functions (Figure 3(b)). The values of λ varied from 0 to 0.14, which is numerically insignificant, thus indicating that the classes or groups formed are distinctive and there is minimum overlapping between them.

PCA employed on the normalized WQ dataset resulted in the formation of two primary principal components (PC 1 and PC 2), with eigenvalues >1. KMO and Bartlett's sphericity test results (Table 4) revealed the measure of sampling adequacy to be 0.742, with 253 degrees of freedom and a chi-square value of approximately 3,743.39 and a value of p close to zero (p < 0.05). This indicates the validity of the PCA due to their statistical significance between the variables (Kaiser 1974). Figure 4 represents the rotation factor matrix (Varimax rotation) of the two PCs. Values of PCs less than ±0.3 were considered to have weak loadings, while those higher than ±0.7 were regarded as PCs having a significant contribution with strong loading values. The PC values lying between ± (0.3–0.7) were considered to have a moderate impact. It was clearly observed that COD, BOD₅, TDS, EC, Cd, Cr, Fe, Pb, Cu, Mn, and Mg exhibited strong positive loadings for PC 1, indicating substantial anthropogenic interference such as the Boragaon landfill and the industries in the wetland. Strong negative loading of DO displayed good results in PC 1 as the values of DO vary inversely with the BOD/COD values. Furthermore, moderate loading values for Cl, ⁠, turbidity, TKN, and OrgN indicate a secondary source of organic contamination to the wetland (possibly from the Basistha River, which drains its water, primarily composing municipal wastewater, into Deepor Beel). PC 2, however, showed a strong positive loading for TSS while having moderate loadings for TH and Ca, indicative of a natural surface water runoff being the major probable contributor.

Figure 5 presents an overall summary of the PCA carried out for the WQ dataset. It provides the parameters (variables) that are associated more with a particular component along with the correlations existing between them. For example, it can be observed that BOD₅, EC, TDS, COD, and all the metals were observed to have a higher affinity towards PC 1. Similar is the case for TA, DO and pH. However, they lie on the opposite quadrant, thus displaying a highly negative correlation. Similarly, TSS, TH, Na⁺, Ca²⁺, ⁠, and F⁻ are affiliated more towards PC 2, with F⁻ and being negatively correlated to other variables. Figure 6 is a spatial representation of the PC loading values (both PC 1 and PC 2). Clearly, PC 1 had higher loading values near the landfill site and the industrial zone, representing sites having higher anthropogenic contamination factor (Figure 6(a)). Similarly, PC 2 showed higher values towards the central region of the wetland (Figure 6(b)), thus indicating natural factors of contamination of Deepor Beel.

The source apportionment and various factor contributions were estimated using the EPA – PMF software, which takes into account the concentration and uncertainty factors as input to the model. The resulting input parameters then generated a signal/noise ratio for each of the WQ parameters. Based on the S/N ratio, the parameters were classified as Strong (S/N > 2), Weak (0.2 < S/N < 2) and Bad (S/N < 0.2) respectively. Four factors; 1, 2, 3, and 4 were considered for the PMF analysis to assess each of their contributions to the water pollution of Deepor Beel, based on the iterations for each factor obtaining the minimum Q(Robust)/Qexp value (Table 5). The model was then simulated, considering all the input parameters and the four factors, providing results, as shown in Figures 7 and 8. It can be seen that Factor 1 is primarily dominated by BOD₅ (69.3%), COD (71%), Cr (78.5%), Fe (71.9%), Mn (70.6%) and Pb (58.6%). This indicates a source rich in organic matter, as well as heavy metals. Thus, leaching from the Boragaon landfill site can be considered as the major source of contamination for Factor 1, since the municipal landfill site leachates are highly organic in nature owing to the various organic wastes (e.g. food wastes) being dumped into the landfill. Also, the inorganic wastes rich in heavy metals such as batteries, scrap and unused metals contribute to the leaching of heavy metals (Chu et al. 1994). Furthermore, DO (79.1%), pH (46.9%), TA (53.6%), TH (42.3%), Turbidity (51.4%), K (48.7%) and Ca (51.2%) were observed to have significant contribution to water contamination of Deepor Beel, pertaining to Factor 2. This is an indicator for a natural cause of pollution such as the surface water runoff. Similarly, for Factor 3, TDS (57.5%), (55.4%), (63.5%), Na⁺ (63%), Mg (83.2%), Cd (73.7%), and Cu (74%) were found to be the predominating parameters. Industrial wastewater effluents enriched with contagious trace metal elements along with cations and anions might be the primary factor responsible for this. Finally, TKN (67.6%), (68.8%) and OrgN (67.2%) were found to be the parameters affecting Factor 4 the most. This indicates contamination from a source primarily rich in nutrients. Thus, discharge from Basistha River (carrying large concentrations of domestic sewage) might be the primary factor responsible. Figure 8(a) and 8(b) represent the Factor Figureprints and the normalized contributions for all the 28 WQ parameters, providing a summary of the PMF results.

This study presented a detailed insight into the geochemistry and assessment of water quality as well as apportioning of various pollution sources contributing to the water contamination of Deepor Beel using different Environmetrics tools. Critical concluding remarks of the study are as follows:

• Hierarchical clustering categorized the 23 sampling locations of the study area into three statistically distinct clusters. The clustering process was further validated using the ‘Stopping Rule’. It was observed that the sites closest to the Boragaon landfill, locations present in the central portion and the industrial zone of the wetland formed separate and distinct clusters based on the site similarities. They were considered as sites of high, low, and moderate contaminations, respectively.

• The discriminant analysis provided the significance of clustering by discriminating the parameters responsible. Nine (EC, TDS, TSS, ⁠, Na⁺, Mg, Cd, Pb and OrgN) out of 28 parameters were found to be responsible for the discrimination, as well as the spatial variability among the three clusters. Low Wilk's λ values also validated that the groups formed were distinctive with minimum overlapping between them.

• The principal component analysis applied to the normalized dataset yielded two principal components, which aided in the identification of potential pollution sources. PC 1 indicated anthropogenic contamination as the primary source, while PC 2 suggested natural phenomena such as surface water runoff to be the primary source. Spatial representation of the sites also revealed similar results.

• The WQ dataset matrix was finally subjected to factorization, generating a model for source apportionment. Leaching from Boragaon landfill site, surface water runoff, discharge of effluents from the industries in the wetland and discharge from Basistha River were found to be the significant factors for the pollution of Deepor Beel.

• The results obtained from this study would be of immense help for restoration and revival of Deepor Beel.

The authors declare that they have no conflict of interest.

The authors extend their sincere thanks to the anonymous reviewers for their valuable suggestions, which helped improve the quality of the manuscript substantially.

All relevant data are included in the paper or its Supplementary Information.