*Article* **Integrated Approach to Quantify the Impact of Land Use and Land Cover Changes on Water Quality of Surma River, Sylhet, Bangladesh**

**Abdul Kadir <sup>1</sup> , Zia Ahmed <sup>1</sup> , Md. Misbah Uddin <sup>2</sup> , Zhixiao Xie <sup>3</sup> and Pankaj Kumar 4,\***


**Abstract:** This study aims to assess the impacts of land use and land cover (LULC) changes on the water quality of the Surma river in Bangladesh. For this, seasonal water quality changes were assessed in comparison to the LULC changes recorded from 2010 to 2019. Obtained results from this study indicated that pH, electrical conductivity (EC), and total dissolved solids (TDS) concentrations were higher during the dry season, while dissolved oxygen (DO), 5-day biological oxygen demand (BOD<sup>5</sup> ), temperature, total suspended solids (TSS), and total solids (TS) concentrations also changed with the season. The analysis of LULC changes within 1000-m buffer zones around the sampling stations revealed that agricultural and vegetation classes decreased; while built-up, waterbody and barren lands increased. Correlation analyses showed that BOD<sup>5</sup> , temperature, EC, TDS, and TSS had a significant relationship (5% level) with LULC types. The regression result indicated that BOD<sup>5</sup> was sensitive to changing waterbody (predictors, R<sup>2</sup> = 0.645), temperature was sensitive to changing waterbodies and agricultural land (R<sup>2</sup> = 0.889); and EC was sensitive to built-up, vegetation, and barren land (R<sup>2</sup> = 0.833). Waterbody, built-up, and agricultural LULC were predictors for TDS (R<sup>2</sup> = 0.993); and waterbody, built-up, and barren LULC were predictors for TSS (R<sup>2</sup> = 0.922). Built-up areas and waterbodies appeared to have the strongest effect on different water quality parameters. Scientific finding from this study will be vital for decision makers in developing more robust land use management plan at the local level.

**Keywords:** water quality; buffer zone; land use/land cover; Bangladesh

### **1. Introduction**

Water is a vital resource for the maintenance of life, ecological functioning, biological diversity, and social well-being. Despite its importance for life, in recent decades, excessive human land use has severely harmed the quality and quantities of available water resources [1]. In particular, it is well known that rivers function as integrators of land-water connections, receiving pollutants from the surrounding landscapes [2], and river water quality could be negatively impacted. Surma River, an important river in Bangladesh, has been collecting pollutants from a wide variety of point and non-point sources along its course from agricultural wastes, industrial effluents, menage wastes, and municipal sewage [3]. Due to population growth, urbanization, and industrialization, the surrounding landscape of the Surma river has been changing, and the riverside has experienced tremendous development in terms of commercial, human settlement, and industrial development. The population and urban sprawl have adverse effects on the quality of the Surma river water; the increased urban area is responsible for generating large amounts of nonpoint source pollution through runoff and degraded the river water quality.

**Citation:** Kadir, A.; Ahmed, Z.; Uddin, M.M.; Xie, Z.; Kumar, P. Integrated Approach to Quantify the Impact of Land Use and Land Cover Changes on Water Quality of Surma River, Sylhet, Bangladesh. *Water* **2022**, *14*, 17. https://doi.org/10.3390/ w14010017

Academic Editor: František Petroviˇc

Received: 16 November 2021 Accepted: 19 December 2021 Published: 22 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Land use and land cover (LULC) within the immediate environments of a waterbody have a direct impact on the physicochemical and microbiological properties of water, and such impact varies with the type, extent, location of human land uses, and the inputs from the watershed. Water quality problems arise when the type and extent of human land use exceed the natural ability of the watershed to mitigate accumulated land-use-related stress [4]. Numerous studies have shown that human activities lead to landscape pattern changes, which in turn had significant impacts on the conditions of river water [2,5,6]. LULC change, especially urbanization, has a major impact on hydrology, affecting water quality and quantity on a range of spatial and temporal scales [7,8]. Zhu [9] found that water quality degrading was particularly affected by alteration from farmland to commercial and residential land, and the expansion in an urban area causes streamflow increase, carrying more sediment, bank erosion, and nutrients in streams. Hossain [10] reported that water quality variables are correlated with LULC change. As a consequence of the spatiotemporal LULC change, the concentration of diffuse pollutants in streams varies as well as vegetation types, and watershed climate are responsible for stream water quality change.

Landscape pattern has a complex, space- and scale-dependent effect on water quality [11,12], and different landscape characteristics play different roles in receiving water at different spatial scales varying from the local to eco-regional scales [13]. Bhaduri et al. [7] noted that the most significant human impacts on the hydrologic system and water resources are caused by land-use changes on local, regional, and global scales, driven by a rise in urban areas. Pollution from nonpoint sources (NPS) is a challenging problem to solve as it comes from a variety of origins difficult to pinpoint, and it occurs in a variety of environments; however, Geographical Information Systems (GIS) software provides a more comprehensive description of land cover patterns and the spatial distribution of NPS pollution [14] and has been commonly used. To explore the landscape pattern's impacts on the lakes and rivers water quality several studies have been conducted [15,16], and it seems that the riparian buffer zone landscape patterns are more powerful in explaining water quality variations [17]. Land use types have an influence on surface water quality which can be analyzed by using statistical methods, remote sensing (RS), and geographic information system (GIS) [18–21]. Li et al. [12] stated that in riparian zones, the landscape category has a significant impact on water quality, and the alterations in the landscape through urban spread have put a lot of pressure on undeveloped land. Ribeiro et al. [22] found that water quality was worse in the sub-basin, also characterized by the presence of more agriculture, permanent conservation area, lesser natural forest, and a greater drainage area; however, the existence of agriculture negatively affect water quality in the riparian area. As proven by numerous studies, water quality is closely correlated to landscape pattern, which includes the landscape structure and spatial configuration [23–26].

So, it is essential to find out the spatiotemporal and future potential effect on water quality by LULC change. Therefore, this study examines the LULC change in the riparian buffer zone of the Surma river, over the period 2010 to 2019. The seasonal surface water quality variation is analyzed next. The study aims to assess the impacts of land use and land cover (LULC) changes on the Surma river surface water quality.

The outcomes of this study will make it helpful to acquire sustainability of land and water resources. Additionally, it will facilitate other researchers' pursuit of studies about LULC change impact on water quality in the study area and similar regions.

### **2. Materials and Methods**

### *2.1. Design of the Study*

We have conducted this research work systematically and scientifically; the major methods and techniques, which were followed very carefully, are illustrated in Figure 1.

We have conducted this research work systematically and scientifically; the major methods and techniques, which were followed very carefully, are illustrated in Figure 1.

**Figure 1.** Research Design. **Figure 1.** Research Design.

**2. Materials and Methods**  *2.1. Design of the Study* 

#### *2.2. Study Area 2.2. Study Area*

In Bangladesh, Surma is an important river, as a part of the Surma–Meghna river system, which originates when the Barak River in northeastern India and then splits into two branches at the Bangladesh border as the Surma (a northern branch that flows west and then runs southwest to the town of Sylhet) and the Kushiyara rivers [27]. Sylhet Sadar Upazila and Dakshin Surma Upazila are two Upazilas of the Sylhet district, which are located in the country's north-eastern region. Upazila is an administrative region in Bangladesh, equivalent to a county of Western countries. The study was conducted in the Surma river portion, situated between Sylhet Sadar Upazila and Dakshin Surma Upazila land, extending 1000 m toward both sides of the riverbank. It is situated between In Bangladesh, Surma is an important river, as a part of the Surma–Meghna river system, which originates when the Barak River in northeastern India and then splits into two branches at the Bangladesh border as the Surma (a northern branch that flows west and then runs southwest to the town of Sylhet) and the Kushiyara rivers [27]. Sylhet Sadar Upazila and Dakshin Surma Upazila are two Upazilas of the Sylhet district, which are located in the country's north-eastern region. Upazila is an administrative region in Bangladesh, equivalent to a county of Western countries. The study was conducted in the Surma river portion, situated between Sylhet Sadar Upazila and Dakshin Surma Upazila land, extending 1000 m toward both sides of the riverbank. It is situated between 24◦51039.100 N to 24◦54037.0700 N latitude and between 91◦49040.900 E to 91◦55051.900 E longitude (Figure 2).

tude (Figure 2).

tude (Figure 2).

**Figure 2.** The study area and sampling stations. **Figure 2.** The study area and sampling stations. **Figure 2.** The study area and sampling stations.

