Decision-Making of LID-BMPs for Adaptive Water Management at the Boise River Watershed in a Changing Global Environment

We conducted a study on water management at the Boise River Watershed in a changing global environment potentially induced by climate variability and urbanization. Environmental ‘hotspots’ associated with water quality and quantity were first identified to select suitable management options, such as Low Impact Development (LID is commonly used for urban storm water management to reduce impacts induced by flash flood in urban environment while improving water quality standard by filtering non-point source pollutants from predominant, impervious land segments in urban settings.) and Best Management Practices (BMPs) for urban and rural land segments, respectively. A decision-making process was employed to evaluate the cost-effectiveness for each management option based on multiple criteria, including water quality, financial challenges, and other environmental concerns. The results show that LID/BMPs were useful to control water quality in the watershed. The effectiveness of LID/BMPs implementation was subject to change with the placement location and consideration objectives associated with economic or environmental aspects. It appears that about 10% of the study area is required to implement water management options (LID/BMP) to improve water quality potentially driven by climate variability and urbanization. We anticipate that this study will make a case toward developing a sustainable water management plan in a changing global environment, especially for the urban–rural interface settings.


Introduction
Water quality is one of the significant environmental issues at the urban-rural interface, such as Boise in the United States (US) [1]. Changing the global environment driven by climate variability, urbanization, and population growth further affects water quality standards across the states. Nutrient loadings from urban runoffs and agricultural return flows create taste and order problems and disrupt recreation activities [2,3]. Population growth and economic development overexploit surface and groundwater resources [4,5], yet it is challenging to manage water resources due to the complicated nature and human dynamics. Exploring adaptive water management, therefore, is a critical exercise to mitigate water quality (WQ) impacts in this changing global environment. Adaptive water management represents systematic processes to improve water management practices and policies by recognition of the importance of the human and social dimensions of resource use.

Study Area
The Boise River Watershed (BRW) was selected for the study area ( Figure 1). The BRW is 10,439 km 2 with a mainstream length of a 164 km stretch and flows into the Snake River near Parma. The mean annual rainfall in this watershed is 274 mm, and the average annual temperature has been 13.70 • C from 1979 to 2016. More than 40% of Idaho's population lives in the greater Boise metropolitan areas, including Boise, Nampa, Meridian, and Caldwell. Concerns of water quality degradation, however, have been discussed in the scientific communities due to recent urbanization and climate variability.
Water 2020, 12, x FOR PEER REVIEW 3 of 17

Study Area
The Boise River Watershed (BRW) was selected for the study area ( Figure 1). The BRW is 10,439 km 2 with a mainstream length of a 164 km stretch and flows into the Snake River near Parma. The mean annual rainfall in this watershed is 274 mm, and the average annual temperature has been 13.70 °C from 1979 to 2016. More than 40% of Idaho's population lives in the greater Boise metropolitan areas, including Boise, Nampa, Meridian, and Caldwell. Concerns of water quality degradation, however, have been discussed in the scientific communities due to recent urbanization and climate variability.

Hydrological Simulation Program-Fortran (HSPF) Model
For a hydrological model, the Hydrological Simulation Program-Fortran (HSPF) from the previous work [25] was selected and applied for this study. HSPF is a watershed scale, process-based, and semi-distributed model. This model can effectively simulate streamflow and water quality associated with land management practices and climate variability at the urban-rural interface, such as BRW [26][27][28][29][30]. The HSPF model consists of the main three modules (PERLND, IMPLND, and RCHERS). Each module has different parameter sets for hydrology, water quality components, and state variables [31]. For HSPF model evaluation, Nash-Sutcliffe Efficiency (NSE), correlation coefficient (R), Mean Square Error-observation standard deviation ratio (RSR), and percentage of bias (PBIAS) are computed for model calibration and validation (see Appendix A).

Meteorological Data
Phase 2 of the North American Land Data Assimilation System (NLDAS-2) data was used to retrieve climate forcing data, including precipitation, temperature, wind speed, solar radiation. The NLDAS-2 is available from the hourly to yearly time step with 1/8th-degree grid spacing for the simulation period from 1 January 1979 to 31 December 2015. Cloud cover data, however, were

