Integral Application of Chemical Mass Balance and Watershed Model to Estimate Point and Nonpoint Source Pollutant Loads in Data-Scarce Little Akaki River, Ethiopia

: The quality of Little Akaki River in Addis Ababa (Ethiopia) is deteriorating signiﬁcantly due to uncontrolled waste released from point and di ﬀ use sources. In this study, pollution load from these sources was quantiﬁed by integrating chemical mass balance analysis (CMB) and the watershed model of pollution load (PLOAD) for chemical oxygen demand, biochemical oxygen demand, total dissolved solid, total nitrogen, nitrate, and phosphate. Water samples monitored bimonthly at 15 main channel monitoring stations and 11-point sources were used for estimation of pollutant load using FLUX32 software in which the ﬂow from the soil and water assessment tool (SWAT) model calibration, measured instantaneous ﬂow, and constituent concentration were used as input. The SWAT simulated the ﬂow quite well with a coe ﬃ cient of determination (R 2 ) of 0.78 and 0.82 and Nash-Sutcli ﬀ (NSE) of 0.76 and 0.80 during calibration and validation, respectively. The uncharacterized nonpoint source load calculated by integrating CMB and PLOAD showed that the contribution of nonpoint source prevails at the middle and downstream segments of the river. Maximum chemical oxygen demand (COD) load from uncharacterized nonpoint sources was calculated at the monitoring station located below the conﬂuence of two rivers (near German Square). On the other hand, high organic pollution load, biochemical oxygen demand (BOD) load, was calculated at a station upstream of Aba Samuel Lake, whereas annual maximum total dissolved solid (TDS), total nitrogen (TN), and phosphate load (PO 4 -P) from the nonpoint source in Little Akaki River (LAR) were found at a river section near Kality Bridge and maximum NO X load was calculated at station near German Square. The integration of the CMB and PLOAD model in this study revealed that the use of area-speciﬁc pollutant export coe ﬃ cients would give relatively accurate results than the use of mean and median ECf values of each land use.


Introduction
Nowadays, urban rivers of developing countries are heavily polluted due to the release of pollutants from the point and nonpoint sources where the determination of accurate pollution load to a river is often difficult due to combined factors of financial, data quality and availability, and technical capability making the river water quality management more challenging [1]. In one way, not only is the determination of waste load from various sources, such as diffuse sources, difficult to quantify as a result of complexity in the generation and uneven distribution of wastewater from various sources [2],

The Study Area
Addis Ababa is the sprawling city and capital, economic, and political center of Ethiopia, found on the border of the greater rift valley at the foothill of the Entoto mountain, with a total land area of 520 square kilometers and population of more than 3 million [19]. The study area (Figure 1) is characterized by a subtropical highland climate with a mean annual maximum and minimum temperature of 24 °C and 12 °C, respectively, and a mean monthly rainfall of 260 mm [20]. The main surface water sources of the city consist of three rivers: Kebena River, Big Akaki River (BAR), and Little Akaki River (LAR), all originating at the foothills of Entoto and draining down to Aba Samuel Lake. The LAR flows from the northwest of the city to the most southern Addis Ababa, before joining the BAR at Aba Samuel Lake, and has a total length of 43 km. The LAR consists of several highly polluted tributaries, mainly located in the middle of the catchment where untreated household waste, including raw sewage, and industrial waste that increases pollution load in the river are discharged. Nearly 65% of the country's industries, ranging from small-scale to large-scale, are concentrated in and around Addis Ababa [21], including food and beverage, textiles, tanneries, rubber, and paper products. The location of most large-scale industries within the vicinity of the LAR that are releasing their wastewater directly to the river without prior treatment (more than 90%) has augmented the pollution load, making the river's water quality unmanageable easily [22]. The river is serving as a natural sewer line for waste originating from various sources, such as domestic, industrial, institutional, and residential areas [23]. The study area is characterized by trachytes, rhyolites, basalts, and several episodes of pyroclastic materials of older volcanic rocks, specifically prevailing in the upper catchments, whereas the western, southwestern, and eastern parts of Addis Ababa are characterized by younger volcanic of trachy-basalt.