*2.3. Water Sampling and Analytical Methods 2.3. Water Sampling and Analytical Methods 2.3. Water Sampling and Analytical Methods* 

This study considered both primary and secondary data on water quality. Water quality data for the dry and wet season of 2010 was collected from the Department of Environment (DoE), Sylhet Divisional Office, Sylhet, Bangladesh. This study considered both primary and secondary data on water quality. Water quality data for the dry and wet season of 2010 was collected from the Department of Environment (DoE), Sylhet Divisional Office, Sylhet, Bangladesh. This study considered both primary and secondary data on water quality. Water quality data for the dry and wet season of 2010 was collected from the Department of Environment (DoE), Sylhet Divisional Office, Sylhet, Bangladesh.

24°51′39.1″ N to 24°54′37.07″ N latitude and between 91°49′40.9″ E to 91°55′51.9″ E longi-

24°51′39.1″ N to 24°54′37.07″ N latitude and between 91°49′40.9″ E to 91°55′51.9″ E longi-

*Water* **2022**, *13*, x FOR PEER REVIEW 4 of 16

Water samples were collected from three (ST-1, ST-2, and ST-3) separate Surma river sampling stations, each with a different geographic location, as shown in Table 1. Following the collection, water samples were kept in an ice box and bought to the laboratory. Sampling was performed in the dry season (2019) and the wet season (2019). Water samples were collected from three (ST-1, ST-2, and ST-3) separate Surma river sampling stations, each with a different geographic location, as shown in Table 1. Following the collection, water samples were kept in an ice box and bought to the laboratory. Sampling was performed in the dry season (2019) and the wet season (2019). Water samples were collected from three (ST-1, ST-2, and ST-3) separate Surma river sampling stations, each with a different geographic location, as shown in Table 1. Following the collection, water samples were kept in an ice box and bought to the laboratory. Sampling was performed in the dry season (2019) and the wet season (2019).

**Table 1.** Sampling stations with location**. Table 1.** Sampling stations with location.


day biological oxygen demand (BOD5), pH, electrical conductivity (EC), temperature, total suspended solids (TSS), total dissolved solids (TDS), and total solids (TS). All the water samples of 2019 were analyzed in the Water Supply and Sewerage Engineering Laboratory of Civil and Environmental Engineering Department, and Environmental Laboratory Eight water quality parameters were selected, including dissolved oxygen (DO), 5 day biological oxygen demand (BOD5), pH, electrical conductivity (EC), temperature, total suspended solids (TSS), total dissolved solids (TDS), and total solids (TS). All the water samples of 2019 were analyzed in the Water Supply and Sewerage Engineering Laboratory of Civil and Environmental Engineering Department, and Environmental Laboratory Eight water quality parameters were selected, including dissolved oxygen (DO), 5-day biological oxygen demand (BOD5), pH, electrical conductivity (EC), temperature, total suspended solids (TSS), total dissolved solids (TDS), and total solids (TS). All the water samples of 2019 were analyzed in the Water Supply and Sewerage Engineering Laboratory of Civil and Environmental Engineering Department, and Environmental Laboratory of Geography and Environment Department, Shahjalal University of Science and Technology, Sylhet, Bangladesh. DO and BOD<sup>5</sup> were measured by the Winkler titration method, pH

Eight water quality parameters were selected, including dissolved oxygen (DO), 5-

and temperature were measured by electrometric methods (pH meter, HANNA-HI 9125), electrical conductivity (EC) was measured using an EC meter (HANNA-HI 98192); and TDS, TSS, and TS were measured in the laboratory following standard methods [28].

### 2.3.1. LULC Change Analysis

For land use and land cover (LULC) change analysis, satellite images of 2010 (Landsat 5 Thematic Mapper) and 2019 (Sentinel 2 MSI), acquired from the United States Geological Survey (USGS), were used to produce LULC classification maps for both years using remote sensing (ERDAS IMAGINE 2014) and geographic information system (ArcGIS 10.8) software. LULC classes are categorized into five major categories including, waterbody (river, canals, pond, lakes, reservoirs), built-up (urban areas, human settlements, road networks. Commercial and industrial areas), agricultural land (cropland, pasture, herb, shrub, fallow land, permeable surface), vegetation (canopy, mixed forest, evergreen forest), and barren land (bare soil, sand, rocks without vegetation). Satellite image preprocessing, as well as geometrical rectification, registration of image, corrections viz. atmospheric and radiometric, were conducted by ERDAS IMAGINE 2014. Supervised classification was conducted to create LULC maps [29]. Accuracy assessment was conducted, indicating that the overall classification accuracy of the 2010 image was 81.65% and a kappa statistics of 0.7524, and overall classification accuracy of the 2019 image was 94.50% and kappa statistics of 0.9024, indicating a very good accuracy of the LULC map. By using ArcGIS 10.8 software, land uses composition within the 1000-m buffer zones around the sampling stations was extracted from the LULC map. A 1000-m buffer scale is stronger than smaller scales in explaining land-use types and their water quality relations [30]. Percentages of these broad LULC types were used to examine the relationship between water quality parameters and LULC types.

### 2.3.2. Statistical Analyses

Descriptive statistics were used to explain the general characteristics of LULC and water quality parameters. Karl Pearson's correlation analysis was used to determine correlations between LULC patterns percentage and water quality parameters (WQPs) at statistical significance at a 5% level. Backward stepwise regression analysis was used to identify the relationship between the percentage of land usage composition within the 1000-m buffer zone and water quality properties. WQPs showing significant correlations with LULC types were considered for the backward stepwise regression analysis. In regression analysis, water quality parameters (BOD5, temperature, EC, TDS, and TSS) were considered dependent variables, while LULC types (waterbody, built-up, vegetation, agricultural land, and barren land) were treated as independent variables. To identify the best combination of land uses for water quality estimation regression equations were compared with R<sup>2</sup> values (value closer to one indicates a greater accuracy of the model). The Statistical Packages for Social Science (IBM SPSS Statistics 20) for windows was used to perform all statistical analyses.

### **3. Results**

### *3.1. Water Quality*

The quality of water throughout the dry and wet seasons in 2010 and 2019 are presented in Tables 2 and 3, respectively, and in Figure 3. DO levels ranged from 3.6 mg/L to 4.4 mg/L in the dry season and 7.6 mg/L to 11.6 mg/L in the wet season in 2019. In the wet season, the highest DO was found at Kazir Bazar, while the lowest was found at Keane Bridge in dry season. The DO level increased in all stations during the wet season. In 2010, the DO level varied from 5.1 mg/L in the wet to 6.2 mg/L in the dry season. The Surma river's mean DO level in 2019 was higher than the DO level in 2010.


**Table 2.** Water quality during the dry and wet season in 2010.

**Table 3.** Water quality during the dry and wet season in 2019.


For BOD5, in 2019, the recorded average concentration of BOD<sup>5</sup> for the Surma river water was 0.97 mg/L and 2.67 mg/L in the dry and wet seasons respectively. The seasonal comparison shows that during the wet season the BOD<sup>5</sup> level was higher than the dry season. The lowest and the highest were found at ST-3. In 2010, the average BOD<sup>5</sup> level was 1.2 mg/L (dry season) and 1.1 mg/L (wet season). The mean BOD<sup>5</sup> level slightly increased in 2019 as compared with 2010.

In 2019, the measured pH amongst different stations varied between 6.97 and 7.96 in the dry season, indicating almost neutral-to-slightly-alkaline water conditions. While in the wet season, pH level varied between 6.47 and 7.31. The average dry season pH was slightly alkaline, while the pH of the wet season was almost neutral. In 2010, Surma river water was almost neutral, and pH varied from 7.4 (dry season) to 6.8 (wet season). The mean value of pH was close to neutral in 2010 and 2019.

The average dry season water temperature was 20.3 ◦C in 2010 and 27.77 ◦C in 2019, while the average wet season water temperature was 29.6 ◦C in 2010 and 26.7 ◦C in 2019. Kazir Bazar had the highest dry season temperature in 2019 and Keane Bridge, as well as Kazir Bazar, had the highest wet season temperature in 2010. The recorded lowest wet season temperature was at ST-1 in the year 2010 and the lowest dry season temperature was at ST-1 and ST-2 in the year 2010. The mean temperature in 2019 is higher than the mean temperature in 2010. The temperature increased during the dry season from 2010 to 2019 but declined during the wet season. The maximum and minimum electrical conductivity (EC) of the Surma river was 322.7 µS/cm at Kazir Bazar in the dry season of 2019 and 54.46 µS/cm at Shahjalal Bridge in the wet season of 2019. During the dry season, EC was high. From 2010 to 2019, the mean EC level decreased.