Hydrological Simulation Program-Fortran (HSPF) Model
For a hydrological model, the Hydrological Simulation Program-Fortran (HSPF) from the previous work [25] was selected and applied for this study. HSPF is a watershed scale, process-based, and semi-distributed model. This model can effectively simulate streamflow and water quality associated with land management practices and climate variability at the urban-rural interface, such as BRW [26][27][28][29][30]. The HSPF model consists of the main three modules (PERLND, IMPLND, and RCHERS). Each module has different parameter sets for hydrology, water quality components, and state variables [31]. For HSPF model evaluation, Nash-Sutcliffe Efficiency (NSE), correlation coefficient (R), Mean Square Error-observation standard deviation ratio (RSR), and percentage of bias (PBIAS) are computed for model calibration and validation (see Appendix A).

Meteorological Data
Phase 2 of the North American Land Data Assimilation System (NLDAS-2) data was used to retrieve climate forcing data, including precipitation, temperature, wind speed, solar radiation. The NLDAS-2 is available from the hourly to yearly time step with 1/8th-degree grid spacing for the simulation period from 1 January 1979 to 31 December 2015. Cloud cover data, however, were derived separately from Climate Data Online (CDO) provided by the National Climatic Data Center (NCDC), because they are not included in NLDAS-2. Dew-point data, T dew , are also calculated separately using an empirical equation [32] based on grid points of NLDAS-2 within the BRW.
where T dew is calculated dew point ( • C), T is air temperature ( • C), A, m, and T n are constants as referred to [32]. P ws is water vapor saturation pressure over water (hPa). Penman-Monteith's equation below [33] is used to compute PET.

Geospatial Data
The Digital Elevation Model (DEM) at 30-m resolution [34] was used to generate topographic relief images, watershed delineation files, and flow directions via automatic watershed delineation processes [35]. The stream network can be derived from DEM and National Hydrography Dataset (NHD) [34] to generate a detailed stream network in higher spatial resolution (1:100,000). For environmental background data, the Land Use Land Cover (LULC) is processed to classify its land segments with respect to the delineated 78 sub-watersheds.

Global Circulation Model (GCM) and Land Use Land Cover (LULC) Data
Future climate data were imported from the Canesm2 model, which represents the best historical climate condition at the BRW [25]. The Cansesm2 model is one of the Canadian CMIP5 models; it combines the physical coupled atmosphere-ocean model (CanCM4) and a terrestrial carbon model (CTEM) based on the Canadian Terrestrial Ecosystem Model (CTEM). Since the Canesm2 model showed that the relatively small bias between the simulated streamflow and the observed flows, which is less than 10%, Canesm2 under Regional Climate Prediction 8.5 climate scenarios (RCP 8.5) was used for HSPF inputs as climate forcing, because it represented severe future climate conditions in the study area. For future LULC conditions, 2100 LULC products developed by the United States Geological Survey (USGS) using the forecast scenarios of land use change (FORE-SCE) model were used [36] with Intergovernmental Panel on Climate Change (IPCC) the Special Report on Emissions Scenarios (SERS) (e.g., A1B, A2, and B1) associated with the global/regional economic, technological, and environmental cooperation, and economic growth. The 2100 LULC coverage, which is a polygon shape in vector format, was reclassified to represent urban, agricultural land, forest land, water/wetland, shrubland, grassland, and barren/mining land to be incorporated into HSPF. To consider rigorous urbanization in the future land use, 2100 LULC under A2 scenario was selected because urban land in projected 2100 LULC under A2 scenario indicates the greatest increase in urban land from 5.40% to 7.92% of whole watershed areas, because it presumes high economic growth and very high population growth globally relative to 2011 urban land (see Table 2). In general, agricultural land and shrubland are converted to urban land when comparing with 2011 LULC. Barren/mining land sharply increased from 20 km 2 in 2011 to 470 km 2 in 2100 under the A2 scenario, which is a 4.31% increase for whole land use relative to 2011 LULC due to the conversion of grasslands and forest. Note that [37] showed that grassland degradation and forest deforestation result in the highest increase in barren/mining land.

