Net Ecosystem Production of a River Relying on Hydrology, Hydrodynamics and Water Quality Monitoring Stations

: Flow and water quality of rivers are highly dynamic. Water quantity and quality are subjected to simultaneous physical, chemical and biological processes making it di ﬃ cult to accurately assess lotic ecosystems. Our study investigated net ecosystem production (NEP) relying on high-frequency data of hydrology, hydrodynamics and water quality. The Kanawha River, West Virginia was investigated along 52.8 km to estimate NEP. Water quality data were collected along the river using three distributed multiprobe sondes that measured water temperature, dissolved oxygen, dissolved oxygen saturation, speciﬁc conductance, turbidity and ORP hourly for 71 days. Flows along the river were predicted by means of the hydrologic and hydrodynamic models in Hydrologic Simulation Program in Fortran (HSPF). It was found that urban local inﬂows were correlated with NEP. However, under hypoxic conditions, local inﬂows were correlated with speciﬁc conductance. Thus, our approach represents an e ﬀ ort for the systematic integration of data derived from models and ﬁeld measurements with the aim of providing an improved assessment of lotic ecosystems.


Introduction
Physicochemical and hydro-morphological properties of a river combined with biological communities [1,2] cause simultaneous physical, chemical and biological processes. Consequently, spatiotemporal variation is under an endless search for equilibrium by means of interactions between biotic and abiotic factors [3,4]. Those interactions can be assessed using the net ecosystem production (NEP) [5], which has been successfully applied to rivers for several decades [6]. NEP is the combined representation of gross primary production and ecosystem respiration. However, NEP must balance reaeration rates with photosynthetic production, respiratory consumption and all processes that can cause changes to dissolved oxygen [7][8][9]. Among the complex interaction of nutrients, biomass and trophic structure, NEP can be used, for example, as a means to explain carbon fluxes [10,11], given the potential of rivers to store, mineralize and transport carbon to coastal areas [9]. To accurately estimate NEP, data associated with flow and water quality must be reliable [12]. Nonetheless, different deployed instruments are needed. To succeed in an analysis over time and space, our study alternatively proposes the use of models to estimate flow as an adequate approach to reduce instrumentation. It also proposes an integration of hydrologic and hydrodynamic models with water quality data which can also be implemented as an automated procedure to estimate NEP. Development, calibration and validation of a hydrologic model to represent a drainage area along a river within the spatiotemporal domain must consider phenomena such as infiltration, evaporation and streamflow [13][14][15]. The hydrologic model took account of information about land use and topography specifications to handle hydro-climatic conditions in order to determine streamflow [16]. Consequently, streamflow can be coupled to a hydrodynamic model to predict flow along the river. In our study, the hydrologic and hydrodynamic models were implemented in the Hydrologic Simulation Program in Fortran (HSPF). Given that there was a set of unknown parameters governing the HSPF model, our study integrated a process for parameters calibration using the non-sorted genetic algorithm II (NSGA-II) [17,18]. This approach increased reliability under uncertainty of new hydrologic scenarios [19][20][21] and supported the use of an optimal solution [20,22,23] as the best set of parameters for the HSPF model. To guarantee satisfactory predictions, HSPEXP+ 2.0 was used to assess the calibrated HSPF model, as in the works of Xie et al. [20] and Lampert and Wu [24]. In this way, the streams module of the HSPF model can provide reliable predictions of the flow and velocity variables at specific locations along the river. Flow and velocity variables were input data to estimate reaeration rates. Nonetheless, an enhanced calculation of reaeration rates can be achieved by means of a set of equations [25] using a standardized Schmidt number [26,27] which was empirically estimated as a function of water temperature [28,29].
To provide an improved assessment of NEP, high-frequency data of water quality is now feasible [30], which seems promising and advantageous with respect to the periodic collection of water samples on a daily or lower frequency. Periodic observations of water quality at high-frequency could incorporate all processes such as cycles of nutrients (e.g., nitrogen, phosphorous), dissolved oxygen balance, sorption/desorption, volatilization, ionization, oxidation, biodegradation, hydrolysis and photolysis [31]. As result of all these processes, dissolved oxygen has been given special attention because it is the key variable for estimating NEP [32]. Dissolved oxygen is also closely related to water mixing, gas exchanges at the air-water interface, water temperature, flow, velocity and irradiance [27]. In addition, dissolved oxygen is subjected to spatial variability according to specifications of land use and local inflows [33][34][35][36][37]. Therefore, our study uses high-frequency data at various locations along the river as a way to contribute to the analysis of spatiotemporal impact of physicochemical properties of water on NEP.
This study was conducted within the Appalachian Region which is subjected to various water related stressors such as mining, urban settlements and industry. The Kanawha River, West Virginia, was chosen because it merges multiple inflows along the river at different rates and locations, disturbing water quantity and quality. Those inflows include tributaries, creeks, combined sewer overflows (CSO) and national pollutant discharge elimination system (NPDES). Therefore, NEP can be used as a proxy to assess those stressors and also to depict their variability over time and space.
In our study, hydrologic and hydrodynamic models and water quality data were integrated using a series of steps to estimate NEP under a high-frequency approach. Those steps were defined in this research as follows: (1) implementing a hydrologic model for drainage area along the river; (2) linking a hydrologic model with a hydrodynamic model; (3) collecting data about water temperature, dissolved oxygen and dissolved oxygen saturation, specific conductance, turbidity and ORP using monitoring stations installed along the river; (4) analyzing NEP under a spatiotemporal approach; and (5) assessing the impact of water quality on NEP and local inflows.