Water Quality Sampling and Flow Simulation
A bimonthly water sample was collected for selected constituents in LAR from April 2018 to March 2019 during dry and wet seasons. A total of 11 physico-chemical parameters were collected and analyzed where dissolved oxygen (DO), water temperature, pH, total dissolved solids (TDS), and conductivity were determined onsite. For the estimation of pollutant load in LAR, water samples were collected from 15 main channel monitoring stations (Table 1) Table 2). The samples were collected using a 1.5 L polyethylene bottle, placed in a cooler box, kept under 4 • C, and immediately transported to the laboratory for analysis. Samples of nitrate and phosphate were prefiltered at site and kept in a cooler box before analysis. The analytical methods and instruments used for analysis are shown in Table 3. All the analytical methods were done according to the standard methods for the examination of water and wastewater [24].  On the other hand, instantaneous flow in LAR was measured at the time of sample collection using the current meter (Dentan CM-1AX, Tokyo, Japan), and the gauge data for soil and water assessment tool (SWAT) model calibration were collected from the Ministry of Water, Irrigation, and Electricity. The SWAT model calibrated on Big Akaki River (BAR) outlet was used to simulate and generate flow at each of the subcatchment outlets (monitoring stations) and later used for pollutant load calculation at monitoring stations, along with the instantaneous flow and constituent concentration. The flow generated by SWAT for each subcatchment outlet was later input to the FLUX32 software for load estimation at each monitoring station.

CMB Analysis and Uncharacterized Nonpoint Source Load
In urban rivers of developing countries, the nonpoint source load estimation is a bit challenging due to uncontrolled and irregular waste release rates and unknown and uneven distribution of diffuse sources entry points [25][26][27][28]. To overcome such complexities, various software nowadays are developed to estimate the pollutant loads in a river by taking advantage of simple mathematical equations, from where the diffuse source loads are estimated. FLUX32 is one such software system, developed by the Minnesota Pollution Control Agency (MPCA) to estimate the pollutant loads carried by tributaries and streams. The software requires two data sets: event-based pollutant concentration and respective instantaneous flow and historical gauge recording or model output of the river flow for the specified period [1]. The software uses six different methods to calculate the pollutant load/flux and the choice of each method depends on the sampling approach and variability of flow and concentration [29]. Accordingly, method six (regression applied to individual daily flows) was selected for load calculation in LAR (Equation (1)).
where Q i = mean flow on day i (m 3 /s), ci = measured constituent concentration (mg/L), a = intercept of ln(c) vs. ln(q) regression, b = slope of ln(c) vs. ln(q) regression, SE 2 = standard error of estimate for ln(c) vs. ln(q) regression and q is instantaneous flow (m 3 /s), W i = pollutant load/flux (kg/yr). In LAR, point source load, such as from industries, was calculated by the product of the average discharge rate of wastewater effluents and the mean concentration, where similar approach was used by Amaya et.al. [30], whereas the load from tributaries was calculated by using FLUX32. This is partly because the point source load is often assumed stable and insignificant change occurs seasonally [31,32]. On the other hand, estimation of nonpoint source loads from available monitored water quality data is quite complex and in LAR it was calculated by using upstream-downstream CMB analysis integrating with the watershed model, PLOAD. Since it was difficult to explicitly estimate the nonpoint source load in LAR directly, the term uncharacterized nonpoint source load was used instead, which might include unidentified point source and unrecognized nonpoint source load, where a similar description was used by Jain et.al. [17]. A simple upstream-downstream mass balance approach could be used as an initial estimation of pollution load from lateral sources [33]. Accordingly; where Q D = river flow at the downstream station, C D = downstream constituent concentration, Q Ui = flow of a river at a river section upstream, C Ui = upstream constituent concentration, Σ Losses = the sum of all losses in the stream, L i = is the net load. The above simple approach was used in LAR lateral diffuse pollutant load (differential load) estimation for two basic reasons. First, the span length between the monitoring stations is very small, indicating that the loss is minimum and hence neglected. Second, the river constituents are assumed completely mixed. In the above simple mathematical mass balance equation (Equation (2)), the term ΣLi does not mean it is only contributed from nonpoint source loads and is not the exact net load at a point, but the combination of all loads and the losses and/or generations [34], which could be due to settlement, resuspension, and decay and the generation due to reaction.

