## 2 Materials and Methods
### Study Framework  
Figure 2 illustrates the methodology adopted in this study. The dynamic problem was specified. Variables and parameters for the river resources system were then identified based on the dynamic problem. In this case, the IWI guided the identification of variables and parameters and tracked the sustainability of the river system. SD methodology was then used to construct the model, where both secondary data and data created by QGIS (Figure 3) were applied for project impact assessment. The model comprised causal loop and stock-flow diagrams (Figures 4 and 5). Vensim PLE (8.2.0) ([[https://vensim.com/](https://vensim.com/)]([https://vensim.com/](https://vensim.com/)), Ventana Systems, Harvard, MA, USA, 1 December 2020) was used to draw the causal loop and stock-flow diagram [46]. The causal diagram provides a comprehensive overview of the resource system. The stock-flow diagram is the main structure for the quantitative calculation of the model. Data was collected from a wide range of sources to develop the model database, which was categorized into four parts: water quality, natural capital, produced capital, and human capital. Data were collected directly from secondary sources and derived from remote sensing (Landsat 8) data. For example, the study employed QGIS to derive vegetation and abandoned land areas from vector and raster data. Administrative boundary vector data was used to mark the boundary of the study setting, obtained from DIVA-GIS website providing data-interpolating variational analysis. Raster data, which included the remote sensing images (Landsat 8), was obtained from the United States Geological Survey Earth Explorer ([[https://earthexplorer.usgs.gov/](https://earthexplorer.usgs.gov/)]([https://earthexplorer.usgs.gov/](https://earthexplorer.usgs.gov/)), Retrieved on 18 February 2021). The semi-automatic classification plugin in QGIS was employed to develop the land use land cover (LULC) map of this study area, after which the area of each class was calculated using raster calculator (unique values report) in QGIS (Figure 3) [47]. Similarly, the data on BOD per person were obtained through the linear interpolation of historical BOD and population data. To obtain this data, we had measurements of BOD for 2010, 2012, 2013, 2015, 2018, and 2019. We estimated the K1 as the outflow rate. We developed a model of BOD over time as a function of the inflow and the outflow. We knew the outflow was the BOD(t) * K1. We needed to estimate the inflow. We estimated the inflow as a smooth curve BODIn(t). We assumed the BODIn(t) was based on the Population(t) [48].  
So, we found the BOD/person(t) that when multiplied by the Population(t) gave us then BODIn(t) that provided the minimum squared error between the measured BOD and the model BOD(t). That gave us values of BOD/person(t) for 2010, 2012, 2013, 2015, 2018 and 2019. Then we used linear interpolation to estimate the values of BOD/person(t) for the year of 2020. Data analysis was conducted by comparing different values of key variables obtained from the model simulation. For example, the three scenarios produced three DO, BOD, and IW values, based on different streamflow, waste control, and migration control values. Finally, sustainability was assessed by comparing water quality improvement and the alternative IW value based on the alternative scenarios.