This section presents the results of various analyses, including delineating CAs, classifying LULC, analyzing LULC changes within CAs, conducting statistical analysis of WPI measurements, performing geo-visual analysis of WPI values for dominant LULC by CA, and analyzing the correlation between LULC composition by CA and measured WPI values using CCA.
WPI- Descriptive statistical and time series analysis
The WPI values for each measurement point during the five-year study are available in Annex A1. Annex A2 includes the table with the descriptive statistics of the measured WPIs. According to NOM 001-SEMARNAT 2021, which establishes permissible limits of contaminants in wastewater discharges in receiving bodies owned by the Mexican state (provided in Table 1), the TSS and FC values exceeded the allowed limits. The WTemp was within the permitted limit (35°C), but in the lower part of the basin, in P6 and P7, the values closest to this limit (33.1°C) were measured in 2019. pH values (8.7 to 9.3 UpH) around the limit (9 UpH) were measured in 2017 in the upper and middle parts of the basin (P1, P2, P3, and P4). The highest concentrations of TN (36.24 UpH, surpassing the limit of 30 mg/L), TP (5.94 mg/L), and BOD (120.94 mg/L) were measured in 2017 in the lowest part of the basin (P7). The TP and BOD measurements did not surpass the corresponding limits (18 mg/L and 180 mg/L, respectively).
The time series of the measured WPIs corresponding to the study period are shown in Fig. 3. The lowest Cond values (up to 388 µS/cm) for the entire study period were measured in the upper part of the basin (P1 and P2). The highest values (> 1500 µS/cm) were observed in March 2019 in P3, while at the other measurement points, the values fluctuated around > 1500 µS/cm for all years. However, the values of this parameter were below the permitted limit throughout the period.
In March 2017, the pH values were above the norm in the middle part of the basin (P3 and P4, 9.2 and 9.3 UpH, respectively), while values of 8.7 were present in the upper part of the basin (P1 and P2).
The TSS exceeded the established limit throughout the study period and for all stations. The values recorded in the upper part of the basin (P1 and P2) were closest to this limit during all measured days. The highest values (close to 100 mg/L) were observed for March 2016 and March 2019 in P3 and for April 2017 in P7, corresponding to the lowest part of the basin. WTemp was within the allowed limit (35°C) for all points and dates. However, the highest temperature (33.1°C) of this WPI was observed in March 2019 in the lower part of the basin (P6 and P7). The WTemp values increased from the highest to the lowest part of the basin, as analyzed via geo-visual analysis in the following section. In March 2017, a WTemp lower than 25°C occurred in the lower part of the basin (P6). The other values of the other seasons and dates fluctuated around 27.5°C.
The levels of BOD and TP found in the basin were within the allowed limits for all measurements. The highest BOD value of 120.94 mg/L was observed in April 2017 at point P7 in the lower part of the basin. The rest of the basin had values close to 0 mg/L. However, the concentration of TN found in April 2017 at point P7 in the lower part of the basin exceeded the limit, with the highest value of 36.24 mg/L. The TN concentrations recorded from the other points and dates were close to 0 mg/L.
The FC values exceeded the allowed limit for practically all dates and parts of the basin, with the highest values observed in the lower part of the basin at point P7 and in the middle part of the basin at point P3 in March 2018. The lowest values of this parameter were observed in the upper part of the basin at points P1 and P2.
CA delineation
Figure 2 displays the CA for each measurement point, including the area for each CA. We identified seven CAs. The CAs in the upper part of the basin have larger areas than those in the middle and lower basin parts.
Validation of the LULC classification
The results of the LULC classification validation are presented in Fig. 4. The overall agreement, relative to the area (Pontius and Santacruz, 2014), for all classified images was 92% or greater for all the studied years. As shown in the figure, natural land cover had the highest overall agreement and percentage of the domain for the studied years, followed by agriculture, built areas, and water. Water had the lowest overall agreement and percentage of the domain, with the highest omission and commission values for all the studied years.
LULC change analysis
Figure 5 and Table 3 show the LULC classification results for each measurement point and corresponding CA over five years and for the whole basin. The basin covers an area of 45349.2 ha. On average, for the five years studied, 84.3% were natural cover, 13.7% were agricultural land, 1.9% were built-up areas, and 0.1% were water. As shown in the table, during the study period, the natural cover decreased by 2.3% (888 ha), while agricultural land increased by 7% (439 ha), water cover increased by 12% (14 ha), and built-up areas showed the greatest increase of 67% (435 ha). This means that agricultural and built-up areas expanded into natural areas. Although the natural cover is predominant in the basin and the decrease in this type was low, it still represents an annual loss of 0.46% compared to the initial year (2016). The observed increase in water areas from 6.6 ha to 20.5 ha showed significant fluctuations in the intermediate years, with values higher than 80 ha, indicating strong seasonal influences on this type of LULC.
The following are the resulting percentage changes in natural coverage for each CA from 2016 to 2020: CA1 decreased by 2.3%, CA2 decreased by 1.0%, CA3 decreased by 2.6%, CA4 decreased by 17.2%, CA5 decreased by 59.5%, CA6 decreased by 98.1%, and CA7 decreased by 54.7%. These values indicate decreases of 536.9 ha, 84.6 ha, 167 ha, 6.1 ha, 74.8 ha, 6.1 ha, and 14.8 ha, respectively. The CAs with the highest percentage of natural coverage, such as CA1, CA2, and CA3, had larger areas and the lowest reduction percentages of this cover type. Conversely, the CAs with the greatest presence of built-up areas, including CA4, CA5, CA6, and CA7, had smaller area decreases but also the greatest proportions of natural coverage reduction.
Agricultural coverage underwent some changes from 2016 to 2020. CA1 increased by 24%, CA2 increased by 6.9%, CA3 decreased by 3.4%, CA4 decreased by 15%, CA5 decreased by 3.7%, CA6 decreased by 34.5%, and CA7 decreased by 2.6% (+ 505 ha, + 65 ha, 90 ha, -17 ha, -17.6 ha, -2.1 ha, and − 4.5 ha, respectively). The areas with the greatest amount of natural coverage (CA1 and CA2) experienced an increase in agricultural coverage, while the others experienced a decrease.
During the entire period, the following changes were observed in the built areas: CA1 increased by 204.6%, CA2 increased by 270.1%, CA3 increased by 72.2%, CA4 increased by 113.4%, CA5 increased by 51.9%, CA6 increased by 12.8%, and CA7 increased by 57.1%. This means that the surface areas of CA1, CA2, CA3, CA4, CA5, CA6, and CA7 increased by 28.9 ha, 18.8 ha, 254.0 ha, 20.09 ha, 88.6 ha, 7.1 ha, and 17.8 ha, respectively. CA1 and CA2 showed the greatest percentage increases compared to their initial state. However, CA3 had the greatest area increase, followed by CA5.
The increases in the water cover during the study period were as follows: CA1 135.2%, CA2 292.3%, CA3 129.9%, CA4 103.9%, CA5 502.7%, CA6 1471.4%, and CA7 1677.8%. The increases in ha for CA1, CA2, CA3, CA4, CA5, CA6, and CA7 were 2.9 ha, 0.8 ha, 2.9 ha, 1.1 ha, 3.8 ha, 1.0 ha, and 1.5 ha, respectively. The greatest percentage increase was observed in CA1 and CA2, followed by CA5, which had the greatest natural coverage. All CAs increased in the range of 0.8 to 3.8 ha.
In general, during the study period, the predominant land cover in the basin was natural. However, this cover type consistently decreased by almost one percent annually. At the same time, built-up and agricultural land cover steadily increased at annual rates of 10% and 33%, respectively. These trends indicate a growing human impact on the basin, which has led to the modifications mentioned above in terms of land use and land cover (LULC).
Table 3
The LULC variation from 2016 to 2020 for the whole basin and the CA (ha).
LULC types
|
Area
|
2016
|
2017
|
2018
|
2019
|
2020
|
Natural
|
CA1
|
22992.6
|
23806.1
|
24027.2
|
22992.6
|
22928.9
|
CA2
|
8231.0
|
8469.6
|
8557.0
|
8230.9
|
8170.9
|
CA3
|
6103.7
|
6343.2
|
6600.7
|
6103.7
|
6179.8
|
CA4
|
15.6
|
15.1
|
30.3
|
15.6
|
19.9
|
CA5
|
77.2
|
67.8
|
113.5
|
77.2
|
50.8
|
CA6
|
0.1
|
4.1
|
1.3
|
0.07
|
0.12
|
CA7
|
26.7
|
15.7
|
22.3
|
26.7
|
12.2
|
Whole basin
|
37446.9
|
38721.5
|
39352.2
|
3744693
|
37361.9
|
Agriculture
|
CA1
|
2542.23
|
1698.4
|
1519.4
|
2542.2
|
2606.4
|
CA2
|
953.77
|
711.7
|
629.6
|
953.8
|
1005.2
|
CA3
|
2693.53
|
2587.4
|
2218.4
|
2693.5
|
2545.1
|
CA4
|
99.26
|
115.9
|
92.3
|
99.3
|
96.1
|
CA5
|
428.74
|
514.7
|
426.5
|
428.7
|
458.3
|
CA6
|
1.8
|
6.9
|
3.3
|
1.8
|
4.0
|
CA7
|
144.46
|
182.4
|
161.8
|
144.5
|
167.6
|
Whole basin
|
6863.79
|
5817.4
|
5051.3
|
6863.8
|
6882.8
|
Built
|
CA1
|
14.2
|
18.4
|
26.4
|
29.7
|
43.1
|
CA2
|
7.0
|
7.6
|
11.7
|
16.9
|
25.8
|
CA3
|
351.9
|
398.9
|
479.5
|
529.3
|
605.9
|
CA4
|
17.7
|
23.4
|
28.8
|
36.2
|
37.8
|
CA5
|
170.8
|
189.2
|
220.6
|
244.3
|
259.4
|
CA6
|
56.0
|
57.3
|
60.3
|
62.9
|
63.1
|
CA7
|
31.1
|
32.1
|
39.3
|
46.2
|
48.9
|
Whole basin
|
648.6
|
726.9
|
866.7
|
965.5
|
1084.0
|
Water
|
CA1
|
2.2
|
60.6
|
10.5
|
18.9
|
5.1
|
CA2
|
0.3
|
14.0
|
4.6
|
1.4
|
1.0
|
CA3
|
2.2
|
5.6
|
36.5
|
8.6
|
5.1
|
CA4
|
1.0
|
1.5
|
4.5
|
4.8
|
2.1
|
CA5
|
0.8
|
1.4
|
12.5
|
22.8
|
4.5
|
CA6
|
0.1
|
0.05
|
3.4
|
3.5
|
1.1
|
CA7
|
0.1
|
0.05
|
6.9
|
12.9
|
1.6
|
Whole basin
|
6.6
|
83.3
|
79.0
|
72.9
|
20.5
|
Geo-visual analysis of WPI measurements and dominant LULC by CA
Figure 6 includes maps showing the average WPI values during the study period in the context of the predominant LULC type according to CA, proximity to urban areas, and location in the basin. According to the figure, the highest levels of Cond were observed at the exits of built areas. These levels did not show accumulation patterns at the exit of the basin. The highest pH values were observed at the outlets of the CA, which were predominantly covered with vegetation, whether natural or agricultural. The highest TSS values were observed at the exit of the CAs, which had the majority of agricultural cover. Additionally, the concentration of TSS also increased at the exit of predominantly urban CAs. WTemp showed a pattern of increasing values downstream, with the lowest values at the outlet of the highest and least anthropized CA and the highest at the exit of the most urbanized CA close to the coastal zone. This may be related to the altitudinal gradient. WPIs such as BOD, TP, TN, and FC had the highest values at the outlet of the entire basin, which indicates that they followed the logic of accumulation through the drainage network toward the lowest part of the basin. Similarly, FC appeared to be associated with urban areas. Overall, the urban areas presented the highest average WPIs.
Correlation analysis of the WPIs and LULC composition by CA
The results of the bivariate Pearson’s correlation analysis, used to define associations between the measured WPI values at each measuring point and the LULC composition by the corresponding CA, are presented in Fig. 7.
Cond and TSS were strongly correlated (r = 1). The same was observed between TP and TN (r = 0.99), TN and DBO (r = 0.96), and TP and DBO (r = 0.97). Water and agriculture were correlated (0.89). Cond and natural cover were correlated (r = 0.83), as were TSS and vegetation (r = 0.83). Agriculture and built covers were correlated (r = 0.79).
CCA Results
The results of the permutation tests of the significance of the applied models, the LC terms, and the CCA axes are shown in Table 4. The statistical significance of the CCA model for CA1 was weak. The model that included the first two CAs was significant (p = 0.0022); however, natural cover and vegetation, as well as the second canonical axis, were not significant. Therefore, the CCA model was considered not significant. However, the CCA models of CAs 1 to 3 were highly significant (p = 0.001), and the axes were significant, except for built-up areas (p = 0.107). Therefore, we analyzed the results of the corresponding biplot in the following section. The CCA models that included the CA 1–4 and CA 1–5 data were significant (p = 0.002 and p = 0.001, respectively). However, in both cases, only the first axis of the CCA was significant, and the terms agriculture and built areas were not significant, so these models cannot be considered statistically valid.
The CCA including CA 1 to 6 data showed statistical significance (p = 0.001) along with its axes and terms, except for agriculture (p = 0.099). This model, which was considered statistically valid, is further analyzed in the following section. Another CCA model, quite similar to the previous one, was the CCA of CAs 1 to 7 (which corresponds to the entire basin). This model, which was also statistically significant (p = 0.001), is discussed in the following section. The term agriculture (p = 0.156) was the only exception. The last CCA model for CAs 4 to 7 was statistically significant (p = 0.001), but only two of its four terms were significant, and its axes were not significant. Therefore, this model is considered to be statistically nonsignificant.
Table 4
Significance (p) of permutation tests for CCA under reduced models
CA
|
Model
|
LC terms
|
CCA-axes
|
P1/CA1
|
0.3218
|
0.3218
|
0.3218
|
P1/CA1 + P2/CA2
|
0.022**
|
Natural: 0.066 .
Agriculture: 0.127
Water: 0.019 *
Built: 0.044 *
|
CCA1: 0.007 **
CCA2: 0.670
|
P1/CA1 + P2/CA2 + P3/CA3
|
0.001 ***
|
Natural: 0.001 ***
Agriculture: 0.002 **
Water: 0.006 **
Built: 0.107
|
CCA1: 0.002**
CCA2: 0.002 **
|
P1/CA1 + P2/CA2 + P3/CA3 + P4/CA4
|
0.002 **
|
Natural: 0.001 ***
Agriculture: 0.331
Water: 0.053 .
Built: 0.638
|
CCA1: 0.004 **
CCA2: 0.360
|
P1/CA1 + P2/CA2 + P3/CA3 + P4/CA4
+P5/CA5
|
0.001 ***
|
Natural: 0.001 ***
Agriculture: 0.228
Water: 0.011 *
Built: 0.352
|
CCA1: 0.002 **
CCA2: 0.449
|
P1/CA1 + P2/CA2 + P3/CA3 + P4/CA4
+P5/CA5 + P6/CA6
|
0.001 ***
|
Natural: 0.001 ***
Agriculture: 0.099 .
Water: 0.048 *
Built: 0.050 *
|
CCA1: 0.001 ***
CCA2: 0. 043 *
|
P1/CA1 + P2/CA2 + P3/CA3 + P4/CA4
+P5/CA5 + P6/CA6 + P7/CA7
|
0.001 ***
|
Natural: 0.001 ***
Agriculture: 0.156
Water: 0.002 **
Built: 0.041 *
|
CCA1: 0.001 ***
CCA2: 0.003 **
|
P4/CA4 + P5/CA5 + P6/CA6 + P7/CA7
|
0.001 ***
|
Natural: 0.001 ***
Agriculture: 0.289
Water: 0.001 ***
Built: 0.110
|
CCA1: 0.001 ***
CCA2: 0.101
|
Key: <0.001: ‘***’, 0.001: ‘**’, 0.01: ‘*’, 0.05: ‘.’, > 0.1: ‘ ’ |
Significant CCA results
The CCA ordination plots of the eight CA groups, as indicated in the methodology section, are shown in Fig. 8. According to Ter Braak (1986), in the diagram, the sites and WPI measurements are represented by their location in the ordination diagrams, the LULC composition by CA is represented by arrows, and the WPI measurement locations and the arrows of the LULC composition jointly reflect the WPI measurement distributions along each of the LULC composition variables, which are considered the environmental variables. In the following paragraphs, only the significant results of the CCA models are analyzed.
CCA for P1/CA1 + P2/CA2 + P3/CA3
The resulting ordination diagram is shown in Fig. 8c. The first CCA axis explained 52.7% of the variance, while the second axis explained 41.9%, for a cumulative proportion of 94.7%. The diagram indicates that Cond and TSS were mutually correlated, as they were collocated. Moreover, the highest WPI measurements were not correlated with P1 or P2. These measurement points correspond to CA1 and CA2, which are the parts of the basin with the highest natural cover area. In contrast, the WPI measurements appeared close to the measurements at P3, which is the outlet of CA3 and had the highest area of human-related LULC cover compared to CA2 and CA1.
The highest weighted average of TP, Cond, and TSS concerning agriculture, natural, and built covers indicate that the highest values of these WPI measurements occurred in CA with those kinds of predominant cover. According to the same criteria, BOD and TN were the next most common, followed by CF, WTemp, and pH, in terms of occurrence in CA, with a predominance of these land cover types. Water appeared to be less correlated with the CCA axis than with the other cover types, but its correlation with the WPI parameters was similar to that of the other cover types.
CCA for P1/CA1 + P2/CA2 + P3/CA3 + P4/CA4 + P5/CA5 + P6/CA6
Figure 8f displays the resulting ordination diagram of this analysis. The first CCA axis explained 65.7% of the variance, while the second axis explained 27.2%, for a cumulative proportion of 93%. In the ordination diagram, Cond and TSS appeared to be mutually correlated given that they are close to each other. For these six aggregated CAs, the highest WPIs are associated with P5, P4, and P6. These measurement points corresponded to the CAs with the highest human-related LULC cover. The highest weighted average of CF, P, Cond, and TSS concerning agriculture, built, and natural covers indicated that the highest values of these WPI measurements occurred in CA with those kinds of predominant cover. BOD, TN, WTemp, and pH were next in terms of occurrence in CA, with a predominance of these land cover types using the same criteria.
Water cover behaved differently in terms of its relationship with the WPI measurements. The highest TSS, Cond, pH, and TP values (in this order) appeared to occur within CA, which had the highest water areas, whereas the lowest values related to this cover type were DBO, WTemp, TN, and CF.
CCA for P1/CA1 + P2/CA2 + P3/CA3 + P4/CA4 + P5/CA5 + P6/CA6 + P7/CA7
The resulting ordination diagram, shown in Fig. 8g, represents the entire basin (the aggregation of all the CAs). The first CCA axis explained 58.9% of the variance, and the second axis explained 29.2%, for a cumulative proportion of 88.2%. As seen in the ordination diagram, Cond and TSS appeared to be mutually correlated because they are collocated. The diagram also shows a general grouping of P1, P2, and P3 on one side and P4, P5, P6, and P7 on the other. The built, agricultural, and natural cover had similar biplot scores, while the water cover remained separate in the plot. The diagram reveals that the highest WPI values of DBO, TP, CF, Cond, and TSS were correlated with P4, P5, P6, and P7, which corresponded to the most anthropized CAs. The highest weighted averages of DBO, TP, TN, CF, Cond, and TSS (in that order) with respect to agriculture, natural, and built covers indicate that the highest values of these WPI measurements occurred in CA with those kinds of predominant cover.