Watershed Description
The study area within the Appalachian Region was defined by 2995 km 2 draining water along the Kanawha River located in West Virginia, USA. Elk, Coal and Pocatalico Rivers ( Figure 1) are tributaries of the Kanawha River; these rivers delimited the drainage area at the location of the flow gages F2, F4 and F5, respectively. The Kanawha River had F1 and F3 flow gages located at the upstream and center of the study area. Locations for the five flow gages are provided in Table 1. Along the Kanawha River, the start and end limits were defined at the highest and lowest elevation of 190 m (38. (a) (b) Figure 1. Description of (a) the study area and (b) drainage area of the Kanawha River, West Virginia, which were delimited by flow gages and water quality monitoring stations specified in Table 1. Figure 1. Description of (a) the study area and (b) drainage area of the Kanawha River, West Virginia, which were delimited by flow gages and water quality monitoring stations specified in Table 1. Table 1. Location of flow gages, weather stations and water quality sondes indicated in Figure 1. Land topography is predominantly hilly dissected terrain. It contains forest, urban, barren and agricultural land representing 86%, 5.7%, 4.8% and 3%, respectively, where barren land was mainly characterized by mining activities, and the rest of the study area was dedicated to wetland. A detailed description of land use is shown in Table 2. Hourly discharge at start of the Kanawha River (F1) had minimum, average and maximum flows of 35.1, 387.2 and 3645 m 3 /s, respectively, during a period of observation of 913 days. The study area along the Kanawha River started with drainage areas at F1 and F2 of 21,680.8 and 2965.5 km 2 , respectively. Inflows of tributaries at F4 and F5 comprised drainage areas of 2232.6 and 616.4 km 2 , respectively. The study area (S1, S2 and S3) was focused on 2995 km 2 , which was bounded by the locations along the Kanawha River of the first flow gage (F1) and three sondes (Q1, Q2 and Q3). This study area provided conditions to conduct an analysis of the flow dynamics and impact of local inflows on water quality. The study area was characterized by annual average minimum and maximum dry bulb temperatures of 6.7 and 18.9 • C, respectively, an annual average precipitation of 1107 mm, and an annual average snowfall of 838 mm. Snow counted as precipitation, which was melted following the heat balance approach relying on precipitation, air temperature, solar radiation, wind velocity and dew point [13]. Climate data about precipitation, dry and dew point temperatures and wind speed were provided by two weather stations ( Figure 1). Solar radiation was not available in the two weather stations, so it was retrieved as an average for the study area from the national solar radiation database [38].

HSPF Model Description
The HSPF model was used for predicting flows of a watershed that joins modeling of watersheds and streams [39][40][41]. In addition, the HSPF model can be used to incorporate the transport of pollutants and nutrients. Information about land use, topography and climate assisted in the estimation of flow in streams ( Table 3) that follows a hydrodynamic pattern according to the terrain slope. Before water arrives at the streams, there are regular paths of water flow within the watershed which can be summarized in Figure 2. For instance, after rainfall there is an immediate interception by the canopy (CEPSC). The remaining water enters soil, which has a capacity to infiltrate and store within the upper zone (UZS) and lower zone (LZS). The excess of water could continue to reach active groundwater (AGWS) and base flow (BASETP) or enter deeper aquifers. were provided by two weather stations ( Figure 1). Solar radiation was not available in the two weather stations, so it was retrieved as an average for the study area from the national solar radiation database [38].

HSPF Model Description
The HSPF model was used for predicting flows of a watershed that joins modeling of watersheds and streams [39][40][41]. In addition, the HSPF model can be used to incorporate the transport of pollutants and nutrients. Information about land use, topography and climate assisted in the estimation of flow in streams that follows a hydrodynamic pattern according to the terrain slope. Before water arrives at the streams, there are regular paths of water flow within the watershed which can be summarized in Figure  2. For instance, after rainfall there is an immediate interception by the canopy (CEPSC). The remaining water enters soil, which has a capacity to infiltrate and store within the upper zone (UZS) and lower zone (LZS). The excess of water could continue to reach active groundwater (AGWS) and base flow (BASETP) or enter deeper aquifers. Concurrently, all of the watershed is subjected to evapotranspiration, where rates depend mostly on solar radiation, air temperature and humidity, whereas underground water is subjected only to evaporation. Each stage follows specific equations to determine flow rates which are defined by a set of parameters according to specifications of the watershed such as average slope (SLSUR) and mean elevation (MLS). Other parameters must be found during a process of adjustment such as groundwater recession flow (KVARY) and infiltration (INFILT). Descriptions of all parameters in the HSPF model are shown in Table 4, where seven values were deduced from the input data related to land use and topography and thirteen parameters must be calibrated in order to reduce error in flow predictions.
Modeling of streams relied on input data from inflows, river network configuration (Table 3) and a one dimensional approach of the river under a fully advective flow. The one dimensional approach required a homogeneous river transect, a representative Manning's coefficient (equal to 0.1) and a slope excluding dam elevation. Then, drainage sections with their corresponding streams were used to generate a HPSF model that was able to predict average flow and water velocity of streams. The governing equations used within the frame of a HSPF model can be found in Duda et al. [13]. Concurrently, all of the watershed is subjected to evapotranspiration, where rates depend mostly on solar radiation, air temperature and humidity, whereas underground water is subjected only to evaporation. Each stage follows specific equations to determine flow rates which are defined by a set of parameters according to specifications of the watershed such as average slope (SLSUR) and mean elevation (MLS). Other parameters must be found during a process of adjustment such as groundwater recession flow (KVARY) and infiltration (INFILT). Descriptions of all parameters in the HSPF model are shown in Table 4, where seven values were deduced from the input data related to land use and topography and thirteen parameters must be calibrated in order to reduce error in flow predictions.
Modeling of streams relied on input data from inflows, river network configuration and a one dimensional approach of the river under a fully advective flow. The one dimensional approach required a homogeneous river transect, a representative Manning's coefficient (equal to 0.1) and a slope excluding dam elevation. Then, drainage sections with their corresponding streams were used to generate a HPSF model that was able to predict average flow and water velocity of streams. The governing equations used within the frame of a HSPF model can be found in Duda et al. [13].
Air temperature below which evapotranspiration (ET) will arbitrarily be reduced below the value obtained from the input time series

Multiobjective Calibration of the HSPF Parameters
Among the various tools available in the search for adequate parameters defining the water dynamics of the HSPF model, we chose NSGA-II to look at the solution of two optimized objectives. The procedure to implement NSGA-II in MATLAB consisted of an iterative evaluation of different scenarios of the hydrologic model. Different scenarios were obtained using initial random values within the range stated in Table 3 regarding the parameters LZSN, INFILT, KVARY, AGWRC, DEEPFR,  BASETP, AGWETP, CEPSC, UZSN, INTFW, IRC, LZETP and NSUR. The iterative evaluation was accomplished for 400 sets of parameters that were evaluated in the HSPF model ( Figure 3). The next generation was deduced by creating 400 new sets of parameters that had a crossover and mutation probability of 0.9 and 0.1, respectively. The NSGA-II considered 1000 generations to define final calibrated values of the 13 parameters. To identify the best set of solutions, NSGA-II implemented two objectives: Nash-Sutcliffe model efficiency (NSE) and the percent bias coefficient (PBIAS) as the criteria to evaluate the error between flow measurements and HSPF model flow predictions. Identification of the optimal solution was accomplished by means of the Pareto front, which has the best solution when the magnitude is minimum for NSE [42] and PBIAS [43].

Field Measurements of Water Quality
The high-frequency monitoring system consisted of three Eureka Manta 2 multiprobe sondes. Each sonde measured water temperature (±0.1 °C), dissolved oxygen (±0.2 mg/L), dissolved oxygen saturation (±1%), specific conductance (±1 of reading), turbidity (±3% of reading) and oxidativereductive potential (±20 mV) with a time step of 1 h. Three locations along the river indicated in Figure 1, defined by Q1, Q2 and Q3, were monitored. That configuration facilitated our estimate of water quality changes as water moved downstream. Differences of sensor measurements served to further assess impact of local inflows comprising point and nonpoint sources of water pollution due to drained water along the river.

Net Ecosystem Production
NEP provides an assessment of rivers that encompasses physical and chemical characteristics. Physical characteristics include slope, width, depth and flow together with chemical characteristics such as nutrients, organic matter and water chemistry. In addition, other factors can be intrinsically intervening in NEP dynamics such as the effects of dams, riparian vegetation and pollution. NEP can also be seen as the balance of autotrophic and heterotrophic elements of the river [7]. Specifically, NEP can be evaluated through Equation (1) of Odum, 1956 [5].
where NEP is the gross primary production minus ecosystem respiration. k is oxygen reaeration coefficient. Cs is dissolved oxygen saturation and C is dissolved oxygen observed. P is the drainage accrual and accounts for all processes happening in the river together with dissolved oxygen of local inflows. Some of those processes include horizontal and vertical advection, photochemical oxidation of organic matter and nonaerobic consumption of oxygen during the time step of observation [7].
Replaces the initial set of parameters

Field Measurements of Water Quality
The high-frequency monitoring system consisted of three Eureka Manta 2 multiprobe sondes. Each sonde measured water temperature (±0.1 • C), dissolved oxygen (±0.2 mg/L), dissolved oxygen saturation (±1%), specific conductance (±1 of reading), turbidity (±3% of reading) and oxidative-reductive potential (±20 mV) with a time step of 1 h. Three locations along the river indicated in Figure 1, defined by Q1, Q2 and Q3, were monitored. That configuration facilitated our estimate of water quality changes as water moved downstream. Differences of sensor measurements served to further assess impact of local inflows comprising point and nonpoint sources of water pollution due to drained water along the river.

Net Ecosystem Production
NEP provides an assessment of rivers that encompasses physical and chemical characteristics. Physical characteristics include slope, width, depth and flow together with chemical characteristics such as nutrients, organic matter and water chemistry. In addition, other factors can be intrinsically intervening in NEP dynamics such as the effects of dams, riparian vegetation and pollution. NEP can also be seen as the balance of autotrophic and heterotrophic elements of the river [7]. Specifically, NEP can be evaluated through Equation (1) of Odum, 1956 [5].
where NEP is the gross primary production minus ecosystem respiration. k is oxygen reaeration coefficient. C s is dissolved oxygen saturation and C is dissolved oxygen observed. P is the drainage accrual and accounts for all processes happening in the river together with dissolved oxygen of local inflows. Some of those processes include horizontal and vertical advection, photochemical oxidation of organic matter and nonaerobic consumption of oxygen during the time step of observation [7]. Estimation of the k value was determined by means of k 600 (Equation (2)) which can be obtained by using one of three candidate Equations (4)-(6) and the Schmidt number. According to Raymond et al. [25], those three equations had the best fit with respect to field measurements. Those three equations also relied on the Schmidt number (Equation (3)) to estimate the mass transfer rates under momentum. The Schmidt number is the ratio of kinematic viscosity to the diffusion coefficient, which in turn can be determined as a function of temperature.
where V is water velocity, S is slope and H is depth of the river. g is the gravity force. Sc is the Schmidt number and T is temperature.

Input Data and Calibration
Data for all flow gages and climate stations were retrieved for the period from 1 October 2015 to 31 March 2018. To match all data, a common time step of one hour was adopted for all variables. Data from two climate stations were averaged instead of segmenting the watershed according to the area of influence, since both stations were in proximity. An example of average precipitation is shown in Figure 4. It should be noted that peaks related to precipitation might not coincide with peaks on measurements of flow gages due to local inflows of tributaries to the Kanawha River. The estimation of evapotranspiration rates was deduced by following the Turc method [44] and adding this data to the HSPF model. Climate data used in this model were compared by means of coefficient of variation (CV) with NASA data sources (i.e., NLDAS and AIRS) in daily time step for precipitation, temperature, dew-point temperature and evapotranspiration; results are presented in Table 5. Data of flow gages originally obtained with a time step of 15 minutes were converted to a 1 h time step using a moving average filter.  The calibration of the HSPF model was conducted for 2413 km 2 , corresponding to S1 and S2 drainage sections. The HSPF model was subjected to an iterative evaluation using the NSGA-II. The best optimal solution was deduced with the minimum Euclidean distance from the origin to the NSE and PBIAS scores. It should be pointed out that NSE and PBIAS were applied to the outflow F3 comprising inflows F1 and F2 together with predicted drained water at S1 and S2. Nonetheless, inflows F1 and F2 greatly contributed to the outflow predictions, given that the drainage area from F1 and F2 to F3 increased by 10%. This means that drainage changed from 24,646.3 to 27,060 km 2 at the location of gage F3. Such conditions enhanced NSE and PBIAS scores, which were 0.96 and 1.97%, respectively, when comparing predictions and observations from 1 October 2015 to 11 January 2018 in flow gage F3 ( Figure 5). NSE and PBIAS scores can be categorized as acceptable [43]; however, those results should be weighed based on the aggregated water between inflow and outflow of the drainage area along the river.
To confirm the adequacy of these parameters (Table 6) as the best optimal solution identified by NSGA-II, the calibrated HSPF model was analyzed by means of the HSPEXP+ 2.0 program in order to fulfill overall criteria regarding water balance (Table 7). It was found that error between predictions and measurements increased as the flow decreased; even so, the criteria were satisfied. Subsequently, the HSPF model predicted a water budget that was distributed as follows: 4.3% to surface flow, 17.9% to interflow, 32.5% to base flow and deep aquifers and 45.3% to evapotranspiration. The rate of evapotranspiration dominated water balance and was driven by BASETP and AGWETP and LZETP parameters. Water also accumulated in soil at rates determined by USZN, LZSN and INFILT parameters; however, a significant volume of water moved down to the base flow and deep aquifers.  The calibration of the HSPF model was conducted for 2413 km 2 , corresponding to S1 and S2 drainage sections. The HSPF model was subjected to an iterative evaluation using the NSGA-II. The best optimal solution was deduced with the minimum Euclidean distance from the origin to the NSE and PBIAS scores. It should be pointed out that NSE and PBIAS were applied to the outflow F3 comprising inflows F1 and F2 together with predicted drained water at S1 and S2. Nonetheless, inflows F1 and F2 greatly contributed to the outflow predictions, given that the drainage area from F1 and F2 to F3 increased by 10%. This means that drainage changed from 24,646.3 to 27,060 km 2 at the location of gage F3. Such conditions enhanced NSE and PBIAS scores, which were 0.96 and 1.97%, respectively, when comparing predictions and observations from 1 October 2015 to 11 January 2018 in flow gage F3 ( Figure 5). NSE and PBIAS scores can be categorized as acceptable [43]; however, those results should be weighed based on the aggregated water between inflow and outflow of the drainage area along the river. The calibration of the HSPF model was conducted for 2413 km 2 , corresponding to S1 and S2 drainage sections. The HSPF model was subjected to an iterative evaluation using the NSGA-II. The best optimal solution was deduced with the minimum Euclidean distance from the origin to the NSE and PBIAS scores. It should be pointed out that NSE and PBIAS were applied to the outflow F3 comprising inflows F1 and F2 together with predicted drained water at S1 and S2. Nonetheless, inflows F1 and F2 greatly contributed to the outflow predictions, given that the drainage area from F1 and F2 to F3 increased by 10%. This means that drainage changed from 24,646.3 to 27,060 km 2 at the location of gage F3. Such conditions enhanced NSE and PBIAS scores, which were 0.96 and 1.97%, respectively, when comparing predictions and observations from 1 October 2015 to 11 January 2018 in flow gage F3 ( Figure 5). NSE and PBIAS scores can be categorized as acceptable [43]; however, those results should be weighed based on the aggregated water between inflow and outflow of the drainage area along the river.
To confirm the adequacy of these parameters (Table 6) as the best optimal solution identified by NSGA-II, the calibrated HSPF model was analyzed by means of the HSPEXP+ 2.0 program in order to fulfill overall criteria regarding water balance (Table 7). It was found that error between predictions and measurements increased as the flow decreased; even so, the criteria were satisfied. Subsequently, the HSPF model predicted a water budget that was distributed as follows: 4.3% to surface flow, 17.9% to interflow, 32.5% to base flow and deep aquifers and 45.3% to evapotranspiration. The rate of evapotranspiration dominated water balance and was driven by BASETP and AGWETP and LZETP parameters. Water also accumulated in soil at rates determined by USZN, LZSN and INFILT parameters; however, a significant volume of water moved down to the base flow and deep aquifers.  To confirm the adequacy of these parameters (Table 6) as the best optimal solution identified by NSGA-II, the calibrated HSPF model was analyzed by means of the HSPEXP+ 2.0 program in order to fulfill overall criteria regarding water balance (Table 7). It was found that error between predictions and measurements increased as the flow decreased; even so, the criteria were satisfied. Subsequently, the HSPF model predicted a water budget that was distributed as follows: 4.3% to surface flow, 17.9% to interflow, 32.5% to base flow and deep aquifers and 45.3% to evapotranspiration. The rate of evapotranspiration dominated water balance and was driven by BASETP and AGWETP and LZETP parameters. Water also accumulated in soil at rates determined by USZN, LZSN and INFILT parameters; however, a significant volume of water moved down to the base flow and deep aquifers. Flow estimations using the calibrated HSPF model at the S1 and S2 drainage sections can be considered reliable as they were validated by flow measurements at gage F3. However, flow estimations at the outlet (Q3) of the S3 drainage area entirely relied on the accuracy of the calibrated HSPF model. The HSPF model validation found a significant contribution of the tributaries to the Kanawha River. For instance, from the total amount of water added within the S2 section, 78% of the water was contributed by the Elk River based on flow measurements at F2. In the same way, from the total amount of water added within the S3 section, 91% of the water was contributed by Coal and Pocatalico Rivers, according to measurements at F4 and F5. Flow dynamics at Q1, Q2 and Q3, during the period from 11 January 2018 to 31 March 2018 ( Figure 6), were based on the combined effects of inflows and drainage areas along the river. S1, S2 and S3 involved local inflows such as rainfalls, CSO and NPDES. The CV for flow data was computed having 0.77, 0.73 and 0.73 for locations Q1, Q2 and Q3, respectively. These CV values verified that flow dynamics were similar only in the Q2 and Q3 locations.  Flow estimations using the calibrated HSPF model at the S1 and S2 drainage sections can be considered reliable as they were validated by flow measurements at gage F3. However, flow estimations at the outlet (Q3) of the S3 drainage area entirely relied on the accuracy of the calibrated HSPF model. The HSPF model validation found a significant contribution of the tributaries to the Kanawha River. For instance, from the total amount of water added within the S2 section, 78% of the water was contributed by the Elk River based on flow measurements at F2. In the same way, from the total amount of water added within the S3 section, 91% of the water was contributed by Coal and Pocatalico Rivers, according to measurements at F4 and F5. Flow dynamics at Q1, Q2 and Q3, during the period from 11 January 2018 to 31 March 2018 (Figure 6), were based on the combined effects of inflows and drainage areas along the river. S1, S2 and S3 involved local inflows such as rainfalls, CSO and NPDES. The CV for flow data was computed having 0.77, 0.73 and 0.73 for locations Q1, Q2 and Q3, respectively. These CV values verified that flow dynamics were similar only in the Q2 and Q3 locations.

Water Quality
Water quality data were collected for 71 days with a time step of 1 h covering the period from 11 January 2018 to 22 March 2018 (Supplementary Data). Erroneous readings were discarded. Measurements for the three locations are shown in Figures 7 and 8 and available in Huber et al. [45]. The CV among all sensor readings (Table 8) showed that minimum and maximum scores were for dissolved oxygen saturation and turbidity, respectively. The same type of sensor readings among the three locations also had minimum and maximum differences of the CV for temperature and dissolved oxygen saturation. In summary, we found that dissolved oxygen saturation had minimum dispersion among all the sensors at the same location and maximum dispersion among the three locations.
Water 2020, 12, x FOR PEER REVIEW 11 of 18

Water Quality
Water quality data were collected for 71 days with a time step of 1 h covering the period from 11 January 2018 to 22 March 2018 (Supplementary Data). Erroneous readings were discarded. Measurements for the three locations are shown in Figures 7 and 8 and available in Huber et al. [45]. The CV among all sensor readings (Table 8) showed that minimum and maximum scores were for dissolved oxygen saturation and turbidity, respectively. The same type of sensor readings among the three locations also had minimum and maximum differences of the CV for temperature and dissolved oxygen saturation. In summary, we found that dissolved oxygen saturation had minimum dispersion among all the sensors at the same location and maximum dispersion among the three locations.

Water Quality
Water quality data were collected for 71 days with a time step of 1 h covering the period from 11 January 2018 to 22 March 2018 (Supplementary Data). Erroneous readings were discarded. Measurements for the three locations are shown in Figures 7 and 8 and available in Huber et al. [45]. The CV among all sensor readings (Table 8) showed that minimum and maximum scores were for dissolved oxygen saturation and turbidity, respectively. The same type of sensor readings among the three locations also had minimum and maximum differences of the CV for temperature and dissolved oxygen saturation. In summary, we found that dissolved oxygen saturation had minimum dispersion among all the sensors at the same location and maximum dispersion among the three locations.

Net Ecosystem Production
The reaeration rates (k) were calculated using Equations (4)-(6) that consequently helped to estimate NEP through Equation (1). The time series of hydrodynamics ( Figure 6) and water quality data ( Figure 7) were used at Q1, Q2 and Q3 locations to estimate NEP (Figure 9). Hydrodynamics of the Kanawha River showed higher flows and lower water velocities as water moved downstream, from Q1 to Q2 and to Q3. From water quality monitoring stations (Figure 7), a significant decay of dissolved oxygen and its saturation were observed in Q2 and Q3. For instance, a length of 23.5 km along the river, the distance between Q1 and Q2, had an average decay from 14.9 to 6.1 mg/L. In the following 29.3 km, the distance between Q2 and Q3, had an additional decay from 6.1 to 4.8 mg/L. Those dissolved oxygen decays reduced the NEP estimations from Q1 to Q2 by 93% and from Q2 to Q3 by 95%.
Water 2020, 12, x FOR PEER REVIEW 12 of 18

Net Ecosystem Production
The reaeration rates (k) were calculated using Equations (4)-(6) that consequently helped to estimate NEP through Equation (1). The time series of hydrodynamics ( Figure 6) and water quality data ( Figure 7) were used at Q1, Q2 and Q3 locations to estimate NEP (Figure 9). Hydrodynamics of the Kanawha River showed higher flows and lower water velocities as water moved downstream, from Q1 to Q2 and to Q3. From water quality monitoring stations (Figure 7), a significant decay of dissolved oxygen and its saturation were observed in Q2 and Q3. For instance, a length of 23.5 km along the river, the distance between Q1 and Q2, had an average decay from 14.9 to 6.1 mg/L. In the following 29.3 km, the distance between Q2 and Q3, had an additional decay from 6.1 to 4.8 mg/L. Those dissolved oxygen decays reduced the NEP estimations from Q1 to Q2 by 93% and from Q2 to Q3 by 95%.

Spatial and Temporal Variability of NEP
Repeatability and reproducibility are issues around NEP estimations when limited datasets are available in either time or space. These issues frequently happen over sub-daily patterns of dissolved oxygen [46]. In practice, NEP estimations are subjected to periods of observations, the choice of location along the river and the choice of the equation to estimate reaeration rates. Our study provided insight into NEP resolution over space and time based on field surveys and estimating reaeration rates through Equation (4).
For large rivers, either temporal or permanent local gradients may be observed as a result of flow regimes generating specific hydrodynamics [47,48] and consequently variability in NEP. Those NEP estimations will be the consequence of the level of turbulence in local mixing and the exchange rate of gases in the air-water boundary layer [49]. Nonetheless, NEP variability is also a consequence of physicochemical properties, such as organic matter [50], nutrient regimes [35], water temperature [51] and flow [10]. Our study conducted a river transect examination of NEP using horizontal and

Spatial and Temporal Variability of NEP
Repeatability and reproducibility are issues around NEP estimations when limited datasets are available in either time or space. These issues frequently happen over sub-daily patterns of dissolved oxygen [46]. In practice, NEP estimations are subjected to periods of observations, the choice of location along the river and the choice of the equation to estimate reaeration rates. Our study provided insight into NEP resolution over space and time based on field surveys and estimating reaeration rates through Equation (4).
For large rivers, either temporal or permanent local gradients may be observed as a result of flow regimes generating specific hydrodynamics [47,48] and consequently variability in NEP. Those NEP estimations will be the consequence of the level of turbulence in local mixing and the exchange rate of gases in the air-water boundary layer [49]. Nonetheless, NEP variability is also a consequence of physicochemical properties, such as organic matter [50], nutrient regimes [35], water temperature [51] and flow [10]. Our study conducted a river transect examination of NEP using horizontal and vertical profiles that occur near Q2. The horizontal profile was based on the Q2a and Q2b locations which are separated by 179 m, whereas the vertical profile was based on the Q2b and Q2c locations which are separated by 1.5 m. Three repetitions were conducted at each location, and water velocity was deduced from the HSPF streams module ( Table 9). The horizontal profile of NEP was not homogeneous since different conditions were observed in the field; however, there was a prevailing lower NEP in Q2b, which mainly occurred due to water temperature and dissolved oxygen measurements. In contrast, we found that the vertical profile generated lower NEP amounts in Q2c with respect to Q2b, which can be inferred as less prevailing irradiance, as it was 2.5 m deep. Still, such a discrepancy in NEP at the vertical profile was lower than the horizontal profile as a consequence of the distance between locations. Thus, these findings illustrated that spatial heterogeneity of NEP is driven by transport-reaction phenomena [52] due to local gradients created by hydrodynamics and their corresponding water quality. NEP in rivers is the result of a dynamic interaction between biotic [53] and abiotic factors [54,55]. Among those biotic factors, autotrophs and heterotrophs are continuously balanced to determine NEP dynamics through the year [56]. Autotrophy impacting NEP along the river is a consequence of nutrient loads from local inflows to the mainstream such as wastewater treatment plants, CSO and NPDES [57]. For instance, it has been found that spatial heterogeneity of NEP can be caused by watersheds comprising urban areas [58]. Our study had an urban area between Q2 and Q3, causing differentiated NEP estimations along the river (Figure 9). A driving variable causing a decay of NEP (Table 10) was dissolved oxygen, which mainly declined due to the various local inflows mixed with the mainstream. It can be interpreted that water residence times of 0.25 ± 0.03 days and 0.62 ± 0.12 days for the Q1-Q2 and Q2-Q3 river sections, respectively, along with local inflows with different water quality did not help to keep the same NEP estimations as observed in Q1. From our study, we can also claim that the balance between autotrophs and heterotrophs was significantly impaired as water moved downstream. We estimated a decay of 1.18 ± 0.38 g[O 2 ]/m 3 /day from NEP Q1 to NEP Q2 and an additional decay of 0.08 ± 0.12 g[O 2 ]/m 3 /day from NEP Q2 to NEP Q3 . The latter one was a consequence of prevailing hypoxic conditions observed in Q3.

Impact of Water Quality on NEP and Local Inflows
Additional water quality data were retrieved that can potentially affect NEP dynamics. In particular, the effects of specific conductance [59] and turbidity [60] on NEP were evaluated using the Spearman coefficient. We found through the Spearman coefficient that NEP and turbidity were positive at the three locations along the river (Table 11). We also found that NEP and specific conductance had negative values at the three locations along the river. Nevertheless, it can be deduced that specific conductance and turbidity could play a significant role in determining NEP if the dissolved oxygen measurements do not reach hypoxic conditions. Table 11. Spearman correlation of net ecosystem production (NEP) with water quality data at three locations of the Kanawha River. Because water quality and consequently NEP changed due to local inflows along the river, we computed volumes of aggregated water to the Kanawha River between paired locations W Q1-Q2 and W Q2-Q3 as well as their corresponding changes on NEP. We also conducted the same calculations for changes in specific conductance, turbidity and ORP measurements. For the aforementioned calculations, we considered average travel times obtained from the HSPF streams module. Then, Spearman correlations were calculated (Table 12). The segment of the river between Q1 and Q2 showed that W Q1-Q2 was mainly correlated with NEP and turbidity. For the segment of the river between Q2 and Q3, we found that W Q2-Q3 was mainly correlated with specific conductance and NEP. From these correlations, we can state that NEP can be used as an indicator to assess water quality of local inflows, as it merges various properties of the river related to hydrodynamics and water quality data along the river. However, a more reliable assessment could be achieved if hypoxic conditions are avoided.

Conclusions
In order to estimate flows along the Kanawha River, our study had to consider a drainage area of 2995 km 2 . A HSPF model was developed and then calibrated by means of NSGA-II in order to identify the best optimal solution. The streams module of the HSPF model served for hydrodynamic modeling, which provided data about flow and average water velocity at Q1, Q2 and Q3 locations along the Kanawha River. In addition, water quality data were collected for 71 days by placing sondes in the same three locations to hourly log dissolved oxygen concentration, dissolved oxygen saturation, water temperature, specific conductance, turbidity and ORP. Flow and average velocity data were used to estimate reaeration rates (k). Then, k values were used together with water quality data to estimate NEP. It was found that NEP greatly depends on the specific location within the river, as it was observed during a river transect examination. Our study also identified a decreasing NEP as water moved downstream, starting from NEP Q1 equal to 1.24 (±0.4) g[O 2 ]/m 3 /day to NEP Q2 equal to 0.09 (±0.07) g[O 2 ]/m 3 /day and to NEP Q3 equal to 0.004 (±0.03) g[O 2 ]/m 3 /day. Such decay was attributed to local inflows (W Q1-Q2 and W Q2-Q3 ), which were computed and correlated with their corresponding changes in water quality and NEP. The best Spearman coefficient (ρ = −0.71) was between W Q1-Q2 and NEP. However, under hypoxic conditions, the best Spearman coefficient (ρ = −0.61) was between W Q2-Q3 and specific conductance. These findings showed that spatial and temporal analyses of NEP were adequately addressed through datasets of hydrology, hydrodynamics and high-frequency data from water quality monitoring stations. Our study can also be useful for further research where assessment of local inflows to the mainstream should be accomplished by means of NEP. These advances encourage us to count more on field surveys, given that the scope of NEP dynamics in rivers depends on multiple scenarios related to flow and water quality conditions.