BEO-Parameter ESTimation (BEOPEST)
An automatic calibration tool, the BEO-Parameter ESTimation (BEOPEST), was applied to calibrate streamflow at six calibration target points using key hydrological parameters (see Table 1). The BEOPEST is a tool to mitigate the computation burden and implement parallelism in the model-independent nonlinear parameter estimation and optimization tool developed by Doherty and Skahill (2006) [37]. The PEST uses a recursive gradient-based optimization technique, linearizing the nonlinear problem by iteratively computing the Jacobian matrix of sensitivities of model observations to parameters. The parameter estimation in PEST is accomplished using the Gauss-Marquardt-Levenberg algorithm (GML) to minimize the user-defined objective function (e.g., minimization of root mean squares between simulated and observed values) [38]. The detailed information of computer parallelism in the HSPF is available in the literature [28].

Watershed Management Tool (WMT)
A watershed management tool was applied to identify critical hotspots (CHSs). Three indices, including Concentration Impact Index (CII), Load Impact Index (LII), and Load per sub-basin area Index (LPSAI) are first used to sort out the final ranked CHSs [39]. The priority scheme (high, medium, and low) with respect to each index is then applied to identify the final CHSs, which meet 10% acceptance level [40]. The most suitable LID/BMP choice is then considered to implement at the urban (impervious dominant areas) and rural land (previous dominant areas) segments, respectively.
Five LID and BMPs options, including (1) Bioretention-LID, (2) filter strip-BMP, (3) grassed swale-LID, (4) wetland-BMP, and (5) detention pond-LID are considered applying at the selected CHSs. The removal efficiency of each LID/BMP is computed using the BMP efficiency editor tool available at winHSPF 3.0 [41]. The reduction rate of sediment, TN, and TP loads is calculated using HSPF model outputs before and after LID/BMP implementation. To determine the best LID/BMP choice, three different schemes of (1) equal weight to each three-factor (EE), (2) more weight on the economic feasibility (EFW), and (3) more weight on environmental concerns (EW) along with specific decision factors were applied, as shown in Table 3. For example, the equal distribution scheme treats all three Water 2020, 12, 2436 7 of 17 decision factors, including water quality (WQ), operational cost (OC), land feasibility (LF) equally with the same weight vector. Thus, one-third (33.33%) of the decision factor is applied to OC and LF, while another one-third is assigned to WQ where the same weight (11%) is applied to each water quality element (Sediment, TN, TP). Two other schemes with more weights on the economic feasibility (EFW) and more weights on environmental concerns (EW) are applied. Unlike water quality, a scheme on the EFW utilizes a different weight vector to emphasize the financial burden when BMP/LID is in place, while the scheme with EW weights more on water quality (WQ) improvement (see Table 3). The MCDA method is then performed to assign the best suitable LID/BMP option for the identified CHSs in sub-watershed scales using the vectors listed in Table 3. The data layers are assessed with one another through a pairwise comparison matrix with regard to their significance by applying the Analytical Hierarchy Process (AHP) technique. AHP is a decision-making process to determine the optimal solution with different weighting [42,43]. It can be employed to formulate and solve the problem hierarchically [44]. Basically, AHP is a hierarchical structure with the goal and criteria, which can be divided into sub-criteria for further analysis. The pairwise comparisons matrix determines preference or the overall priorities for the final rank using the normalized principal priority vector (Eigenvector). From this matrix, a consistency check is performed, while the computed Consistency Index (CI) relative to the eigenvalue should be less than 10% for satisfactory results [40], as defined by: where r max is the maximal eigenvalue, n is the number of rows or columns in the pairwise matrix, CR is the consistency ratio (CR), which should be less than 10% for the acceptable consistency, and RI is the random index [40]. The total cost for LID/BMP implementations and cost-per-unit pollutant load reduction (e.g., sediment, TN, TP loads) are estimated as shown in Table 4. Total cost is the sum of construction and maintenance cost for the selected LID/BMP choice. Note that cost-per-unit pollutant load reduction is used to calculate cost-per-kg pollutant reduction for sediment, TN, and TP, on average.