Watershed Model Selection
The study of pollutant loads contributing to the pollution of a river is vital for better water quality management. Sekhar et.al. [34] suggested that the catchment pollution management plan should follow a complete study of three components: point sources, nonpoint (background) sources, and natural processes. Though the estimation of point source load is relatively easy, it is challenging to quantify diffuse source loads, specifically in developing countries, where the estimation is often based on simple Sustainability 2020, 12, 7084 6 of 18 empirical equations with limited hydro-meteorological and water quality data. To fill such gaps, many watershed models have been developed and studies were conducted to determine pollutant loads from diffuse sources at catchment scale such as the hydrological simulation program-FORTRAN (HSPF) [35], agricultural nonpoint source pollution model (AGNPS) [36], pollution load (PLOAD) [37], soil and water assessment tool (SWAT) [38], and storm water management model (SWMM) [25]. However, most of the models developed so far are complex and require a large number of data, and hence are not feasible for data-scarce areas like Ethiopia. However, watershed level nonpoint source pollution management could be achieved by using a simple but reliable and relatively accurate model with a reasonable and limited budget. In LAR, hence, we used the PLOAD model due to its versatility, simple data usage, and ease of applicability for the study area [1], integrating with the CMB analysis approach based on monitored water quality data. Most researchers prefer the use of PLOAD due to its cheaper and faster water pollution management of water bodies [2] and the capability and adaptability of the model in different watersheds [39]. On the other hand, Zinabu et.al. [1] recommended the use of the PLOAD model in Ethiopia for nonpoint source pollution management.
The PLOAD is a BASINS (better assessment science integrating point and nonpoint sources) model plugin used to estimate nonpoint source load at catchment level interpreted as an annual load [40]. The model integrates point source and GIS-based land-use data to estimate the nonpoint sources' load contribution from each land use using two approaches: the export coefficient and simple method. Both approaches can be applied based on the data availability and applicability on a watershed, but generally the simple method is used in smaller watersheds, usually less than 1 square mile, while the export coefficient method is used in a mixed land uses [41] for the estimation of constituents such as total suspended solid (TSS), TDS, BOD, COD, NO x (nitrate + nitrite), total Kjeldahl Nitrogen (TKN), ammonia, feacal coliforms (FC), lead, and zinc [40]. In LAR, we used the export coefficient approach where pollutant loads in PLOAD are calculated by where L p = pollutant load (kg/yr), L PU = pollutant export coefficient for each land use (kg/ha/yr), A U = area by certain land use, ha.

