Modeling Biological Oxygen Demand Load Capacity in a Data-Scarce Basin with Important Anthropogenic Interventions

: Most water bodies are currently used as receptors for pollutants coming mainly from the industrial and domestic sectors. The Biob í o river is subjected to multiple anthropogenic pressures such as industrial water supply, drinking water, hydroelectric power generation, agriculture, and the ﬁnal receptor body of a large amount of industrial and urban waste, pressures that will intensify due to the decrease in water ﬂow as a result of climate change. In this context, organic contamination has been found mainly from sewage discharges and oxidizable waste discharges generated by industrial processes. In this sense, the objective of this research is to determine the Biological Oxygen Demand Loading Capacity (LC) in a basin with a low density of water quality data subjected to strong anthropogenic pressures. To estimate the carrying capacity in a section of the Biob í o River, the water quality model River and Stream Water Quality Model- Qual2K version 2.11b8, developed by Chapra, was used. This model solves the Streeter–Phelps equation, proposing an analytical expression to relate the dissolved oxygen (DO) and biochemical oxygen demand (BOD 5 ) variables. These variables were modeled for different critical scenarios of minimum ﬂows in return periods of 5, 50, and 100 years, determining that the studied section of the Biob í o river would have a high carrying capacity to not be affected by its organic matter pollution.