**Figure 3.** Seasonal and spatial variations of WQPs in 2010 and 2019. **Figure 3.** Seasonal and spatial variations of WQPs in 2010 and 2019.

The results of TDS showed that in 2019 the mean TDS concentrations were 149.47 mg/L (dry season) and 34.11 mg/L (wet season). In the dry and wet seasons, maximum TDS was found at Kazir Bazar (ST-3) and minimum TDS was found at ST-1. The mean TDS concentration in 2010 was significantly higher than the mean TDS concentration in 2019. The average maximum (616.67 mg/L) and minimum (49.9 mg/L) TSS concentrations were recorded during the wet and the dry season, respectively, in 2019. In 2019, the mean concentration of TSS was higher than the concentration of the year 2010. From 2010 to 2019, the TSS level increased, on the other hand, the TDS level decreases. In 2019, the recorded total solids (TS) in the dry season varied from 188.7 mg/L at Shahjalal Bridge (ST-1) to 209.9 mg/L at Kazir Bazar (ST-3). The TS ranged from 507.9 mg/L (minimum) at ST-3 to 757.23 mg/L (maximum) at ST-1. The seasonal comparison of the average TS level showed that the wet season's concentration was higher than the dry season's. In 2010, the average TSS concentration varied from 446.6 mg/L during the wet season to 550 mg/L during the dry season with a mean concentration of 498.3 mg/L that was greater than the mean concentration in 2019. The results of TDS showed that in 2019 the mean TDS concentrations were 149.47 mg/L (dry season) and 34.11 mg/L (wet season). In the dry and wet seasons, maximum TDS was found at Kazir Bazar (ST-3) and minimum TDS was found at ST-1. The mean TDS concentration in 2010 was significantly higher than the mean TDS concentration in 2019. The average maximum (616.67 mg/L) and minimum (49.9 mg/L) TSS concentrations were recorded during the wet and the dry season, respectively, in 2019. In 2019, the mean concentration of TSS was higher than the concentration of the year 2010. From 2010 to 2019, the TSS level increased, on the other hand, the TDS level decreases. In 2019, the recorded total solids (TS) in the dry season varied from 188.7 mg/L at Shahjalal Bridge (ST-1) to 209.9 mg/L at Kazir Bazar (ST-3). The TS ranged from 507.9 mg/L (minimum) at ST-3 to 757.23 mg/L (maximum) at ST-1. The seasonal comparison of the average TS level showed that the wet season's concentration was higher than the dry season's. In 2010, the average TSS concentration varied from 446.6 mg/L during the wet season to 550 mg/L during the dry season with a mean concentration of 498.3 mg/L that was greater than the mean concentration in 2019.

### *3.2. Descriptive Statistics of the Water Quality Parameters (WQPs)*

*Water* **2022**, *13*, x FOR PEER REVIEW 8 of 16

The descriptive statistics of the WQPs of the Surma River from the year 2010 to 2019 are explained in Table 4. *3.2. Descriptive Statistics of the Water Quality Parameters (WQPs)*  The descriptive statistics of the WQPs of the Surma River from the year 2010 to 2019


**Table 4.** Descriptive statistics of the WQPs. are explained in Table 4.

### *3.3. LULC Change in the Buffer Zone 3.3. LULC Change in the Buffer Zone*

Land use change data, as the percentage of land of 1000-m buffer zone of monitoring stations, were derived from LULC change analysis from 2010 to 2019, and they were linked with water quality data. The results show that, in 2019, the central part was dominated by built-up area, while in 2010, built-up area was not clustered at a specific zone (Figure 4). Land use change data, as the percentage of land of 1000-m buffer zone of monitoring stations, were derived from LULC change analysis from 2010 to 2019, and they were linked with water quality data. The results show that, in 2019, the central part was dominated by built-up area, while in 2010, built-up area was not clustered at a specific zone (Figure 4).

**Figure 4.** LULC change during the period 2010–2019. **Figure 4.** LULC change during the period 2010–2019.

The LULC change in the buffer zones showed that the agricultural land area decreased while the built-up area significantly increased from 2010 to 2019. During this period, the increase of barren land and decrease of vegetation area was also observed in the The LULC change in the buffer zones showed that the agricultural land area decreased while the built-up area significantly increased from 2010 to 2019. During this period, the

increase of barren land and decrease of vegetation area was also observed in the buffer zones. In 2019 the highest built-up area was found at the buffer zone of station ST-2.

The analysis illustrated that in 2010, at ST-1 (32.4%) and ST-2 (39.9%), agricultural land usage was dominant in the 1000-m buffer zones (Table 5). On the other hand, vegetation cover (40.6%) was the dominant land area in 2010 at ST-3. In 2019, the built-up area increased in all zones and become the dominant land-use type.


**Table 5.** LULC change within the 1000-m buffer zone of sampling stations.

By comparing land usage distribution from 2010 to 2019, all zones had shown a notable decline in agricultural land area. For the built-up area, the ST-2 zone showed a maximum increase (22.8%), while the ST-1 zone showed a minimum increase (11.7%), and at the ST-3zone, 21.1% area increased.

### *3.4. Descriptive Statistics of the LULC Types*

The descriptive statistics of LULC types within a 1000-m buffer zone at stations of the Surma River area from 2010 to 2019 are reported in Table 6.


**Table 6.** Descriptive statistics of the LULC types.

In 2010, waterbody ranged from 5.39% to 6.28% with a mean value of 5.84% ± 0.45%. In 2019, the range and mean of waterbody are 6.91% to 7.5% and 7.15% ± 0.31% respectively. The built-up area, from 2010 to 2019, ranged from 18.71% to 30.12% and from 38.81% to 48.57%, with a mean value of 24.89% ± 5.76% and 43.07% ± 5.00%, respectively. Agricultural land use ranged from 32.28% to 39.87% in 2010, with a mean value of 34.84% ± 4.36%. On the other hand, in 2019, agricultural land use ranged from 18.06% to 23.83% with a mean value of 20.07% ± 3.26%.

Vegetation area in 2010 and 2019 ranged from 26.72% to 40.62% and from 20.67% to 32.26% with a mean value of 31.73% ± 7.72% and 25.1% ± 6.26%, respectively. Barren land ranged from 2.17% to 3.36% in 2010, with a mean value of 2.7% ± 0.61%. On the other hand, in 2019, barren land cover ranged from a minimum of 3.58% to a maximum of 6.18% with a mean value of 4.62% ± 1.37%.

### *3.5. Water Quality Relation with LULC* Correlation Analysis

A correlation analysis revealed that LULC patterns were correlated significantly with one or more parameters of water quality within the 1000-m buffer zone scale (Table 7). The analysis result reveals that only BOD5, temperature, EC, TDS, and TSS had a strong and significant relationship at the 5% level of significance with the different LULC types. Other parameters viz. DO, pH, and TS had both positive and negative relationships with the LULC types but were not significant. BOD<sup>5</sup> was found to have a positive significant relationship with waterbody (0.82), which means that, if the waterbody should increase, the BOD<sup>5</sup> level will also increase at a significant level. At the same time, it was also found that temperature had a positive significant relationship with waterbody (0.82) but it had a negatively significant relationship with agricultural land (−0.90). This result reveals that, with the increase in agricultural land area, temperature will decrease at a significant level or vice-versa.


**Table 7.** Linear relationship (Pearson correlation, r) between WQPs and LULC types.

NB: Bold letters indicates a significant relationship at the 5% level.

Electric conductivity (EC) had about a perfect positive significant relationship with vegetation (0.92) but a negative relationship with the land use types of built-up area (−0.82) and barren land (−0.84). Total dissolved solids (TDS) had also a significant relationship at the 5% level of significance, with three types of LULC types whereas total suspended solid (TSS) had four types of LULC types. Pearson's correlation result revealed that with the increase of waterbody and built-up area TSS value will increase significantly (a positive relationship) but the TDS value will decrease significantly (a negative relationship). A different result was also found for agricultural land; agricultural land had inverse relationships with the changes in TDS and TSS. The result reveals that, if agricultural land increases, TDS (0.92) will also increase significantly at the 5% level of significance, but TSS (−0.83) will decrease.

Vegetation cover was only correlated with the change of electric conductivity, which had about a perfect positive significant correlation.

### *3.6. Fulfillment of Distributional Assumption of Dependent Variables*