Pollutant Export Coefficient
The export coefficient (ECf) is the total amount of pollutant load transported from certain land use per unit area over a specified period of time [42]. When estimating catchment nonpoint source contribution by ECf, each land use is assumed to contribute to the pollutant load per land area and is hence expressed in kg/ha/yr. The watershed shape file was delineated by using ArcSWAT, where the land use in the study area is mostly dominated by urban and agricultural set-ups where informal settlements prevail ( Figure 2). Accordingly, urban land use has the highest percentage coverage with 51.8%, followed by agricultural land (25.72%), forest (10.18%), rangeland (7.2%), bare land (4.63%), and water (0.46%).
Mathematically, the pollutant load using export coefficient with an inclusion of precipitation induced pollution can be expressed by where L i,j is calculated load of constituent i at the subcatchment outlet j (kg/yr); n is the number of land uses contributing; E k,i is the export coefficient of land use k for the constituent i (kg/ha/yr); A k,j is the area of land use k for the subcatchment j; P i,j is precipitation-induced constituent i load at a subcatchment j (kg/yr). P i,j is assumed negligible in LAR. Mathematically, the pollutant load using export coefficient with an inclusion of precipitation induced pollution can be expressed by Li,j= Ek,i×Ak,j+Pi,j n k=1 (4) where Li,j is calculated load of constituent i at the subcatchment outlet j (kg/yr); n is the number of land uses contributing; Ek,i is the export coefficient of land use k for the constituent i (kg/ha/yr); Ak,j is the area of land use k for the subcatchment j; Pi,j is precipitation-induced constituent i load at a subcatchment j (kg/yr). Pi,j is assumed negligible in LAR. Availability of the local pollutant export coefficient is a prerequisite for accurate determination of pollutants loads in a watershed. However, the study area did not have an established pollutant ECf and determination for LAR depends on the data on another watershed elsewhere with nearly similar hydrological, topographical, land use, and climatic set-up. To account for catchment variability and to select appropriate ECf values, we evaluated coefficients globally (Table 4). For pollutant ECf in PLOAD, we reviewed literature values from Ethiopia [1], Canada [43], USA [41], China [26], Taiwan [44], New Zealand [45], Philippines [30], Egypt [46], Lithuania [47], and Japan [11]. Table 4 summarizes the ECf for pollutants and nutrients in literature and selected for LAR for preliminary estimation.  Availability of the local pollutant export coefficient is a prerequisite for accurate determination of pollutants loads in a watershed. However, the study area did not have an established pollutant ECf and determination for LAR depends on the data on another watershed elsewhere with nearly similar hydrological, topographical, land use, and climatic set-up. To account for catchment variability and to select appropriate ECf values, we evaluated coefficients globally (Table 4). For pollutant ECf in PLOAD, we reviewed literature values from Ethiopia [1], Canada [43], USA [41], China [26], Taiwan [44], New Zealand [45], Philippines [30], Egypt [46], Lithuania [47], and Japan [11]. Table 4 summarizes the ECf for pollutants and nutrients in literature and selected for LAR for preliminary estimation.

Calibration and Validation of PLOAD
PLOAD uses GIS-based input data such as land use, watershed boundary, pollutant loading rate (ECf), rainfall depth and optional best management practices (BMPs), terrain imperviousness, and point sources load based on the type of the approach used for estimation. In PLOAD, nonpoint sources from each land use were calculated based on the ECf for each land use. In LAR, once the ECf for each land use was assigned, the PLOAD was calibrated by using the uncharacterized nonpoint source load calculated by CMB analysis and validated using another set of data using an optimized ECf and measured nonpoint load through CMB. The PLOAD model performance was then evaluated by comparing the measured pollutant load (CMB analysis) with the model output until the total percentage error between the measured and model-predicted value became zero or close to zero. Since the model has no direct calibration option, an Excel (2016)-based optimization on Excel Solver was used. In the Excel Solver, the objective to be optimized was set to minimize the percentage total relative error with a possibility of zero value. We selected a GRG nonlinear optimization in Solver due to its faster performance, which uses the local optimum solution. Accordingly, the performance of the model was checked by; where ES is an error of estimation, MPL is measured pollutant load, PPL is PLOAD predicted load.

Results and Discussion
The results in this section are presented in a way that the pollutant load for selected segments of the LAR and monitoring stations of the catchment outlets are represented. The discussion mainly focuses on the major pollution hotspots in the watershed and the pollutant contribution of various land uses were quantified.

Point Sources Load in LAR
The pollutant load from point sources in LAR was much smaller than the tributaries load due to the relatively higher flow rate and pollution level of the tributaries than the point sources. However, stations M3 to M11 (Figure 1) were heavily loaded by point source pollution that contributes significant pollutant loads to LAR including, the soft drink industry, wine industry, abattoir, tobacco factory, and hospitals. Similarly, the heavily polluted Mesalemya stream that joins the main river upstream of outlet M3 and a tributary that crosses densely populated urban center, Merkato, and receives many wastewaters from industries joining the main river at M4 have highly augmented the pollutant load in LAR. Besides, the load contributed by the Addis Ababa Abattoir near Kera Beg Tera (M5) was found to be very high due to significantly higher water consumption from the slaughterhouse that generates wastewater with a higher flow rate. Table 5 summarizes the load contributed by point sources to the LAR. Almost all of the point sources near LAR discharging the wastewater directly to the river have either no treatment plant or couldn't fully operate. From the point source load summary on Table 5, it can be apparently seen that the contribution of stations T6, AA_Kerawoch, AA_WWTPE, and M8 were quite significant for the LAR organic and nutrient pollution.