Introduction
Global freshwater demand has shown an average growth of 1% per year since the 1980s, mainly due to economic and population growth [1]. The sectors with the highest global water consumption are: agriculture, accounting for 69% of annual water withdrawals, the industrial sector with 19%, and domestic consumption with 12%, and this is expected to increase by 20-30% by 2050 [1]. In this sense, the overexploitation of water resources has generated severe environmental problems, e.g., water pollution, loss of aquatic habitat, fragmentation of ecosystems, and others, all associated with the anthropic pressures reducing the quantity and quality of available water [2][3][4][5].
The Latin American and Caribbean region has 31% of the world's freshwater resources [5]. However, water pollution is a problem that is increasing due to continued urbanization and economic development [1]. Although the agricultural sector is the largest Water 2021, 13, 2379 2 of 12 consumer of water, the urban and industrial sectors are the largest emitters of toxic substances to water bodies, with wastewater discharges being the primary source of pollution and loss of water quality [3,5,6].
The quality of water bodies depends on natural processes such as climatological conditions (e.g., precipitation, temperature), soil types, vegetation cover, erosion, runoff, and the loss of water due to anthropogenic pressures; exploitation of water resources, changes in land use (agriculture), industrial activities, and urban settlements [7][8][9].
Additionally, the main effect of climate change in freshwater bodies is to increase the frequency of extreme events, especially droughts, increasing the concentrations of pollutants present in surface water bodies [10,11]. Droughts, reduced flow rates and increased water temperature cause a decrease in dissolved oxygen, altering the dilution and self-purification capacity in rivers expressed as an increase in suspended biochemical oxygen demand [12].
Dissolved oxygen is an indicator of the quality and health of aquatic ecosystems [13][14][15][16][17][18][19][20][21], as well as of the self-purification processes of natural systems [17,18,[22][23][24]. Self-purification is a slow process, which depends on the: (a) flow rate, (b) water movement, and (c) depth and surface area of the water body. In this sense, if the concentration of the pollutant is very high, the water flow is required to travel long distances before achieving acceptable concentrations [22,[24][25][26].
Load capacity (LC) is a concept frequently used in river water quality management, which is defined as the maximum quantity of a particular substance that can exist in a receptor body without altering the ecological characteristics of the system [27,28]. The US EPA, through the Clean Water Act (CWA), proposes the Total Maximum Daily Load (TMDL) as an indicator of LC [29][30][31]. This concept is defined as the maximum daily concentration possible to discharge into a water body to meet specific environmental quality standards at a given location. Additionally, it is possible to quantify this load in other temporal scales such as monthly, seasonal, or annual [29,30,[32][33][34][35].
Over the years, water quality models have been developed as a fundamental tool in the prediction and control of pollutants due to their ability to represent biological, chemical, and physical processes and changes in aquatic ecosystems [32,36]. The most used models are QUAL2K [16,32], MIKE11 [37], HEC-RAS [38,39], and SWAT [40,41], which predict the behavior of river water quality using hydrology and hydraulics variables and pollutants characteristics [42]. However, the application of these models faces significant challenges due to the large amount of detailed and spatialized data required, e.g., land use, time series of pollutant concentration, daily flows series, and samples along the river. Information frequently is not available in developing countries, which are often limited and can have significant uncertainties [42][43][44][45].
In this context, the objective of this research is to determine the Biological Oxygen Demand LC in a basin with a low density of water quality data subjected to strong anthropogenic pressures. The Biobío river basin in Chile was selected for this purpose, due to its low density of information and the multiple anthropogenic pressures to which it is subjected, such as industrial water supply, drinking water, hydroelectric power generation, agriculture, and final receptor body of a large amount of industrial and urban waste.

Study Area Study Basin
The Biobío river basin is located in south-central Chile, specifically between coordinates 36 • 45 -38 • 49 S and 71 • 00 -73 • 20 W (Figure 1). It covers an area of 24,260 km 2 , being the third largest basin in the country. The main channel is 380 km and is the second-longest river in Chile [46]. It is also located between two political-administrative regions: the Biobío Region (72% of the total area) and the Araucanía Region (28%) [47]. The hydrological regime of the Biobío basin is pluvio-nival, with an average flow at its source being 30 m 3 /s [47] and average flow at its mouth being 960 m 3 /s, with a maximum and minimum flow average of 1600 m 3 /s (July) and 160 m 3 /s (March) [48]. In addition, it is subject to the influence of different environments and geographical factors, causing the dynamics of the system to be highly variable from the beginning of its course to its mouth [47][48][49].
Water 2021, 13, x FOR PEER REVIEW 3 of 12 influence of different environments and geographical factors, causing the dynamics of the system to be highly variable from the beginning of its course to its mouth [47][48][49]. This study considers a 159.6 km of the Biobío River between the cities of Rucalhue and Hualqui (Figure 1). The selection of the section considered the existence of different point sources (sewage treatment plants and industrial effluents). The fluviometric and water quality records were obtained from the Chilean Water Authority (DGA) [50] and the Biobío River Monitoring Program (PMBB) of the EULA (European-Latin America) Center-Chile.

Water Quality Legislation in Chile
In 2004, Chile began issuing Secondary Environmental Quality Standards (NSCA) to ensure the protection and conservation of aquatic ecosystems of inland surface waters and reduce pollution and maintain their quality as much as possible. Currently, five Secondary Standards have been issued, two for lake bodies (Lago Llanquihue and Lago Villarrica), and three for river basins (the Serrano river basin, Maipo river basin, and Biobío river basin) [51].
The secondary environmental quality standard for the Biobío river basin aims to conserve or preserve aquatic ecosystems and their ecosystem services by maintaining or improving water quality in the basin. Therefore, environmental quality standards were established for 19 pollutants, including DO and BOD5. Fourteen monitoring areas were implemented to comply with the standards, six along the Biobío River and eight in its main tributaries [51,52]. This study considers a 159.6 km of the Biobío River between the cities of Rucalhue and Hualqui ( Figure 1). The selection of the section considered the existence of different point sources (sewage treatment plants and industrial effluents). The fluviometric and water quality records were obtained from the Chilean Water Authority (DGA) [50] and the Biobío River Monitoring Program (PMBB) of the EULA (European-Latin America) Center-Chile.

Water Quality Legislation in Chile
In 2004, Chile began issuing Secondary Environmental Quality Standards (NSCA) to ensure the protection and conservation of aquatic ecosystems of inland surface waters and reduce pollution and maintain their quality as much as possible. Currently, five Secondary Standards have been issued, two for lake bodies (Lago Llanquihue and Lago Villarrica), and three for river basins (the Serrano river basin, Maipo river basin, and Biobío river basin) [51].
The secondary environmental quality standard for the Biobío river basin aims to conserve or preserve aquatic ecosystems and their ecosystem services by maintaining or improving water quality in the basin. Therefore, environmental quality standards were established for 19 pollutants, including DO and BOD 5 . Fourteen monitoring areas were implemented to comply with the standards, six along the Biobío River and eight in its main tributaries [51,52].

Water Quality Model
The most widely used method for determining the variation of DO, BOD, and loading capacity concentrations in a river course is the Streeter and Phelps (1925) model [53]. The empirical equations proposed by these researchers have been widely used to evaluate the impact of waste discharges on dissolved oxygen concentrations and assimilative capacity in rivers [14,16,[53][54][55][56][57][58][59]. Currently, a large number of water quality models solve the equations developed by Streeter and Phelps, the most widely applied being the River and Stream Water Quality Model-Qual2K [16,17,32,35,60,61]. This one-dimensional mathematical model assumes complete mixing of the flow in both vertical and transverse directions and divides the reach of interest into a specific number of computational elements with homogeneous hydraulic characteristics [62][63][64][65].
In the study section, the 159.6 km river reach was divided into seven sections with similar hydraulic characteristics, i.e., slope, width, and elevation ( Figure 2). For this purpose, ArcGis 10.1 software was used with a 30 × 30-m-resolution digital elevation model (DEM) that was obtained from the Shuttle Radar Topography Mission [66]. Given the complexity of the system, the average velocity and depth of water flow in the study section were determined using the River Analysis System (Hec-Ras) 4.0 model developed by the Hydrologic Engineering Center of the U.S. Army Corps of Engineers [16].
Water 2021, 13, x FOR PEER REVIEW 4 of 12

Water Quality Model
The most widely used method for determining the variation of DO, BOD, and loading capacity concentrations in a river course is the Streeter and Phelps (1925) model [53]. The empirical equations proposed by these researchers have been widely used to evaluate the impact of waste discharges on dissolved oxygen concentrations and assimilative capacity in rivers [14,16,[53][54][55][56][57][58][59]. Currently, a large number of water quality models solve the equations developed by Streeter and Phelps, the most widely applied being the River and Stream Water Quality Model-Qual2K [16,17,32,35,60,61]. This one-dimensional mathematical model assumes complete mixing of the flow in both vertical and transverse directions and divides the reach of interest into a specific number of computational elements with homogeneous hydraulic characteristics [62][63][64][65].
In the study section, the 159.6 km river reach was divided into seven sections with similar hydraulic characteristics, i.e., slope, width, and elevation ( Figure 2). For this purpose, ArcGis 10.1 software was used with a 30 × 30-m-resolution digital elevation model (DEM) that was obtained from the Shuttle Radar Topography Mission [66]. Given the complexity of the system, the average velocity and depth of water flow in the study section were determined using the River Analysis System (Hec-Ras) 4.0 model developed by the Hydrologic Engineering Center of the U.S. Army Corps of Engineers [16].  The calibration and validation of the model were carried out with data from the monthly DO concentrations of the PMBB water quality stations. The calibration process was performed for March 2009, calibrating the parameters k a (Reaeration Coefficient), k d (Deoxygenation Coefficient), and the temperature correction coefficients θ a and θ d [65]. The validation period was for March 2012. The goodness of fit (GoF) indicators Coefficient of Determination (R 2 ), Nash-Sutcliffe Efficiency (NSE), Percent Bias (PBIAS), and the Agreement Index (d) [19,[67][68][69] were used for calibration and validation periods. Finally, The LC of BOD 5 was evaluated for the extreme drought flows with return periods of 5, 50, and 100 years. The fluviometric record of the Biobío River Station in Rucalhue (BNA code 08317001-8) [50] from 1960 to 2013 was used for this analysis [70].
To determine the maximum BOD concentration that can be purified by the Biobío River in extreme minimum flow scenarios, a minimum DO concentration of 5 mg/L was considered. This value is determined in Chilean Standard NCh 1.333 [71], which establishes water quality criteria according to physical, chemical, and biological aspects, depending on the user assigned to the water body. Table 1 shows the GoFs indicators results in the calibration and validation phases. The model performance during calibration can be classified as very good, in the case of R 2 , NSE, and d indicators, and satisfactory in the case of PBIAS, according to Moriasi et al. (2015) [68]. In the validation period, R 2 decreases to good and PBIAS changes to underestimation; according to Moriasi et al. (2015) [68], the other GoF remained in the same category. These results show that the efficiency criteria of the model are high, which indicate that the Qual2K acceptably represents the behavior of the variables studied for the study section of the Biobío River but tends to underestimate the DO concentration. Finally, Table 2 shows the values obtained for the parameters k a , k d , θ a, and θ d for the seven sections used. Values obtained by calibration for the different governing parameters for the system ( Table 2) show that the k a coefficient presents a high range of variation among the seven sections of the studied reach. This is because the surface reaeration process is one of the primary and variable sources of oxygen in the river system and directly affects DO concentrations. Langbein and Durum (1967) [72] and Bowie et al. (1985) [73] indicated that defining a range for this coefficient is very complex since it depends on factors such as wind speed, temperature, and hydraulic parameters. Similarly, Langbein and Durum (1967) [72] estimated that k a variations are mainly attributable to changes in the hydraulic characteristics of rivers, e.g., width, slope, area, and length of the riverbed, which generate variations in depth and flow velocity, causing the river ecosystem to adopt the conditions for DO concentrations to increase or decrease. They also note that in rivers with high depths, the k a value is lower than rivers with shallow depths, concluding that the higher the k a coefficient, the more intense the reaeration process. On the other hand, Link (1998) [74] determined that wastewater discharge into rivers decreases their capacity for reaeration since contaminated water causes the diffusion of DO to be lower. Therefore, the wide range obtained for this coefficient is not surprising since the Biobío River has significant hydraulic differences along its course.

Results and Discussion
The k d coefficient presented lower variations between sections. However, this parameter did not show significant effects within the DO modeling. This behavior is due to the high dilution capacity of the Biobío River, which is why BOD 5 is diluted very quickly when it encounters the water layer, generating a lower oxygen consumption due to BOD degradation. It should also be considered that effluent treatment decreases the oxidizable fraction of BOD, causing its oxidation in the water body to be slower. The temperature correction factor θ a was the same for all sections. This factor presented a high incidence in the DO concentrations, which shows that the DO is sensitive to temperature. The temperature is a factor that affects the reaeration process; increasing the factor decreases the oxygen concentration. As in the previous case, the factor θ d had the same value for all the sections. This factor had a low incidence on dissolved oxygen concentrations since it affects oxygen consumption by BOD degradation, specifically the k d coefficient.
To determine the extreme drought flow scenarios for return periods of 5, 50, and 100 years, a frequency analysis was performed ( Figure 3) for a series of minimum monthly flows available at the Biobío River Station in Rucalhue. The Anderson-Darling goodness of fit test was applied with a 95% confidence interval [75], as shown in Table 3, using as a criterion the smallest difference between the critical value and the A 2 test statistic. The probability distribution that best fit the data series of minimum monthly flows was the Generalized Extreme Value, which gives as results a flow of 82.59 m 3 /s for a 5-year return period, 66.82 m 3 /s for a 50-year return period, and of 64.32 m 3 /s for a 100-year return period. defining a range for this coefficient is very complex since it depends on factors such as wind speed, temperature, and hydraulic parameters. Similarly, Langbein and Durum (1967) [72] estimated that ka variations are mainly attributable to changes in the hydraulic characteristics of rivers, e.g., width, slope, area, and length of the riverbed, which generate variations in depth and flow velocity, causing the river ecosystem to adopt the conditions for DO concentrations to increase or decrease. They also note that in rivers with high depths, the ka value is lower than rivers with shallow depths, concluding that the higher the ka coefficient, the more intense the reaeration process. On the other hand, Link (1998) [74] determined that wastewater discharge into rivers decreases their capacity for reaeration since contaminated water causes the diffusion of DO to be lower. Therefore, the wide range obtained for this coefficient is not surprising since the Biobío River has significant hydraulic differences along its course.
The kd coefficient presented lower variations between sections. However, this parameter did not show significant effects within the DO modeling. This behavior is due to the high dilution capacity of the Biobío River, which is why BOD5 is diluted very quickly when it encounters the water layer, generating a lower oxygen consumption due to BOD degradation. It should also be considered that effluent treatment decreases the oxidizable fraction of BOD, causing its oxidation in the water body to be slower. The temperature correction factor θa was the same for all sections. This factor presented a high incidence in the DO concentrations, which shows that the DO is sensitive to temperature. The temperature is a factor that affects the reaeration process; increasing the factor decreases the oxygen concentration. As in the previous case, the factor θd had the same value for all the sections. This factor had a low incidence on dissolved oxygen concentrations since it affects oxygen consumption by BOD degradation, specifically the kd coefficient.
To determine the extreme drought flow scenarios for return periods of 5, 50, and 100 years, a frequency analysis was performed ( Figure 3) for a series of minimum monthly flows available at the Biobío River Station in Rucalhue. The Anderson-Darling goodness of fit test was applied with a 95% confidence interval [75], as shown in Table 3, using as a criterion the smallest difference between the critical value and the A 2 test statistic. The probability distribution that best fit the data series of minimum monthly flows was the Generalized Extreme Value, which gives as results a flow of 82.59 m 3 /s for a 5-year return period, 66.82 m 3 /s for a 50-year return period, and of 64.32 m 3 /s for a 100-year return period.   Figure 4 shows the simulations made with Qual2K by the three extreme minimum flow scenarios. Figure 4a shows the DO and BOD concentrations by minimum flow with a return period of 5 years (82.59 m 3 /s). The BOD concentration increases downstream, reaching its maximum concentration at 46.3 km, and is equivalent to 1.49 mg/L, while the dissolved oxygen concentration decreases downstream, reaching its minimum concentration at 46.3 km, equal to 8.63 mg/L. Figure  The results show that the highest BOD purification capacity is reached for a minimum flow with a return period of 5 years, equivalent to 82.59 m 3 /s with a concentration of 2078 mg/L. As the return period increases, the BOD concentration that can be purified decreases, reaching the lowest purification capacity in a return period of 100 years, with a minimum flow of 64.32 m 3 /s, filtering a maximum concentration of BOD of 1580 mg/L. For the three simulated scenarios, the minimum DO concentration is reached at 104 km ( Table 4). The high load capacity of the Biobío river can be attributed to its high flows, which means that it has a high DBO dilution capacity without considerably altering the aquatic ecosystem, in agreement with reported by [76]. When comparing the results presented in Figure 4 with the environmental quality levels established in the Secondary Environmental Quality Standard for the Protection of Surface Continental Waters of the Biobío River Basin, in the section from Rucalhue to upstream of the Vergara River (191 to 106.5 km) the DO and BOD concentrations obtained for the return periods of 5, 50, and 100 years comply with the environmental quality levels established by monitoring station BI-30, setting a minimum DO concentration of 9 mg/L and a maximum BOD concentration of 2 mg/L. From upstream of the confluence of the Vergara River to upstream of the confluence with the Gomero River is the BI-40 monitoring station, which establishes a minimum DO concentration of 9 mg/L and a maximum BOD concentration of 2 mg/L. Therefore, from 106.5 to 66.3 km, the DO concentrations obtained for the return periods do not comply with the environmental quality levels, unlike the BOD concentrations, which comply with the levels established by the regulations. Finally, upstream from the Gomero River to the mechanical bridge is the monitoring area of station BI-50, which establishes a minimum DO concentration of 8.7 mg/L and a maximum BOD concentration of 2 mg/L. Therefore, from 66.3 to 54 km, the environmental quality standards for DO and BOD are met; however, from 54 to 32 km, the standard levels for DO concentrations are not met. The combination of these types of tools is useful in designing water quality management plans and developing decontamination plans to understand the local situation and control surface water pollution in multiple-use watersheds.
Water 2021, 13, x FOR PEER REVIEW 8 of 12 Finally, upstream from the Gomero River to the mechanical bridge is the monitoring area of station BI-50, which establishes a minimum DO concentration of 8.7 mg/L and a maximum BOD concentration of 2 mg/L. Therefore, from 66.3 to 54 km, the environmental quality standards for DO and BOD are met; however, from 54 to 32 km, the standard levels for DO concentrations are not met. The combination of these types of tools is useful in designing water quality management plans and developing decontamination plans to understand the local situation and control surface water pollution in multiple-use watersheds. (a)

Conclusions
According to the efficiency criteria used, the predictions made by the Qual2K model are satisfactory, both in the calibration and validation phases, so this model is a good approach by estimated the Biobío River BOD load capacity.
It was found that the most influential parameters in the model response were the k a coefficient and the θ a factor, demonstrating that the reaeration process and temperature determine the DO concentrations in the river.
Based on the different extreme minimum flow scenarios, the model response indicates that the section studied has a high loading capacity for BOD, due principally to its high flows in the summer months. Therefore, under the modeled conditions, the Biobío River does not present organic matter pollution and has a high available capacity for BOD assimilation. According to the values obtained, the highest BOD loading capacity is reached in the 5-year return period for a minimum flow of 82.59 m 3 /s and a BOD concentration of 2078 mg/L, and the lowest loading capacity is obtained in a 100-year return period, equivalent to 1580 mg/L of BOD and minimum flow of 64.32 m 3 /s. This study shows that the Qual2K model coupled to Hec-Ras is an effective tool for calculating the permissible BOD load capacity in the Biobío River section. The application of these models can be extended to similar problems with other contaminants such as nitrogen and phosphorus.
The results of the present investigation demonstrate that the calibration and validation of this type of tool can be of great use in the future in the planning and control of water pollution in the Biobío River, in addition to providing information on the size of BOD concentrations that may have future discharges into this receiving water body.