The median of BOD<sup>5</sup> was found at 1.38, which was very similar to the average BOD<sup>5</sup> (1.48). The similarity between mean BOD<sup>5</sup> and median BOD<sup>5</sup> seems to follow the symmetrical distribution of BOD5, which means that the dependent variable, BOD5, fulfills the distributional assumption for the classical model.

The similarity between the mean (26.12) and median (26.30) temperatures indicates that temperature seems to follow the normal distribution (Figure 5). There were no observed outliers in temperature, thus we can proceed to the classical model with it. The variable EC seems to follow a bell-shaped distribution, because its second quartile (194.63) lay

in the middle of the first (185.25) and third (203.56) quartile. With the fulfillment of the distributional assumption of EC, we can proceed to a classical model. *Water* **2022**, *13*, x FOR PEER REVIEW 11 of 16

**Figure 5.** Box and whisker plots of dependent variables. **Figure 5.** Box and whisker plots of dependent variables.

*3.7. Regression Analysis* 

The median of TDS was found as 224.83 and lay in the middle of the first quartile (93.88) and third quartile (391.25). That's why the small departure of median TDS from mean TDS does not affect the symmetrical shape of TDS. Thus, we can proceed to a classical model for TDS. The median of TDS was found as 224.83 and lay in the middle of the first quartile (93.88) and third quartile (391.25). That's why the small departure of median TDS from mean TDS does not affect the symmetrical shape of TDS. Thus, we can proceed to a classical model for TDS.

The small departure of median TSS from the middle position of the box may have a small effect, an asymmetric shape, and also there was a difference between the mean TSS (220.63) and the median TSS (189.48). This means that TSS may follow a slightly skewed distribution. So, we can, but barely, proceed to the classical model with TSS. The small departure of median TSS from the middle position of the box may have a small effect, an asymmetric shape, and also there was a difference between the mean TSS (220.63) and the median TSS (189.48). This means that TSS may follow a slightly skewed distribution. So, we can, but barely, proceed to the classical model with TSS.

Backward stepwise regression identified the relationship between WQP and LULC types that determine the combination of land uses for water quality estimation (Table 8). For the case of BOD5, only the waterbody was used as a predictor, which was found significant in a Pearson's correlation analysis. Waterbody was not found as an especially

### *3.7. Regression Analysis*

Backward stepwise regression identified the relationship between WQP and LULC types that determine the combination of land uses for water quality estimation (Table 8). For the case of BOD5, only the waterbody was used as a predictor, which was found significant in a Pearson's correlation analysis. Waterbody was not found as an especially strong predictor, as the adjusted r-squared value was 0.645. For temperature, waterbody and agricultural land were used as predictors (R<sup>2</sup> = 0.889). Similarly, for electric conductivity, built-up, vegetation, and barren land were used as predictors (R<sup>2</sup> = 0.833). For total dissolved solid, waterbody, built-up, and agricultural land were used as predictors (R<sup>2</sup> = 0.993). For TSS, waterbody, built-up, and barren land were used as predictors at R<sup>2</sup> of 0.922. From the regression analysis, it could be concluded that BOD<sup>5</sup> showed sensitivity to changing waterbody, whereas temperature was sensitive to changing waterbody and agricultural land. EC showed sensitivity on built-up, vegetation, and barren land. TDS showed sensitivity on the waterbody, built-up and agricultural land whereas, TSS for waterbody, built-up and barren land. In this study, different parameters of water quality viz. biological oxygen demand, electric conductivity, TDS, and TSS, tended to be most affected by built-up and waterbody land usage types.

**Table 8.** Linear regression models of LULC types on the WQPs.


Y\_2019 = year 2019 (dummy or indicator); Temp = temperature; W = waterbody, A = agricultural land; Bu = builtup; V = vegetation; Ba = barren land.

The equation for BOD<sup>5</sup> explains, that for a one-unit change of waterbody, BOD<sup>5</sup> would increase 0.068 times. For temperature, it would decrease, both for waterbody and agricultural land. In the same way, the equation for electric conductivity explains that, if the built-up area increases by one unit, EC will decrease 0.802 times, and for barren land, it will decrease 7.369 times. Yet, EC could increase 1.332 times if vegetation covers were increased by one unit. For TDS and TSS, the changes were the same, as is quantified and described in their equations.

### **4. Discussion**

The water quality of any river is sturdily influenced by landscape characteristics, including land use land cover types and their spatial patterns [31]. LULC change mainly depends on how humans alter the natural landscapes and socioeconomic growth through space and time. Land covers involve the physical features of the Earth's surface that are occupied by vegetation, water, soil, and other characteristics of the Earth's surface created through human activities, whereas land used by human beings for habitats concerning economic activities is referred to as land cover [32].

River water qualities near the urban area changed due to several factors and LULC change was the most significant among them. LULC change has had a major impact on water quality. In the urbanization process, the built-up area increases rapidly. As a result, the quality of surrounding river water deteriorates. The rapid expansion of human

settlements in the urbanized area discharge huge volumes of sewage water, as the point source of pollution, which includes a high level of nutrients and metals. Therefore, water physicochemical parameters such as pH, dissolved oxygen, biochemical oxygen demand, and chemical oxygen demand have significantly changed. The impact of land-use changes on water quality is usually studied by analyzing the relationships between land use and water quality parameters. Water quality differs according to location, time, weather, and pollution sources [33,34], and contamination is generally determined by studying the physical and chemical properties of the water bodies [35].

Findings from this study indicate that LULC types are considerably associated with one or more water quality parameters in the 1000-m buffer zone scale. It also reveals that only BOD, temperature, EC, TDS, and TSS had a strong and significant relationship with the different LULC types. In comparison, Kerala observed the water quality parameters of the Chalakudy river and compared them with diverse land use patterns over four seasons [36]. They found that urban land use was associated with poor water quality during the study period when there were changes in land use and land cover patterns [37]. Moreover, land use changes in the surrounding area of cities can modify the surface properties of watersheds that influence runoff quality and quantity. The impact of LULC changes on water quality involves analyzing the relationship between land use and water quality indicators [38]. In this research, among five types of land use and land cover, waterbody and built-up were the most significant variables in predicting water quality parameters. The Pearson correlation suggested that BOD is positively correlated with waterbody and built-up area but negatively correlated with agricultural land. A similar study by Tong and Chen [39] observed the water quality of a watershed in relation to land use change in Ohio State, USA, and their results indicate that BOD was positively correlated with residential and commercial lands but had only a non-significant correlation with agricultural land. Xiao et al. [40] conducted a multi-scale analysis of the relationship between urban river water qualities with landscape patterns in different seasons in Huzhou City, China, and their findings point out that, at a different scale, their relationships varied with the composition of land-use types—but, built-up land was most significant. These results suggest that with the development of different types of built-up land, water quality parameters change and exacerbate contamination.

Regression analysis of the present study showed that BOD was sensitive to changing waterbody, whereas temperature was sensitive to changing waterbody and agricultural land. Urban land is a mixture of different land uses types, such as residential, industrial, commercial, and other built-up areas. In contrast with other land use, wastewater is generated more in urban areas, and also urbanization increased coverage by impervious surfaces, which influences storm flow speed and runoff volumes [41]. Runoff and the huge volume of storm and drainage water are mainly responsible for this relationship between BOD and waterbody sensitivity because a higher BOD value indicates that a greater amount of organic matter is present; thus, storm flows and drainage water deliver more pollutants to the surrounding urban catchment, especially in a river. Due to the huge volume of drainage and stormwater, the pollutant quantity increased for this reason in the wet seasons, while BOD was high at a few stations.

The present study's Pearson's correlation results showed that, with the increase of waterbody and built-up, TSS will increase significantly (a positive relationship) and backward stepwise regression identified a relationship between WQP and LULC types, indicating that TDS were sensitive to waterbody, built-up, and agricultural land, whereas TSS was found sensitive for waterbody, built-up and barren land. Wang and Zhang [42] studied the relationship between landscape types and water quality index (WQI), in a multi-scale analysis in the Ebinur Lake oasis; their findings revealed that, for different buffers, both positive and negative relationships exist between certain land use and land cover (LULC) types and the water quality index, but there was considerable correlation between water quality index and landscape index. Li et al. [43] found a relationship between land use/cover and water quality using correlation and regression analyses in the Liao River basin, China, also

indicating that BOD5, COD, sediment, and hardness were considerably associated with land use.