Flow Simulation and Pollutants Flux in LAR
The SWAT calibrated at BAR was used to generate flow in LAR subcatchment outlets where the model output along with the instantaneous flow and constituent concentration was used in FLUX32 software to calculate the pollutant flux (load). Accordingly, the SWAT simulated the flow quite well in Figures 3 and 4 with an R 2 , NSE, and RSR value of 0.78, 0.76, and 0.49 during calibration and 0.82, 0.8, and 0.45 during validation, respectively. The model performance indicators above (R 2 , NSE and RSR) were good enough to interpret the model output for any purpose. From Figure 3, it can be seen that the deviation between the SWAT model simulated and observed peaks. This could primarily be due to the model performance. The hydrological model performance determined by the Nash-Sutcliff (NSE) was found to be 0.76 and 0.78 during calibration and validation in the study area, which is good enough to interpret the model output [50] and the deviation between the observed and simulated flow is interpreted by the error. Similar results were reported elsewhere with similar trends between the model simulation and observation in the works of Abbaspour et.al. [51], Rostamian et.al. [52], and Shawul et.al. [53]. Sometimes the time lag between the small rainfall event and the main rain event could dictate the variation between the model simulation and observed values. This explanation is supported by the study conducted by Li et.al. [35], who used hydrological and water quality model, HSPF, where a similar trend with this study between the model simulated and observed flow was observed. The spatial location of rain meters and the heterogeneity among rainfall stations could also determine the deviation. Though there was deviation between model simulated and observed values at some peak points, the model performance indicators (specifically the Nash-Sutcliff) could suggest that the model output can be interpreted with a good accuracy. The SWAT-generated subcatchment outlet flow was used to calculate the load in FLUX32. Accordingly, the flow-weighted concentration calculated by method 6 (Equation (1)) was < ±20% of all other methods in FLUX32. The residual plot of bias (as slope) for flow, date, and month at each catchment outlet in LAR was in the range of 0-0.05, which is quite acceptable. Similarly, the plot of slope significance was in the range of 0.88-0.99 ≈ 1. The coefficient of variation (CV) is recommended to be in the range of 0-0.2 during flow-weighted load calculation and in LAR, the CV has resulted in the range of 0.03-0.101, which is quite good. During pollutant load calculation in LAR using FLUX32, the presence of the outlier was checked statistically by testing the significance level, p ≤ 0.05. Sustainability 2020, 12, x FOR PEER REVIEW 10 of 18

Chemical Mass Balance and Pollutants Differential Load in LAR
Once the loads were calculated at each of the subcatchment outlets by FLUX32, a CMB analysis was performed following an upstream-downstream mass balance approach (Equation (2)) to determine the differential load. A similar approach was also followed by [33,34,54]. From the calculated CMB analysis at each monitoring station (Table 6), we found that the prevalence of differential uncharacterized nonpoint source load at the middle and downstream segment of the LAR with highest calculated differential load were BOD, COD, and TDS, which had shown highest loads at segment outlets M9, M10, M12, M13, and M14, and the influence of uncharacterized nonpoint sources were found to be significantly high ( Figure 5). The areas were dominated by urban set-ups (industrial with both large and small scale), residential settlement (usually informal), and agricultural (predominantly small scale) land uses. The organic pollution contribution from the nonpoint sources prevailed in the study area where the maximum calculated differential BOD load was observed at M9 with 1833.2 t/yr, where the station is located downstream of two highly polluted but large streams that make a confluence. On the other hand, the maximum COD differential load from the nonpoint source was observed at the same station with a calculated load of 5588.7 t/yr contributed by the subcatchment having an area of 826.419 ha, followed by M13, M12, and M10, which had contributed