HSPF Calibration
A total of six calibration target points were calibrated from 2001 to 2015 and validated from 1981 to 2000 for streamflow. The initial two years were applied as warm-up periods for calibration (1999 to 2000) and validation (1979 to 1980). Water quality calibration for sediment, TN, and TP were only conducted at the watershed outlet due to the limited observed water quality data. Figure 2 shows hydrographs' comparison of the monthly observed and simulated streamflow at the calibration target points, and Table 5 shows the performance measures of water quality calibration at the watershed outlet. Overall, the simulated streamflow closely matched the observed flows. The magnitude of peak flow indicates somewhat different results owing to not only water diversion for irrigation near calibration target point 5 and 6, but also inherent model limitation from reservoir operation rules at the calibration target point 5. Based on the recommend statistics of hydrological model performance at the monthly time step [50], HSPF shows the acceptable performance for calibration and validation at the watershed outlet based on R values for streamflow (above 0.70), sediment load (above 0.8), and nutrient loads (above 0.8) (see Table 5). NSE, RSR, and PBIAS statistics also indicate the fair and good model performance (streamflow, sediment load, and TN load) after calibration.   . Future precipitation tends to increase substantially in summer by 69.16% and then followed by winter with 34.14% for the F3 period as opposed to the baseline condition. Mean annual precipitation, however, decreases by 2.51% for F1, while it increases by 8.30% and 27.85% for F2 and F3, respectively. Future temperature shows the increasing trend for all three periods as opposed to the baseline. A consistent increase in temperature is observed throughout the future periods under the RCP8.5 scenario. Especially, summer and winter temperatures highly increase over time, while mean annual temperature shows somewhat of a variation (the increasing trend in 3.03 • C, 4.83 • C, and 6.47 • C for F1, F2, and F3, respectively). flow indicates somewhat different results owing to not only water diversion for irrigation near calibration target point 5 and 6, but also inherent model limitation from reservoir operation rules at the calibration target point 5. Based on the recommend statistics of hydrological model performance at the monthly time step [50], HSPF shows the acceptable performance for calibration and validation at the watershed outlet based on R values for streamflow (above 0.70), sediment load (above 0.8), and nutrient loads (above 0.8) (see Table 5). NSE, RSR, and PBIAS statistics also indicate the fair and good model performance (streamflow, sediment load, and TN load) after calibration.    Figure 3 shows the projected changes in seasonal precipitation and temperature for the three future time windows (F1: 2021-2045, F2: 2046-2070, F3: 2071-2095) relative to the baseline period . Future precipitation tends to increase substantially in summer by 69.16% and then followed by winter with 34.14% for the F3 period as opposed to the baseline condition. Mean annual precipitation, however, decreases by 2.51% for F1, while it increases by 8.30% and 27.85% for F2 and F3, respectively. Future temperature shows the increasing trend for all three periods as opposed to the baseline. A consistent increase in temperature is observed throughout the future periods under the RCP8.5 scenario. Especially, summer and winter temperatures highly increase over time, while mean annual temperature shows somewhat of a variation (the increasing trend in 3.03 °C, 4.83 °C, and 6.47 °C for F1, F2, and F3, respectively).   Table 6 indicates annual mean variations of streamflow and water quality loads induced by the combined climate variability under RCP8.5 and 2100 LULC associated with the A2 scenario after HSPF calibrations. Annual mean streamflow decreases by 24.14% and 11.19% for F1 and F2, respectively, while it increases by 0.98% for F3 when the combined climate and land use changes are taken into account. The relative change in annual mean sediment load, TN load, and TP load associated with climate variability and land use change show the increasing trend over time for F1, F2, and F3. The period of F3, in particular, shows the most substantial increase in sediment load, TN, and TP by 49.30%, 84.77%, and 21.16%, respectively, as opposed to the baseline. The simulated results  Table 6 indicates annual mean variations of streamflow and water quality loads induced by the combined climate variability under RCP8.5 and 2100 LULC associated with the A2 scenario after HSPF calibrations. Annual mean streamflow decreases by 24.14% and 11.19% for F1 and F2, respectively, while it increases by 0.98% for F3 when the combined climate and land use changes are taken into account. The relative change in annual mean sediment load, TN load, and TP load associated with climate variability and land use change show the increasing trend over time for F1, F2, and F3. The period of F3, in particular, shows the most substantial increase in sediment load, TN, and TP by 49.30%, 84.77%, and 21.16%, respectively, as opposed to the baseline. The simulated results under F3, therefore, are used for further analysis to evaluate water management in the study area.