Although, in reviewing the literature, we found similar studies from different parts of the world, but, considering that the present work is the first such kind of investigation for a data-scarce region such as our study area, these results hold a lot of merit to both scientific communities as well as policy makers. This scientific evidence will lay a foundation for designing more robust adaptation and mitigation measures for water resource management in a timely manner.

### **5. Conclusions**

The study revealed the concentrations of several physicochemical parameters of river water in relation to LULC changes. However, the result of the present investigation indicates that LULC changes and seasonal variations (influences the concentration) have a significant impact on water quality parameters. The results of the LULC change analysis indicate builtup, waterbody, and barren land increased and agricultural land and vegetation decreased. Built-up area is dominant in LULC types and the change of LULC pattern within 1000 m buffer zones had a significant impact on the water quality parameters of the Surma river. LULC information, in relation to the surrounding river water quality in urban areas, is very important for planning, monitoring, and management of river water because, in urbanized and densely populated cities, river water is also used for drinking purposes, after treatment, and also for recreational purposes. LULC change causes severe environmental problems worldwide and poses a threat to water quality. Spatiotemporal information about LULC change patterns with water quality helps in finding a solution to this problem. The study provides useful tools for future study, which, combined with LULC change and its relationship with different water quality parameters, can help decision-makers in formulating of rules and guidelines about sustainable land use, especially in city areas, and aid in minimizing negative impacts on water quality.

**Author Contributions:** Conceptualization: A.K., Z.A., M.M.U., Z.X., P.K.; methodology: A.K., Z.A.; formal analysis: A.K., Z.A.; writing—original draft preparation: A.K., Z.A., M.M.U., Z.X., P.K.; writing—review and editing: A.K., Z.A., M.M.U., Z.X., P.K.; funding acquisition: P.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** Research work was supported by the Ministry of Science and Technology, Bangladesh, under the National Science and Technology (NST) Fellowship (FY 2018-19). Also, the publication is supported by the Asia-Pacific Network for Global Change Research (APN) under Collaborative Regional Research Programme (CRRP) with project reference number CRRP2019-01MY-Kumar.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Cheng-Wei Hung <sup>1</sup> and Lin-Han Chiang Hsieh 2,\***


**Abstract:** Eutrophication is an environmental pollution problem that occurs in natural water bodies. Regression analyses with interaction terms are carried out to identify the factors influencing the Shimen, Mingde, and Fongshan Reservoirs in Taiwan. The results indicate that the main factor influencing these reservoirs is total phosphorus. In the Shimen and Mingde Reservoirs, the influence of total phosphorus, when interacting with other factors, on water quality trophic state is more serious than that of total phosphorus per se. This implies that the actual influence of total phosphorus on the eutrophic condition could be underestimated. Furthermore, there was no deterministic causality between climate and water quality variables. In addition, time lagged effects, or the influence of their interaction with other variables, were considered separately in this study to further determine the actual relationships between water trophic state and influencing factors. The influencing patterns for three reservoirs are different, because the type, size, and background environment of each reservoir are different. This is as expected, since it is difficult to predict eutrophication in reservoirs with a universal index or equation. However, the multiple linear regression model used in this study could be a suitable quick-to-use, case-by-case model option for this problem.

**Keywords:** eutrophication; variable interactions; multiple linear regression; reservoir

### **1. Introduction**

Constructing reservoirs is one of the most effective ways of storing water in Taiwan. Surface runoff may be caught during periods of high flow and provide water for people's livelihood, industry, and agriculture during periods of water shortage.

Eutrophication, the most challenging water pollution problem in water bodies, will eventually become an issue in many reservoirs [1]. Eutrophication negatively affects the water quality, safety, ecological integrity, and sustainability of global water resources [2–4]. It has long been believed that excessive phosphorus is the main reason for eutrophication [5]. However, population density, urbanization, and agricultural activities are also factors that influence the water quality of freshwater systems [6–10]. Since the 1940s, a substantial population increase, land-use intensification, and the use of agricultural fertilizers from developed countries [11], as well as the use of detergents containing phosphate compounds since the 1950s, have accelerated the eutrophication of waterbodies [12].

Eutrophication influences the water volume and quality in reservoirs. Regarding water volume, algae distributed on the water surface causes water hypoxia and a decrease in water transparency [13,14], leading to a substantial death of aquatic organisms, which are then deposited on the bottom of the reservoir, which, in turn, reduces the reservoir's capacity over time. Regarding water quality, the proliferation of algae causes algal blooms and releases algal poison, which both influences water quality conditions such as dissolved oxygen, transparency, odor, and pH value, and it also causes problems during the filtration of drinking water, increasing health risks to people [15,16]. Eutrophication also causes the

**Citation:** Hung, C.-W.; Chiang Hsieh, L.-H. Analysis of Factors Influencing the Trophic State of Drinking Water Reservoirs in Taiwan. *Water* **2021**, *13*, 3228. https://doi.org/10.3390/ w13223228

Academic Editors: Pankaj Kumar and Ram Avtar

Received: 13 October 2021 Accepted: 11 November 2021 Published: 14 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

depletion of dissolved oxygen in water, which may have potentially harmful influences on the habitats of fish and macroinvertebrates [17,18]. In addition, eutrophication results in the release of nutrients such as phosphorus and ammonium into the water [19–21], potentially producing toxic heavy metal ions [22].

At present, reservoirs in Taiwan are facing a crisis of deteriorating water quality. According to the Environmental Water Quality Monitoring Annual Report of the Environmental Protection Agency, among the 26 reservoirs on the main island of Taiwan, 7 reservoirs are in a eutrophic state, 18 are in a mesotrophic state, and only 1 is in an oligotrophic state.

However, the Carlson trophic state index (CTSI) for assessing the water quality of reservoirs is not completely suitable in Taiwan. Taiwan is affected by frequent heavy rainfall, and particularly after typhoons and heavy rains, large amounts of sand, soil, and rock flow into the reservoirs, leading to a large increase in the concentration of suspended solids and a decrease in water transparency. Thus, it is quite likely that a water body with a high CTSI value, but without a large amount of algae distributed on the water surface, may be identified as eutrophic.

This study aims to investigate the correlation between weather and water quality factors and their degree of influence on the trophic state of reservoirs both as single variables and as interactions of variables. We also discuss the suitability of CTSI for assessing water quality in reservoirs in Taiwan. We used weather and water quality data from 2017 to 2019 from three main reservoirs in Taiwan: Shimen Reservoir, Mingde Reservoir, and Fongshan Reservoir. Chlorophyll a was used as an indicator to illustrate the degree of eutrophication, and data were analyzed using multiple linear regressions (MLR) including time lags and variable interactions.

### **2. Materials and Methods**

### *2.1. Characteristics of Reservoir*

The Shimen Reservoir is a stable source of water supply in northern Taiwan and is the third largest reservoir in Taiwan (Figure 1). It is a multi-objective water conservancy project that combines benefits such as irrigation, power generation, water supply, flood control, sightseeing, and recreation. The Shimen Reservoir has a catchment area of 763 km<sup>2</sup> , a full water level of 8 km<sup>2</sup> , a total storage capacity of 309 million m<sup>3</sup> , and an effective storage capacity of 197 million m<sup>3</sup> [23]. The irrigation area of the Shimen Reservoir includes three counties: Hsinchu, Taoyuan, and New Taipei City. It provides for a daily consumption of 800,000 m<sup>3</sup> of livelihood water and also the Shimen Power Plant with 230 million kWh of power generation per year [24]. *Water* **2021**, *13*, x FOR PEER REVIEW 3 of 15

**Figure 1.** Location of Shimen Reservoir (including water catchment area). **Figure 1.** Location of Shimen Reservoir (including water catchment area).

**Figure 2.** Location of Mingde Reservoir (including water catchment area).

**Figure 3.** Location of Fongshan Reservoir (off-stream).

The Mingde Reservoir provides more water for consumption in the Miaoli County, where there are more mountains and few fields (Figure 2). The Mingde Reservoir has a catchment area of 61 km<sup>2</sup> , a full water level of 1.7 km<sup>2</sup> , a total storage capacity of 17.7 million m<sup>3</sup> , and an effective storage capacity of 12.2 million m<sup>3</sup> . The irrigation area of the Mingde Reservoir is 13 km<sup>2</sup> , and it provides 27,000 m<sup>3</sup> of water daily [25]. **Figure 1.** Location of Shimen Reservoir (including water catchment area).

*Water* **2021**, *13*, x FOR PEER REVIEW 3 of 15

*Water* **2021**, *13*, x FOR PEER REVIEW 3 of 15