Chemical Mass Balance and Pollutants Differential Load in LAR
Once the loads were calculated at each of the subcatchment outlets by FLUX32, a CMB analysis was performed following an upstream-downstream mass balance approach (Equation (2)) to determine the differential load. A similar approach was also followed by [33,34,54]. From the calculated CMB analysis at each monitoring station (Table 6), we found that the prevalence of differential uncharacterized nonpoint source load at the middle and downstream segment of the LAR with highest calculated differential load were BOD, COD, and TDS, which had shown highest loads at segment outlets M9, M10, M12, M13, and M14, and the influence of uncharacterized nonpoint sources were found to be significantly high ( Figure 5). The areas were dominated by urban set-ups (industrial with both large and small scale), residential settlement (usually informal), and agricultural (predominantly small scale) land uses. The organic pollution contribution from the nonpoint sources prevailed in the study area where the maximum calculated differential BOD load was observed at M9 with 1833.2 t/yr, where the station is located downstream of two highly polluted but large streams that make a confluence. On the other hand, the maximum COD differential load from the nonpoint source was observed at the same station with a calculated load of 5588.7 t/yr contributed by the subcatchment having an area of 826.419 ha, followed by M13, M12, and M10, which had contributed

Chemical Mass Balance and Pollutants Differential Load in LAR
Once the loads were calculated at each of the subcatchment outlets by FLUX32, a CMB analysis was performed following an upstream-downstream mass balance approach (Equation (2)) to determine the differential load. A similar approach was also followed by [33,34,54]. From the calculated CMB analysis at each monitoring station (Table 6), we found that the prevalence of differential uncharacterized nonpoint source load at the middle and downstream segment of the LAR with highest calculated differential load were BOD, COD, and TDS, which had shown highest loads at segment outlets M9, M10, M12, M13, and M14, and the influence of uncharacterized nonpoint sources were found to be significantly high ( Figure 5). The areas were dominated by urban set-ups (industrial with both large and small scale), residential settlement (usually informal), and agricultural (predominantly small scale) land uses. The organic pollution contribution from the nonpoint sources prevailed in the study area where the maximum calculated differential BOD load was observed at M9 with 1833.2 t/yr, where the station is located downstream of two highly polluted but large streams that make a confluence. On the other hand, the maximum COD differential load from the nonpoint source was observed at the same station with a calculated load of 5588.7 t/yr contributed by the subcatchment having an area of 826.419 ha, followed by M13, M12, and M10, which had contributed an annual load of 3217.2 t, 2306.9 t, and 1859.55 t, respectively. Similarly, high BOD depletion was observed at downstream stations M10 (1867.49 t/yr), M11 (738.49 t/yr), and M14 (1007.76 t/yr), where the impacts of nonpoint sources were assumed to be insignificant, showing that the river is going through a high rate of recovery, possibly due to morphometric effects and the improved self-purification capacity of the river [55]. Pa = parameters; differential load = incremental load of the downstream station relative to the upstream station/s; NO x * = reported as NO 3 +NO 2 ; † deficit (sink).
High PO 4 -P differential load was calculated at most downstream stations characterized by small-scale urban agricultural activities prevailing at M12 (23.63 t/yr), followed by M10 (12.92 t/yr) and M9 (12.42 t/yr). Besides, the area is characterized by a large number of small-scale industries, dumping sites for informal solid waste including bio-waste, wastewater treatment plant effluent, and animal remains. On the other hand, from Table 6, it can be seen that high differential TN load from uncharacterized nonpoint source was calculated at monitoring stations M9 (157.2 t/yr), M10 (674.79 t/yr), and M11 (97.02 t/yr), whereas stations M4, M12, and M13 were identified as areas with TN sink. COD, BOD, TDS, and TN were found to be dominant nonpoint source load contributions and were found in large quantities prevailing in the middle and downstream segments, indicating the increased impact of washouts from agricultural and urban land uses. The CMB analysis in LAR revealed that most of the organic waste load were concentrated at the middle segment of the river, whereas the sink for these pollutants was found far downstream. The CMB analysis in LAR also showed that NO x , TN, and BOD had sinks at most of the LAR monitoring stations calculated at outlets. Stations M9-M14 are the recognized sink areas for organic pollutants and nutrients located in the middle and downstream segments of LAR. A similar finding was reported by Elósegui et al. [56], where most of the nutrient load in river Agüera in Northern Spain was retained in the middle of the river. The station M11 is a place where the organic pollutants are highly degraded and hence the area was identified as a major nutrient sink. It is a place where a public protected park (Biheretsige Park) is located, which could reduce the impact of nonpoint source load contributions. Relative to the downstream and middle catchments, the upstream segments have low nonpoint source load contributions, where a similar result was also found by Jain et.al. [33] on Hindon River, India. This could be due to relatively low flow, less anthropogenic influence in the area, and the presence of buffer zones such as grass strips and urban forests.
From the CMB analysis ( Figure 6), it can be seen that most of the constituents in the upstream section of the LAR (M1-M5) had minimum differential load, showing reduced impact of nonpoint source pollution. Conversely, differential loads of BOD, TN, and NO x have shown most of the sinks were observed at the middle and downstream segment of LAR, where riversides are protected by grasses and plants, predominantly at M10 and M11. Maximum loads of nutrients, such as TN, NO x , and PO4-P, were recorded at M10, which was also identified as a major sink for BOD. Similarly, M14 was found as an area where both TN and NO x have shown a positive differential load. The station is found downstream of the discharge point of wastewater treatment plant effluent and is also characterized by small-scale urban agriculture. an annual load of 3217.2 t, 2306.9 t, and 1859.55 t, respectively. Similarly, high BOD depletion was observed at downstream stations M10 (1867.49 t/yr), M11 (738.49 t/yr), and M14 (1007.76 t/yr), where the impacts of nonpoint sources were assumed to be insignificant, showing that the river is going through a high rate of recovery, possibly due to morphometric effects and the improved selfpurification capacity of the river [55].