The Identified Critical Hotspots (CHSs)
Three indices, including Concentration Impact Index (CII), Load Impact Index (LII), and Load per sub-basin area Index (LPSAI) are used to select critical hotspots (CHSs) in the BRW. As shown in Figure 4, the CII method selects 16 sub-watersheds indicating CHSs (low-, medium-, and high-priority) that are mainly located in the lower Boise watershed below the Lucky Peak dam. It is expected that potential land use change from forestry to marginal/agricultural lands might contribute additional nutrient load (nitrogen) to the Boise River. Similarly, the identified CHSs using the LII method are mostly located in the lower Boise Watershed near the Lucky Peak dam ( Figure 5). Perhaps, sediment deposition to the reservoir results in water quality impairment in the middle segment of the Boise River. Unlike CII and LII, LPSAI's results, however, show very differently in the sense that the identified CHSs are located sparsely across the watershed (see Figure 6). Since the LPSAI method classifies the sub-watersheds with a relatively large amount of pollutants in a small area, significant loads are possibly inherited from upstream to confluences where nutrient converges from small stream segments.

The Suitable LID/BMP Choice for CHSs
For the identified CHSs above, a decision process to select the best suitable LID/BMP option was carried out based on three different management schemes (EE, EFW, EW) using the AHP structure. The pairwise comparison matrices for CII, LII, and LPSAI are generated as a function of the reduction rate of sediment and nutrient concentration at each sub-watershed level. Table 7 shows one example of the pairwise comparison matrix for sub-watershed 78 based on sediment reduction using the CII method (not showing all results in this paper). The pairwise comparison matrices of CR (Equation (5)) are calculated 0.091%, 0.001%, 0.01% for sediment, TP, and TN, respectively. Since these CR values are less than 10%, this performance is satisfactory [40]. The weight vector of land feasibility varies depending on slope, drainage area, and land type for an individual sub-watershed. Table 8 shows the selected CHSs associated with their suitable LID/BMP options using three different indices. For example, if a water manager is particularly interested in nutrient concentration (CII) rather than nutrient load (LII or LPSAI), but financial challenge (FC) is inevitable, he or she can select only the filter strip option to be installed at the respective sub-watershed (see Table 8). In contrast, bioretention-LID would be preferable if both nutrient concentration (CII) and nutrient load (LII or LPSAI) are critical, because the environmental concern is evolving in their community. Basically, the water manager can pick and choose the best LID/BMP option to meet their needs.
Perhaps, sediment deposition to the reservoir results in water quality impairment in the middle segment of the Boise River. Unlike CII and LII, LPSAI's results, however, show very differently in the sense that the identified CHSs are located sparsely across the watershed (see Figure 6). Since the LPSAI method classifies the sub-watersheds with a relatively large amount of pollutants in a small area, significant loads are possibly inherited from upstream to confluences where nutrient converges from small stream segments.

The Suitable LID/BMP Choice for CHSs
For the identified CHSs above, a decision process to select the best suitable LID/BMP option was carried out based on three different management schemes (EE, EFW, EW) using the AHP structure. The pairwise comparison matrices for CII, LII, and LPSAI are generated as a function of the reduction rate of sediment and nutrient concentration at each sub-watershed level. Table 7 shows one example