**Figure 2.** Location of Mingde Reservoir (including water catchment area). **Figure 2.** Location of Mingde Reservoir (including water catchment area).

The Fongshan Reservoir is an off-site reservoir located in Kaohsiung City and provides Kaohsiung with a large additional water supply because of the high population concentration and the rapid development of industry and commerce, which increases water consumption in the area (Figure 3). It has a catchment area of 2.75 km<sup>2</sup> , a full water level of 0.75 km<sup>2</sup> , a total storage capacity of 9.2 million m<sup>3</sup> , and an effective storage capacity of 8.5 million m<sup>3</sup> [26]. The Fongshan Reservoir supplies 1.6 million tons of water daily, of which 350,000 tons, 22%, caters for industrial consumption [27]. **Figure 2.** Location of Mingde Reservoir (including water catchment area).

**Figure 3.** Location of Fongshan Reservoir (off-stream). **Figure 3.** Location of Fongshan Reservoir (off-stream).

### *2.2. Dataset*

The Taiwan Environmental Protection Administration (EPA) changed reservoir water quality monitoring from quarterly monitoring in previous years to monthly monitoring from January 2017. In this study, we used monthly weather and water quality data from the Shimen, Mingde, and Fongshan Reservoirs from January 2017 to December 2019.

Weather data included the daily statistics of rainfall (mm) and inflow (mm) data in the catchment area from 2017 to 2019 and was downloaded from the disaster prevention information service website of the Water Resources Agency (WRA). The data collection period corresponds to that of the EPA water quality monitoring data from each reservoir. In addition, monthly water temperature (WT) data from 2017 to 2019 were collected from the national water quality monitoring information website of the EPA [28].

Water quality data were collected from the monthly statistics on the national water quality monitoring information website of the EPA and included chlorophyll a (Chl-a), dissolved oxygen (DO), transparency (SD), total phosphorus (TP), pH, conductivity, suspended solids (SS), chemical oxygen demand (COD), and ammonia nitrogen (AN) sampled from 2017 to 2019 [28].

In order to investigate whether the different seasons affect the degree of influence, March to May were designated as Season 1 (S1), June to August were designated as Season 2 (S2), September to November were designated as Season 3 (S3), and December to February were designated as Season 4 (S4). S4 was used as a control, and S1 to S3 were analyzed to identify the degree of influence that each variable has on water quality in the different seasons.

### *2.3. Methodology*

### 2.3.1. Regression Analysis

We used multiple linear regression (MLR) using regression analysis (Equation (1)). In order to identify how much each factor influences the eutrophication of the reservoirs, we used rapidly adjusting variables in regression models to analyze weather and water quality factors.

$$Y\_{\bar{i}} = \beta\_0 + \beta\_1 X\_{1\bar{i}} + \beta\_2 X\_{2\bar{i}} + \dots + \beta\_k X\_{k\bar{i}} + \varepsilon\_{\bar{i}} \tag{1}$$

where *Y<sup>i</sup>* is the *i*-th observation value in the dependent variable, which represents the concentration of chlorophyll a in this study. The independent variables *X*1*<sup>i</sup>* to *Xki* are weather (rainfall, inflow, and water temperature) and water quality (Chl-a, DO, etc.) factors. *β*<sup>0</sup> is an intercept term, *β*<sup>1</sup> to *β<sup>k</sup>* are the slope terms, and also the unknown coefficients corresponding to the independent variables *X*1*<sup>i</sup>* to *Xki*. *ε<sup>i</sup>* is a random error term.

### 2.3.2. Time-Lag

The reason of applying time-lag variables in this study is that the current weather or water quality factors do not necessarily have an immediate influence on the Chl-a concentration. These time lag situations might be one week, one month, or even more, so it is not suitable to use the monitoring current data in the regression analysis.

Basic analysis, Lag 1, and Lag 2 data in the regression model were selected using the following steps: Firstly, we select the significant data of the three types of data. Secondly, if more than one of the three types of data were significant, we selected the data collected closest to the monitoring date, i.e., basic analysis data take precedence over Lag 1 data, which takes precedence over Lag 2 data. Thirdly, if there was no significance in the three types of data, basic analysis data were selected as a representative term.

### 2.3.3. Interaction Terms

We also test whether the influence of weather and water quality factors on the water trophic state is affected by their interaction. The traditional ordinary least square (OLS) formula can be illustrated as shown in Equation (2):

$$\mathcal{Y} = \beta\_1 \mathcal{X}\_1 + \dots + \beta\_4 \mathcal{X}\_4 \tag{2}$$

Regarding the interaction of factors, we need to consider whether factors are correlated. We tested the pair-by-pair interaction of factors that are theoretically related using an MLR equation and then analyzed whether the correlation was still significant after the interaction. If it was significant, the two factors were grouped to form an interaction term that was added in the equation, and the correlation was analyzed as shown in Equations (3) and (4). *β*5*X*1*X*<sup>2</sup> and *β*6*X*3*X*<sup>4</sup> are the interaction terms added on the basis of an MLR.

$$Y\_1 = \beta\_1 X\_1 + \dots + \beta\_4 X\_4 + \beta\_5 X\_1 X\_2 \tag{3}$$

$$Y\_2 = \beta\_1 X\_1 + \dots + \beta\_4 X\_4 + \beta\_6 X\_3 X\_4 \tag{4}$$

If the individual analysis results of interaction terms *β*5*X*1*X*<sup>2</sup> and *β*6*X*3*X*<sup>4</sup> in the regressions of *Y*<sup>1</sup> and *Y*<sup>2</sup> were both significant, then the two factors were analyzed in the same regression formula for *Y*3, as shown in Equation (5).

$$Y\_3 = \beta\_1 X\_1 + \dots + \beta\_4 X\_4 + \beta\_5 X\_1 X\_2 + \beta\_6 X\_3 X\_4 \tag{5}$$

If the regression results of *β*5*X*1*X*<sup>2</sup> and *β*6*X*3*X*<sup>4</sup> in *Y*<sup>3</sup> were both significantly correlated, the combination was retained, and the method was repeated on other interaction terms. At most, two interaction groups were included in the regression formula of each reservoir, and the factors in the two groups did not overlap with each other.

We used three conditions for selecting two groups of interaction terms with simultaneous significant correlations: Firstly, we considered TP and AN that represent the nutrient factors in the two interaction groups. Secondly, we considered WT, rainfall, and inflow that are representative of the weather factors, and WT had priority over rainfall, and rainfall had priority over inflow. Thirdly, we considered the R<sup>2</sup> value that could be explained by applying each group to the regression formula as a final selection step.

### 2.3.4. Interrelationship of Interaction Terms

Equation (2) was simplified by reducing variables and adding an interaction term, as shown in Equation (6).

$$Y = \beta\_1 X\_1 + \beta\_2 X\_2 + \beta\_3 X\_1 X\_2 \tag{6}$$

Equation (6) was rewritten as Equation (7), in which other variables are fixed, *X*<sup>1</sup> increases by 1 unit, and *X*<sup>2</sup> does not increase, and the dependent variable becomes *Y*1.

$$Y\_1 = \beta\_1(X\_1 + 1) + \beta\_2 X\_2 + \beta\_3(X\_1 + 1)X\_2 \tag{7}$$

Equation (7) minus Equation (6) provides Equation (8), which represents a situation in which other variables are fixed, *X*<sup>1</sup> increases by 1 unit and *X*<sup>2</sup> does not increase, and the unit amount ∆*Y*<sup>1</sup> is the dependent variable that will increase.

$$
\Delta Y\_1 = Y\_1 - Y = \beta\_1 + \beta\_3 X\_2 \tag{8}
$$

Similarly, if other variables are fixed but *X*<sup>1</sup> does not increase and *X*<sup>2</sup> increases by 1 unit, the dependent variable becomes *Y*<sup>2</sup> (Equation (9)). Equation (9) minus Equation (6) provides Equation (10), which represents the unit amount ∆*Y*2—the dependent variable that will increase.

$$Y\_2 = \beta\_2 + \beta\_3 X\_1 \tag{9}$$

$$
\Delta \mathbf{Y}\_2 = \mathbf{Y}\_2 - \mathbf{Y} = \beta\_2 + \beta\_3 \mathbf{X}\_1 \tag{10}
$$

When both *X*<sup>1</sup> and *X*<sup>2</sup> increase by 1 unit, and other variables are fixed, Equation (11) is obtained. Equation (11) minus Equation (6) provides Equation (12), which represents the unit amount ∆*Y*3—the dependent variable that will increase.