Integration of CMB-PLOAD and Uncharacterized Nonpoint Source Load
Loads calculated at selected LAR subcatchment outlets (monitoring stations) by CMB were used for PLOAD calibration and the selected ECf for LAR was used as an initial estimation for calibration. During the calibration of PLOAD in Excel Solver, the ECf were used as independent variables and the ranges were constraints by setting the upper and lower bounds of the ECf based on the literature for optimization. At the initial stage of calibration and preoptimization, the total percentage error between the model-predicted and measured load at all monitoring stations for COD, BOD, TDS, NOX,

Integration of CMB-PLOAD and Uncharacterized Nonpoint Source Load
Loads calculated at selected LAR subcatchment outlets (monitoring stations) by CMB were used for PLOAD calibration and the selected ECf for LAR was used as an initial estimation for calibration. During the calibration of PLOAD in Excel Solver, the ECf were used as independent variables and the ranges were constraints by setting the upper and lower bounds of the ECf based on the literature for optimization. At the initial stage of calibration and preoptimization, the total percentage error between the model-predicted and measured load at all monitoring stations for COD, BOD, TDS, NO X , PO 4 -P, and TN were 3680.6%, 2965.64%, 767.97%, 446.8%, 750.40%, and 899.36% respectively. After the optimization of the ECf in Solver, the average percentage relative error dropped to 12.16%, 3.98%, 0.53%, 26.86%, 8.9%, and 29.22%, respectively. On the other hand, the sum of errors in COD was found to be higher than BOD, which resulted in a total sum of error at 3680.6% preoptimization, where much was attributed by station M5, at which point the PLOAD underestimated both constituents at the station.
The PLOAD prediction for TDS was relatively accurate, where the total relative error at all monitoring stations before optimization was found to be 767.97%, which sharply dropped to 0.53% after optimization. Unlike other parameters, the PLOAD model has overestimated the TDS load in the upstream section of the LAR (M1-M5, and M9), whereas the model underestimated the PO 4 -P load in LAR at the downstream segment of the river. High total errors were recorded at smaller catchments of LAR than the bigger catchments, and Zinabu et.al. [1] also came up with a similar finding in Kombolcha catchment, Ethiopia. From the optimized ECf, we found that urban and agricultural land use ECf varies greatly with varying urban land-use types and has shown a significant variation spatially. The water body land-use showed minimum change over the loading of constituents across the monitoring stations, which could be due to less area coverage in the watershed.
The PLOAD was rerun using the median and mean value of the optimized and calibrated ECf. Though the median ECf gave minimum total percentage error relative to the mean value, high variation in the ranges of the ECf made the load vary greatly, and hence the loads calculated by both mean and median ECf were not acceptable. The optimized ECf values revealed that the use of average and median value resulted in an underestimation of pollutant load at some catchment outlets, whereas overestimation on others. Thus, the use of area-specific ECf was relatively acceptable and was found effective in better estimation of pollutant loads in LAR. The study conducted by Shrestha et al. [11] also recommended the use of area-specific, and development of local, ECf for effective pollutant load calculations. The optimized values can then be interpreted well for LAR and used for the management of the nonpoint source pollution load in the catchment. Similarly, the PLOAD was validated for a different data set without change in the optimized export coefficient and the error calculated from the model was acceptable enough for further interpretation. Accordingly, the percentage error between the PLOAD predicted and measured values for COD, BOD, TDS, NO x , PO 4 -P, and TN during validation were found to be 16.41%, 7.06%, 1.77%, 13.23%, 5.4%, and 18.83%, respectively.
The calibrated pollutant ECfs showed that the urban land uses had significantly varying export coefficients. Accordingly, the pollutant loading rate for urban land use ranged from 42.1-2083 kg/ha/yr, 63.3-2012 kg/ha/yr, and 0.49-1.6 kg/ha/yr for COD, TDS, and PO 4 -P, respectively, varying with the subtype of urban land use (such as residential, commercial, industrial) and location. The urban land uses dominated by residential settlements and industries do have high loading rates, whereas other urban land uses have relatively lesser rates. Similarly, the urban land use pollutant ECf for BOD and NO x also vary greatly, ranging from 120-1950.55 kg/ha/yr and 0.1-47.32 kg/ha/yr respectively. The high variation in ECf spatially is due to the high difference in the impact of nonpoint sources among various urban land uses. The upstream urban land uses have less pollutant loading than the middle and downstream catchments (

Conclusions
In this study, conjunctive application of chemical mass balance and watershed model, PLOAD, was used to estimate the nonpoint source load in the data-scarce Little Akaki River, Ethiopia, which was found effective. The approach proved to be more efficient in the study area, which ultimately focused on the determination of organic pollutants and nutrient loads based on a continuously monitored water quality data in the river. The following major conclusions were drawn from the research findings.

•
The impact of nonpoint sources in the upstream segment of the LAR catchment was relatively less than the downstream and middle segments, primarily due to the reduced impacts of unrecognized point sources, less urban settlements, better land-use protection and management. Moreover, lesser flow rate in the upper segment of the river could be playing a critical role for the lower diffuse source load in the area.

•
The integration of CMB and catchment nonpoint source pollution models such as PLOAD could be an effective and alternative pollutant load estimation approach in data scarce areas.

•
The nonpoint source pollutant load was found to be very high in areas where urban land uses prevail, followed by agricultural and barren land uses, indicating the nonpoint source pollution management focus areas. Mitigation measures involving these land uses is recommended. • Area-specific (local) pollutant export coefficients were found to be more effective and accurate load estimation approach than the use of mean and median export coefficients, which ultimately give a lower error during pollutants load calculation. Despite higher accuracy of CMB for the estimation of differential uncharacterized nonpoint source load estimation, integrating with a simple watershed model was found to be a good alternative for a more accurate representation of the diffuse source load. • Under-and overestimation of the PLOAD for pollutant load estimation was observed at some catchment outlets, which when integrated with CMB analysis gave a promising result. But adaptation of global export coefficient to the local condition with different hydro-climatic setup is the limitation of the integral modeling approach.

•
The pollutant and nutrient export coefficients developed in LAR catchment could be transferred to other catchment elsewhere in the country for similar application for preliminary nonpoint source pollutant load management. The accuracy and effectiveness of the CMB for nonpoint source load estimation highly depends on a number of factors, such as frequency of data collection, distance between the monitoring stations, and identification of the major point sources to the river.
It could be concluded that the integral application of chemical mass balance and watershed models such as PLOAD could be a better option for the estimation of nonpoint source pollutant loads in areas with few monitored water quality data. Future studies incorporating the vast and long-term monitoring program at larger catchment scale would be helpful for better pollution load management in the river.