Evaluation of LID/BMP Choice
The cost-effectiveness of the management scheme associated with three different indices is summarized in Table 9 to provide useful insights. The pollutant load reduction and efficiency of LID/BMPs vary among their types and location, so the pollutant load reduction per unit area (kg/ha) and total cost per unit pollutant load reduction ($1000/kg) associated with evaluation criteria are great of interest. Overall, filter strip-BMP is the most economical option at the sub-watershed levels compared with other LID/BMP options, while the detention pond-LID and bioretention-LID options are identified as intermediate effectiveness of LID/BMPs due to higher construction cost and pollutant load reduction rate. Thus, when the EE scheme is applied, the results show that LPSAI is the greatest reduction of the average pollutant load by 4.21 kg/ha, while CII is the most cost-effective by 195,728 dollars/kg due to the lower total cost-per-unit pollutant load reduction for TN and TP. Although LPSAI is the greatest reduction of the average pollutant in this particular case, the final decision should be made by the water manager or other decision-maker to select the best LID/BMP out of three options, including Grassed swale-LID, Filter strip-BMP, and Detention pond-LID (see Tables 8 and 9). Similarly, when LPSAI is used, the water manager considers that environmental concern is a higher priority so that he/she can select either detention pond-LID or filter strip-BMP. If that is the case, the average pollutant load reduction can be achieved by 3.76 kg/ha, while the cost-effectiveness is reported as $234,736/kg.

Summary and Conclusions
The study investigated the potential consequences of climate variability and urbanization associated with land use change on water quality and quantity at the Boise River Watershed over the next few decades. To adapt such consequences, potential water management options, such as LID and BMP, are considered to be implemented in urban (impervious land segment) and rural (pervious land segment) settings. A decision-making process was then proceeded to identify CHSs, to select proper LID/BMP at CHS, and to evaluate the selected LID/BMP options associated with decision factors, including water quality, operational costs, and land feasibility.
A spatial target method using three different indices (CII, LII, and LPSAL) and AHP was then used for CHS identification and LID/BMP selection, respectively. About 8 sub-watersheds out of 78 watersheds are identified as high-priority CHSs using each index, and suitable LID/BMP options are selected and applied accordingly based on different management schemes (economic feasibility, environmental concern). Each LID/BMP option is then evaluated to provide useful insights for water managers to adapt potential consequences in a changing global environment (climate variability and urbanization).
The result shows that only filter strip-BMP is the best option when the economic feasibility is the most critical at the Boise River Watershed, regardless of index values. However, when the environmental concern is the best interest of multiple stakeholders, a water manager may choose the detention pond-LID option as well. It appears that more flexible LID/BMP options are available when all three determining factors (water quality, operational cost, and land feasibility) are treated equally and obviously; a water manager will have more degree of freedom during his decision processes if that is the case. Therefore, these results can strengthen science-based decision-making and provide useful insights to policy-makers in cost-effective water resource management. We anticipate that this informative study can build a case toward developing adaptive water management at the urban-rural interface in a changing global environment. However, in the future, it will be necessary to improve WMT by expanding the number of scenarios through conducting a survey of farmers, water managers, or policy-makers in the watershed to reflect more reliable criteria and their preference for LID/BMPs implementation.

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

Appendix A
where Q Oi and Q Si are observed and simulated streamflow at the time step, respectively. Q Oi and Q Si are mean observed and simulated streamflow for the simulation period. N is the total number of values within the simulation period. R is the correlation coefficient between the predicted and observed values. It ranges from 0.0 to 1.0. A higher value indicates better agreement between predicted and observed data. Ref. [51] indicated that R values greater than 0.7 show acceptable model performance. NSE is the percentage of the observed variance and determines the efficiency criterion for model verification [52]. It is calculated from minus infinity to 1.0. Higher positive values indicate better agreement between observed and simulated values. RSR is a standardized Root Mean Square Error (RMSE) based on observed standard deviation recommended by [53]. The zero value shows the optimal model performance. PBIAS calculates the average tendency of the simulated values to be larger or smaller than observed counterparts [54]. A lower PBIAS value (e.g., close to zero) indicates better performance. A positive PBIAS indicates underestimated bias, while negative PBIAS values show the overestimated bias.