$$Y\_3 = \beta\_1(X\_1 + 1) + \beta\_2(X\_2 + 1) + \beta\_3(X\_1 + 1)(X\_2 + 1)\tag{11}$$

$$
\Delta Y\_3 = Y\_3 - Y = \beta\_1 + \beta\_2 + \beta\_3 + \beta\_3 (X\_1 + X\_2) \tag{12}
$$

Equations (8), (10), and (12) were integrated (as shown in Table 1), where ∆*X*<sup>1</sup> and ∆*X*<sup>2</sup> are the unit amounts by which the independent variables, *X*<sup>1</sup> and *X*2, are increased. Table 1 shows that when *X*<sup>1</sup> and *X*<sup>2</sup> are increased by 1 unit separately or simultaneously, the unit amount ∆*Y*, which is the dependent variable, will increase.

**Table 1.** Interrelationship of interaction terms.


### *2.4. Model*

In this study, Chl-a is set as the dependent variable, and the eleven weather and water quality factors, WT, DO, SD, TP, pH, conductivity, SS, COD, AN, rainfall, and inflow, are set as independent variables. STATA (version 13) is used in this study to perform correlation analysis and regression analysis.

The research model of this study is based on MLR. Using the result of a Hausman Test, we chose to include a random effects model in the MLR model to avoid or reduce ignoring differences between the data (which results in the omission of variables and leads to estimation errors) and also to reduce the occurrence of collinearity problems between variables.

The MLR model built in this study includes the random effects model and performs a time-lag analysis as well as adding specific interaction terms separately based on the difference in the reservoir data.

### **3. Results and Discussion**

*3.1. Descriptive Statistics*

This study lists the descriptive statistics of the Shimen, Mingde, and Fongshan Reservoirs (Table 2). The minimums of TP, SS, COD, and AN were below the detection limit (ND).


**Table 2.** Descriptive statistics from the Shimen, Mingde, and Fongshan Reservoir data.

SR: Shimen Reservoir, MR: Mingde Reservoir, FR: Fongshan Reservoir, ND: The monitoring data is smaller than the detection limit

The average values of Chl-a, TP, Conductivity, SS, COD, and AN in the Fongshan Reservoir are several times more than the average values of the other two reservoirs. The concentration of Chl-a is 12 times that of the Shimen Reservoir and 2.4 times that of the Mingde Reservoir. The concentration of TP is 40 times that of the Shimen and Mingde Reservoirs. The concentration of AN is 52.5 times that of the Shimen Reservoir and 56.2 times that of the Mingde Reservoir.

The maximum Chl-a, TP, Conductivity, COD, and AN were recorded in the Fongshan Reservoir, the maximum SS was recorded in the Shimen Reservoir, and the maximum pH was recorded in the Mingde Reservoir. There were large variations in the daily rainfall and inflow data. It is possible that there was no rain on the day monitoring was conducted but heavy rainfalls the next day; therefore, we do not discuss the maximum and minimum rainfall and inflow.

The minimum DO of the Fongshan Reservoir was about three times lower than those of the Shimen and Mingde Reservoirs. However, its minimum TP concentration was still much larger than that of the Shimen and Mingde Reservoirs. The Fongshan Reservoir was the only reservoir with a pH of less than 7. The minimum conductivity of the Fongshan Reservoir was more than two times that of the Shimen and Mingde Reservoirs. The minimum SS of the Fongshan Reservoir was about two times that of the Mingde Reservoir.

### *3.2. Correlation Analysis*

Results from the correlation analysis of the Shimen, Mingde, and Fongshan Reservoirs are show in Tables 3 and 4. The absolute value of the correlation coefficient was above 0.6, which can be regarded as a high correlation between the two factors.


**Table 3.** Correlation coefficients between factors from correlation analyses (1/2).

**Table 4.** Correlation coefficients between factors from correlation analyses (2/2).


In the Shimen Reservoir, WT and pH were highly correlated with a correlation coefficient of 0.700 (Table 3).

In the Mingde Reservoir, the pairs of WT and conductivity, DO and pH, pH and conductivity, and rainfall and inflow were highly correlated with correlation coefficients of −0.644, 0.634, −0.649, and 0.826 respectively (Tables 3 and 4).

In the Fongshan Reservoir, the pairs of Chl-a and COD, DO and pH, SD and SS, TP and Conductivity, TP and AN, and conductivity and AN were highly correlated with correlation coefficients of 0.602, 0.700, −0.600, 0.854, 0.827, and 0.696, respectively (Tables 3 and 4).

The correlation coefficient of rainfall and inflow in the Mingde Reservoir as well as the correlation coefficient of TP with conductivity and TP with AN in the Fongshan Reservoir exceeded 0.8. However, when performing an MLR analysis, a high correlation coefficient between the independent variables causes the problem of collinearity in the regression results. This problem means that the lower the correlation between the independent variables, the more it reflects the relationship with dependent variables.

*3.3. Analysis Result*

Tables 5–7 show the regression analysis result after including time-lag and interaction terms from the Shimen, Mingde, and Fongshan Reservoir data, respectively. This study divides independent variables into "Basic" analysis without time lag, "Lag 1" data with one month lagged, and "Lag 2" data with two months lagged. For example, there is a time lag between an increase in water temperature, the growth of algae, the flow time of rainfall runoff to the reservoir, etc.


**Table 5.** Results from regression analyses of the Shimen Reservoir data.

\* *p* < 0.05.

According to the three conditions listed in Section 2.3.3, 'WT × COD' and 'TP × AN' were selected as interaction terms for the Shimen Reservoir regression model; 'WT × pH' and 'DO × TP' were the interaction terms for the Mingde Reservoir regression model. Since there was no significant interaction term for the Fongshan Reservoir, we used the result of the time-lag analysis as the final analysis result. The coefficient and correlation of WT, TP, COD, and AN in the Shimen and Mingde Reservoirs should be considered using interaction terms.


**Table 6.** Results from regression analyses of the Mingde Reservoir data.

\* *p* < 0.05.

**Table 7.** Results from regression analyses of the Fongshan Reservoir data.


\* *p* < 0.05.

For example, when using the interaction term 'TP × AN' for the Shimen Reservoir in Equation (6), the concentration of Chl-a is designated as the dependent variable *Y*, the concentration of TP is designated as the independent variable *X*1, and the concentration of AN is designated as the independent variable *X*2. *β*1, *β*2, and *β*<sup>3</sup> are the coefficients of TP, AN, and the 'TP × AN' interaction term, respectively.

Assuming the value of TP is 0.02 mg/L and AN is 0.03 mg/L, then 'TP × AN' is 0.0006 after multiplying the two. Inserting the above values and the coefficients '−26.874<sup>0</sup> , '−54.684<sup>0</sup> , and '3817.1980 into Equation (6) results in a concentration of 0.112 (mg/L) Chl-a when other variables are fixed and the concentration of TP and AN are 0.02 and 0.03 mg/L (Equation (13)).

$$[(-26.874 \times 0.02)] + [(-54.684) \times 0.03] + (3817.198 \times 0.0006) = 0.112\tag{13}$$

When TP and AN both increase by 1 unit, the concentration of TP becomes 1.02 mg/L and AN becomes 1.03 mg/L, and, thus, the value of 'TP × AN' will become 1.0506. Then, the concentration of Chl-a is calculated as 3926.5 mg/L. This means that the unit value of Chl-a increases when other variables are fixed and TP and AN are both increased by 1 unit (Equation (14)).

$$(-26.874) + (-54.684) + 3817.198 + 3817.198 \times (0.02 + 0.03) = 3926.5\tag{14}$$

These results show that when evaluating the water quality trophic state, time-lag and the additive relationship of interactions between factors should also be taken into account to evaluate more accurately the potential for water eutrophication.

Note that the final model for each reservoir is different. This reflects the fact that the eutrophication process in reservoirs is a complex of many different factors, and the process is highly case-by-case. As shown in the introduction, influencing factors sometimes work in opposite directions in different cases from different studies. The result of this study further confirms that even reservoirs in Taiwan show very different patterns when it comes to the factors influencing the trophic state.

### *3.4. Standardization Coefficient*

The regression coefficient of each factor was multiplied by the standard deviation of its data to form a 'standardized coefficient'. The basic amounts and units of each factor were different, and it is difficult to predict the dependent variable values based on only the regression coefficient value of each factor. The relative importance of each factor can be compared by standardizing the coefficients, making it is easier to intuitively understand the degree of influence of each independent variable on the dependent variable.

Tables 8–10 show the standardization coefficient of the Shimen, Mingde, and Fongshan Reservoirs. These results show the degree of influence of each factor on Chl-a while taking into account time-lag and variable interactions.


**Table 8.** Standardization coefficients of variables measured at the Shimen Reservoir.

\* *p* < 0.05. \*\* 'Std. Coef.': Is the factor regression coefficient multiplied by the standard deviation. All other factors are fixed, the degree of influence of each increase by one standard deviation of the factor affects the concentration of Chl-a.

At the Shimen Reservoir (Table 8), the influence of interactions must be considered when analyzing WT, TP, COD, and AN, so these variables will not be discussed separately. Other highly influential factors are conductivity (standardization coefficient = −1.024), secondly pH (standardization coefficient = −0.630), and lastly SS with a lag of one month (standardization coefficient = 0.209).

At the Mingde Reservoir (Table 9), the influence of interactions must be considered when analyzing WT, DO, TP, and pH, so these variables will not be discussed separately. Other factors with a high influence were inflow (standardization coefficient = 11.295), followed by rainfall (standardization coefficient = −11.166), and lastly SD with two months lag (standardization coefficient = 0.548).


**Table 9.** Standardization coefficients of variables measured at the Mingde Reservoir.

\* *p* < 0.05. \*\* 'Std. Coef.': Is the factor regression coefficient multiplied by the standard deviation. All other factors are fixed, the degree of influence of each increase by one standard deviation of the factor affects the concentration of Chl-a.

**Table 10.** Standardization coefficients of variables measured at the Fongshan Reservoir.


\* *p* < 0.05. \*\* 'Std. Coef.': Is the factor regression coefficient multiplied by the standard deviation. All other factors are fixed, the degree of influence of each increase by one standard deviation of the factor affects the concentration of Chl-a.

Data from the Fongshan Reservoir (Table 10) showed that COD had the greatest influence on Chl-a (standardization coefficient = 21.247), which was followed by conductivity (standardization coefficient = −20.725) and lastly SS (standardization coefficient = −3.651).

To summarize, at the Shimen Reservoir, Chl-a was significantly and immediately affected by WT, pH, Conductivity, COD, and AN, and significantly, but not immediately, affected by SS and rainfall. DO, SD, TP, and inflow did not significantly affect Chl-a at the Shimen Reservoir. Chl-a at the Mingde Reservoir was significantly and immediately affected by WT, DO, TP, pH, COD, AN, rainfall, and inflow, while the effect of SD was significant but not immediate. Conductivity and SS did not significantly affect Chl-a at the Mingde Reservoir. WT, SD, Conductivity, SS, and COD significantly and immediately affected Chl-a at the Fongshan Reservoir; TP and pH also significantly affected Chl-a, but not immediately; while the effects of DO, AN, rainfall, and inflow were not significant.

### *3.5. Correlation of Factors*

Table 11 illustrates the correlation results of each weather and water quality factor. Table 11 clearly indicates that the correlation of WT in the Shimen and Mingde Reservoirs needs to be considered within an interaction and is significantly positive in the Fongshan Reservoir. The correlation of DO in the Mingde Reservoir needs to be considered within an interaction, and it is not significantly correlated in the Shimen and Fongshan Reservoirs. The correlation of SD is significantly positive in the Mingde Reservoir and significantly negative in the Fongshan Reservoir, but it is not significantly correlated in the Shimen Reservoir. The correlation of TP in the Shimen and Mingde Reservoirs needs to be considered within an interaction and is significantly positive in the Fongshan Reservoir. The

correlation of pH in the Mingde Reservoir needs to be considered within an interaction and is significantly negative in both the Shimen and Fongshan Reservoirs. The correlations of conductivity in the Shimen and Fongshan Reservoirs are significantly negative but not significant in the Mingde Reservoir. The correlation of SS is significantly positive in the Shimen Reservoir and significantly negative in the Fongshan Reservoir, but there is no significant correlation in the Mingde Reservoir. The correlation of COD in the Shimen Reservoir needs to be considered within an interaction, but it is significantly positive in both the Mingde and Fongshan Reservoirs. The correlation of AN in the Shimen Reservoir needs to be considered within an interaction, and it is significantly positive in the Mingde Reservoir but not significantly correlated in the Fongshan Reservoir. The correlation of rainfall is significantly positive in the Shimen Reservoir and significantly negative in the Mingde Reservoir, but there is no significant correlation in the Fongshan Reservoir. The correlation of inflow is significantly positive in the Mingde Reservoir, but there is no significant correlation in the Shimen and Fongshan Reservoirs.

**Table 11.** Correlation of factors from the Shimen, Mingde, and Fongshan Reservoirs.


d I c: The correlation of factors needs to take into account interaction relationships. d X c: Factors were not significantly correlated. d + c: Factors have a significant positive correlation. d − c: Factors have a significant negative correlation. d N c: The variable is not relevant for the reservoir.

There was no deterministic causality between climate and water quality variables. For example, the pH in the Fongshan Reservoir is negatively correlated with Chl-a, but Zang (2011) shows that Chl-a is positively correlated with pH and DO [29]. The same case as Blumberg (1990) finds that WT is negatively correlated with DO [30], but Chen (2007) shows that WT is positively correlated with DO [31]. In another, case Watson (2016) shows that Chl-a is negatively correlated with DO [32], but Zang (2011) shows a positive correlation [29].

The interaction combinations for the Shimen Reservoir are 'WT × COD' and 'TP × AN', and for the Mingde Reservoir, they are 'WT × pH' and 'DO × TP'. There were no significant correlation interactions for the Fongshan Reservoir. Note that temperature and total phosphorus are the only two factors that have a positive influence among all three reservoirs. However, the effect of temperature negatively interacts with COD in Shimen and with pH in Mingde, making the effect of temperature on the trophic state actually more minor than expected. On the other hand, the interaction terms related to the total phosphorus in Shimen and Mingde are magnifying the effect. This result further supports that total phosphorus is the main factor for the trophic state.

To summarize, these results indicate that the influencing factors of the trophic state in reservoirs defer from case to case; thus, it is difficult to find a one-size-fits-all equation to be perfectly suitable in all cases.

### **4. Conclusions**

The main factor influencing the three reservoirs is total phosphorus. At the Shimen and Mingde Reservoirs, in particular, the interactive effect of TP with other factors on the water quality trophic state was greater than that of TP alone, indicating that more attention should be paid to the interaction effect between the influencing factors. However, there is no significant interaction effect found to further aggravate the trophic state between weather and water quality factors. In the case of these three reservoirs in Taiwan, an additional deterioration of eutrophication from the climate-change-related interaction effect is not a concern.

The analysis of characteristics influenced by time lags and the analysis of the interactions between factors provide a deeper understanding of the correlation between each factor and the degree to which they influence the water quality trophic state. Furthermore, the length of the time lag and the significant combinations of influencing factors vary from reservoir to reservoir, indicating that the patterns of eutrophication might differ according to different reservoir conditions. These results imply that factors influencing the tropic state in a reservoir might vary by reservoir type, geological and meteorological conditions, as well as other potential factors. In other words, forming a model that describes the tropic state for a reservoir is highly case sensitive. The perfect solution of a one-size-fits-all model might not exist. Researchers should carefully review all possible factors before finalizing a model.

In this study, the R<sup>2</sup> values of the MLR model developed for the three reservoirs were all above 0.5, indicating that the regression model for each reservoir explains more than half of the cause of the water quality trophic state. The results indicate that the regression model developed during this study and the methods used are both feasible for assessing the water quality trophic state.

**Author Contributions:** Conceptualization, L.-H.C.H.; methodology, C.-W.H. and L.-H.C.H.; software, C.-W.H.; validation, L.-H.C.H.; writing—original draft preparation, C.-W.H.; writing—review and editing, L.-H.C.H.; visualization, C.-W.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** Ministry of Science and Technology, Taiwan: 109-2221-E-033-004-MY2.

**Data Availability Statement:** Environmental Protection Administration, Environmental Water Quality Information (https://wq.epa.gov.tw/EWQP/en/Default.aspx, accessed on 31 August 2021); Water Resources Agency, Disaster Prevention Information Service (https://fhy.wra.gov.tw/fhy/ Monitor/Reservoir, accessed on 31 August 2021).

**Acknowledgments:** The authors would like to thank Uni-edit (www.uni-edit.net, 31 August 2021) for editing and proofreading this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

