© EDP Sciences, 2020

1 Introduction

Stream metabolism represents the structural and functional characteristics of a stream ecosystem. It is one of the fundamental indicators of nutrient and organic matter cycling and stream health, based on physical, chemical and biological properties of streams (hydrology, geomorphology, climate, water chemistry, aquatic and riparian vegetation, etc.) (Mulholland et al., 2005; Izagirre et al., 2008; Williamson et al., 2008; Bernot et al., 2010). Gross primary production (GPP), net ecosystem metabolism (NEM) and community respiration (R_c) are the main components of stream metabolism. Primary productivity is defined as the rate of which organic matter (biomass) is formed from inorganic carbon by photosynthesizing organisms, thus representing the conversion of solar energy to reduced chemical energy (Bott, 2006).

Diel oxygen techniques (DOTs) have been commonly used to estimate metabolism components in streams, lakes and lagoons and are based on measuring diel variations of dissolved oxygen (DO) concentrations due to photosynthesis, respiration and reaeration in the water bodies (Odum, 1956; Seeley, 1969; Wang et al., 2003; Vallino et al., 2005; Staehr and Sand-Jensen, 2007; Van de Bogert et al., 2007; Ciavatta et al., 2008; Hanson et al., 2008; Staehr et al., 2010; Karakaya et al., 2011). In the related literature, different DOTs such as whole-stream method (also known as open channel diel oxygen method: single or two-station method) and chamber methods have been used to measure the stream metabolism affected by land use (agriculture, forest, etc.), livestock production, groundwater inputs, etc. (Ganning and Wulff, 1970; Grimm and Fisher, 1984; Marzolf et al., 1994; Young and Huryn, 1998; Fellows et al., 2001; Mulholland et al., 2001; Hall and Tank, 2005; Houser et al., 2005; Roberts et al., 2007; Gücker et al., 2009; Bernot et al., 2010; Riley and Dodds, 2013; Yates et al., 2013; Roley et al., 2014; Houser et al., 2015). Effects of wastewater treatment plant (WTP) discharges on stream metabolism are also important as WTPs are prevalent across the world and their discharges are sources of nutrients and organic carbon. There are many studies related to the influence of WTP discharges, such as Kicklighter (1987), Gücker et al. (2006), Ruggiero et al. (2006), Sanchez-Perez et al. (2009), Wassenaar et al. (2010), Chen (2013), Aristi et al. (2015), Chesworth (2016). They concluded that WTP discharges increased R_c while different GPP results were obtained. For example, Gücker et al. (2006) found that GPP increased after the discharges, while Kicklighter (1987) determined as it decreased. Sunlight absorbed by the high density wastewater was indicated as the reason of these different results of GPP in the related literature.

In this regard, the first purpose of this study was to investigate the effects of the urban WTP discharges which consist of both treated and untreated raw sewage on the metabolism of Büyüksu Stream (Bolu, Turkey). The previous studies only focused on the effects of WTP discharges, land use, etc. on stream metabolism by using DOTs and there are no comprehensively model studies between metabolism rates and environmental variables. Thus, the second purpose was to model the metabolism components (GPP, NEM, and R_c) as a function of measured environmental variables including air temperature (T_air), atmospheric pressure (P_atm), relative humidity (RH) and water quality variables. Within this first study related to the effects of WTP on the stream metabolism in Turkey, both the results were compared with the previous studies and it was sought to fill the gap in the statistical modelling of metabolism rates.

2 Material and method

2.1 Study site

Büyüksu Stream is located in Bolu Province in the western Black Sea region of Turkey, is fed by the two main tributaries of Abant Creek and Mudurnu Creek, and has a drainage catchment area of 1112.5 km² (see Fig. 1). The basin of Büyüksu Stream consists of approximately 3.1% residential and urban areas, 28.3% agricultural areas, 0.24% wetlands, and 68.36% forests, prairies, etc. According to the long-term meteorological data between 1927 and 2016, Bolu has a cool temperate climate with snowy winters and warm summers with cool nights. The mean annual temperature is 10.5 °C, the mean annual maximum and minimum temperature are 17.1 °C (max. 39.8 °C) and 4.5 °C (min. −34 °C), respectively. The mean annual precipitation is 545.3 mm, while the mean annual number of days with precipitation is 137.7, and the mean annual sunshine hours are 65.6 h (total of mean daily hours of every month) (Turkish State Meteorological Service, 2017).

Bolu City WTP (Activated Sludge Process) has two separate and continuous discharges, which are treated and untreated, received by Büyüksu Stream. Because the capacity of the WTP is not sufficient to treat all the city's wastewaters (The WTP can only treat 50% of the entire wastewater). The distance between two discharges is approximately 15 m.

Measurements were carried out in the following four reaches: before (BD) and after discharges (AD) of the WTP in Büyüksu Stream, and the springs of Mudurnu Creek and Abant Creek between August 2015 and December 2016. The springs (headwaters) of Mudurnu Creek and Abant Creek were sampled to represent reference conditions for water quality in the basin. Coordinates and locations of four reaches were given in Table 1 and Figure 1, respectively. The springs of Abant Creek and Mudurnu Creek are approximately 41.2 km and 31.4 km away to BD reach, respectively. The distance between BD and AD reaches is 0.5 km, and AD reach is 0.2 km downstream of the WTP discharges.

Fig. 1 Locations of the four reaches in the basin of Büyüksu Stream.

2.2 Measurements of environmental variables

Metabolism components were estimated using the two-station method developed by Odum (1956) for 1–6 days in every month between August 2015 and December 2016 (17 months). DO (mg L⁻¹) and water temperature (T_w, °C) measurements were performed at the upstream (US) and downstream (DS) of each stream reach with one-minute intervals (dt = 1 min) for at least 36 hours (2–7 days) by using oxygen data loggers (MiniDOT, PME, Vista, CA, USA). The reach length was selected as 150 m for each study reach, according to Bales and Nardi (2007). The data loggers were placed in protection cages throughout the measurements. While measuring DO and T_w, air temperature (T_air, °C), atmospheric pressure (P_atm, bar) and relative humidity (RH, %) were simultaneously measured (dt = 1 min) by using data loggers (RHT50, Extech Instruments, USA).

Water samples were collected at 15-minute intervals during two hours for composite samples (one composite sample consists of 9 water samples mixed) representing average water quality while stream flow rate (Q, m³ s⁻¹), stream velocity (V, m s⁻¹), stream depth (D, m), and stream width (W, m) were measured by using an acoustic doppler velocimeter (SonTek FlowTracker Handheld ADV, California, USA) on the day of both deployment and collection of the DO loggers in each reach and every month. Linear interpolations based on time were applied for Q, V, and D values measured at the start and end of the deployment of the DO loggers instead of average of two measurements for more accurate estimations of GPP, NEM, and R_c.

In the water samples, pH, specific conductivity (SC, µS cm⁻¹), biochemical oxygen demand (BOD₅, mg L⁻¹), total nitrogen (TN, mg L⁻¹), total phosphorus (TP, mg L⁻¹), orthophosphate (ortho-PO₄-P, mg L⁻¹), ammonium nitrogen (NH₄-N, mg L⁻¹), nitrate nitrogen (NO₃-N, mg L⁻¹), chlorophyll a (chl-a, µg L⁻¹) and turbidity (NTU) were measured. pH and SC were measured by using a multi-parameter probe (Hach HQ40d portable meter, Hach Company, Loveland, CO, USA). BOD₅ was measured by using a respirometric pressure system (WTW Oxitop IS6, Germany). TN was measured by Hach LCK 138 Laton cuvette test (ınorganically and organically bonded nitrogen is oxidized to nitrate by digestion with peroxo-disulphate, the nitrate ions react with 2.6-dimethylphenol in a solution of sulphuric and phosphoric acid to form a nitrophenol). TP was measured by Hach LCK 349 phosphate cuvette test (Phosphate ions react with molybdate and antimony ions in an acidic solution to form an antimonyl phosphomolybdate complex, which is reduced by ascorbic acid to phosphomolybdenum blue). Ortho-PO₄-P was measured by using Hach phosphorus reactive powder pillows (ascorbic acid method adapted from APHA, 1999). NH₄-N was measured by using Hach nitrogen reactive powder pillows (Salicylate method adapted from Reardon et al., 1966). NO₃-N was measured by Hach LCK 339 nitrate cuvette test (Nitrate ions in solutions containing sulphuric and phosphoric acids react with 2.6-dimethylphenol to form 4-nitro-2.6-dimethylphenol). Chl-a was measured by spectrophotometric method (APHA, 1999). DR 5000 UV/VIS spectrophotometer (Hach Lange, Germany) was used in the measurements of TN, TP, ortho-PO₄-P, NH₄-N, NO₃-N, chl-a. Turbidity was measured by using a portable turbidimeter (HF Scientific Micro TPI, Fort Myers, USA). Sampling days of year (DOYs) in the study were given in Table 1.

2.3 Estimations of metabolism components

Metabolism components were calculated by using the following equations (Bales and Nardi, 2007):$$F_r(t)\hspace{0.17em}=\hspace{0.17em}{D}_{\mathrm{a}\mathrm{v}\mathrm{g}}\times \hspace{0.17em}{K}_2(T_w)\hspace{0.17em}\times \hspace{0.17em}Q\hspace{0.17em}\times \hspace{0.17em}{T}_t\times \hspace{0.17em}{C}_f$$(1)

where F_r(t) is the reaeration flux (mg O₂ reach⁻¹ min⁻¹) at time t, D_avg is the reach-averaged DO deficit (mg L⁻¹) as DO_sat (t) (saturated DO concentration) minus DO(t), K₂ (T_w) is the reaeration rate coefficient (min⁻¹) at T_w (°C), Q is stream flow rate (L s⁻¹), T_t is travel time (min), and C_f is the unit conversion factor (1 min = 60 s). In this study, K₂ at T_w = 20 °C was estimated by using the equation (2) (Owens et al., 1964) because it is suitable for characteristics (depth, velocity, etc.) of the sampled reaches.$$K_2=5.35\times V^{0.67}\times D^{-1.85}{T}_w=20°\mathrm{C}\left(0.12\le D\le 3.35m\mathrm{a}\mathrm{n}\mathrm{d}0.03\le V\le 1.52m{s}^{-1}\right)$$(2)

K₂ at any T_w was calculated by using the following equation by Elmore and West (1961):$$K_2(T_w)=K_2(T_w=20°\mathrm{C})\times {1.024}^{(T_w-20)}$$(3)

NEM was calculated for each measurement interval (one minute) as follows:$$\text{NEM}(t)=\left(\left[{\text{DO}}_{\mathit{\text{d}}}(t)-{\text{DO}}_{\mathit{\text{u}}}(t-T_t)\right]\times Q\times C_f\right)-F_r(t)$$(4)

where NEM(t) is the net metabolism flux (mg O₂ reach⁻¹ min⁻¹) at time t, DO_d(t) is downstream DO concentration (mg L⁻¹) at time t, DO_u(t-T_t) is upstream DO concentration (mg L⁻¹) at time t–T_t. The nighttime NEM was assumed to be equal to R_c as no GPP occurs at night. Daytime R_c (R_d) was estimated using the approach proposed by Mulholland et al. (2001) in which mean respiration rates of one-hour pre-dawn and one-hour post-sunset (dusk) periods were calculated and assumed to be respiration rates at the start and end of daylight, respectively. Then, R_d at each measurement interval during daylight was extrapolated between the start and end rates as follows:$$R_d(t)=R_{\mathrm{p}\mathrm{r}\mathrm{e}-\mathrm{d}\mathrm{a}\mathrm{w}\mathrm{n}}+\left[(R_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}-\mathrm{s}\mathrm{u}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{t}}-R_{\mathrm{p}\mathrm{r}\mathrm{e}-\mathrm{d}\mathrm{a}\mathrm{w}\mathrm{n}})/(n_d-1)\times N_i\right]$$(5)where R_d(t) is the daylight respiration flux (mg O₂ reach⁻¹ min⁻¹) at given time t, R_pre-dawn is the mean respiration rate (mg O₂ reach⁻¹ min⁻¹) calculated from individual net metabolism fluxes for the one-hour period before dawn. R_post-sunset is the mean respiration rate (mg O₂ reach⁻¹ min⁻¹) calculated from individual net metabolism fluxes for the one-hour period following sunset. n_d is the total number of measurement intervals during daylight hours including the first and last daylight readings. N_i is the serial number of the daylight measurement interval for calculating R_d(t), with N_i = 1 = 0. A 24-hour total R_c was calculated as g O₂ reach⁻¹ day⁻¹ from the two nighttime periods and the daylight period as the following:$$R_c={\displaystyle \sum _{t=1}^{t=nd}}{R}_d(t)+(S_{N1}+S_{N2})\times \left[\left(24-T_d\right)/\left({}_{\mathrm{}}+T_{N2}\right)\right]$$(6)where S_N1 is the sum of the net metabolism flux for the first night-time period, S_N2 is the sum of the net metabolism flux for the second night-time period, T_d is a duration of the daytime period in hours, T_N1 is a duration of the first night-time period in hours, T_N2 is a duration of the second night-time period in hours. GPP was calculated for each measurement interval (one minute) as follows:$$\mathrm{G}\mathrm{P}\mathrm{P}(t)=\mathrm{N}\mathrm{E}\mathrm{M}(t)-R_c(t)$$(7)

where GPP(t) is gross primary productivity flux (g O₂ reach⁻¹ min⁻¹) at time t, R_c(t) is community respiration flux (g O₂ reach⁻¹ min⁻¹) at time t. GPP, NEM and R_c estimates were converted to g O₂ m⁻² day⁻¹ by using surface areas of reaches. The surface area was calculated by multiplying width and length of the reach. Measured DO data were corrected for the issues of sensor fouling and calibration drift during the deployment using the method by Wagner et al. (2006).

2.4 Statistical analyses

All statistical analyses were performed by using Minitab 17.0 software (Minitab Inc. State College, PA, USA). Pearson's correlation matrix among the metabolism rates and the associated environmental variables in each reach was performed to detect the direction and strength of significant linear relationships. With the parametric tests applied in this study, the following diagnostic assumptions were performed: normality test of residuals of regression models through Q–Q plot and Anderson–Darling statistics, homoscedasticity test through plot of residuals versus fitted values, multicollinearity test through variance inflation factor (according to VIF > 40), and autocorrelation test through Durbin–Watson (DW) statistics. Autocorrelation test was applied according to the following decision rules: no autocorrelation when DW ≈2; a strong positive autocorrelation when DW ≈0; and a strong negative autocorrelation when DW ≈4.

A 17-month-entire dataset was randomly split into training and testing subsets. Based on the parameterization (training) dataset, the best-fit multiple non-linear regression (MNLR) models of NEM and R_c (g O₂ m⁻² day⁻¹) were obtained using the stepwise procedure by which significant predictors were selected from among the environmental variables and their interactions through order 2. Alpha-to-enter and alpha-to-remove values of the stepwise procedure were selected as 0.05. The adjusted coefficient of determination (R ²_adj), and prediction coefficient of determination (R ²_pred) were used to evaluate the goodness-of-fit and predictive power of the MNLR models, respectively. In order to represent spatiotemporal heterogeneity, DOY and naturality gradient (NG) variables (reaches of Abant, Mudurnu, BD and AD) were forced as the predictors into the models. NG variables were selected as categorical predictors in the models. MNLR models were also validated by using a validation (testing) dataset randomly selected from the entire dataset and apart from the parameterization (training) dataset.

Tukey's multiple comparison tests following one-way analysis of variance (ANOVA) were applied to determine significant mean GPP, NEM, R_c differences along the reach disturbance gradient (NG) by using all data between August 2015 and December 2016. Thus, effects of the WTP discharges on metabolism rates could be determined by comparing the reaches. Also, the mean metabolism components were compared in terms of the seasons (winter, spring, summer, autumn). Significance level (α) was selected as 0.05 for all comparisons.

3 Results

3.1 Linear relationships between the environmental variables and the metabolism rates

According to mean monthly environmental variables and metabolism rates between August 2015 and December 2016, AD reach had a maximum R_c rate with highest values of SC, BOD₅, TN, TP, NH₄-N, ortho-PO₄-P, turbidity and flow rate in each month (see Tab. 2 and Appendix Tab. 1). The averages of NEM rates were negative in each reach for all months except August 2015 and April 2016 in BD reach.

Results of Pearson's correlation matrix applied among the environmental variables and the metabolism rates showed that there was a discrepancy in correlations between the environmental variables and the metabolism rates in four reaches. In each reach, these environmental drivers correlated with metabolism rates differently (negative or positive, low or high) (see Appendix Tabs. 3–6). This indicates that the linear relationships could spatiotemporally change. However, it might be expected that high concentrations of BOD₅ (r = 0.349, p = 0.001), NH₄-N (r = 0.478, p < 0.001) and turbidity (r = 0.232, p = 0.037) measured in AD reach have positive and negative correlations to R_c and GPP, respectively (BOD₅: r = −0.154, p = 0.171; NH₄-N: r = −0.228, p = 0.041; turbidity: r = −0.053, p = 0.641) (see Appendix Tab. 6). Because, R_c might be increased by high DO consumption due to the decomposition of organic matter and the oxidation of NH₄-N (nitrification) arising from the discharges of the WTP. On the other hand, GPP might be decreased since high concentrations of turbidity and suspended solids absorbed sunlight and restricted a sufficient solar energy for the photosynthesis. Also, high concentration of ortho-PO₄-P measured in AD reach has a positive correlation to R_c (r = 0.418, p < 0.001). Ortho-PO₄-P and NH₄-N are nutrients for heterotrophic microorganisms as well as autotrophs (Kirchman, 1994). It might be a factor for increasing R_c. Similarly, a positive correlation between T_w and R_c (r = 0.563, p < 0.001) indicates that T_w might stimulate the activity of heterotrophic community in AD reach. This case was seen for GPP in Abant and BD reaches as T_w might stimulate the activity of autotrophs (Abant: r = 0.795, p < 0.001; BD: r = 0.221, p = 0.047).

3.2 MNLR modelling of metabolism rates

Descriptive statistics of the predictors and responses used for the MNLR models were given in Table 2. The MNLR models were derived for R_c and NEM as equations (8) and (10), and they were given in Tables 3 and 4, respectively. The best-fit MNLR model of R_c was parameterized with R ²_adj of 87.14% and R ²_pred of 81.66%. Thus, DOY, NG, P_atm, DOY × T_air, DOY × P_atm, DOY × BOD₅, DOY × SC, RH × BOD₅, T_w × BOD₅, BOD₅ × SC, BOD₅ × pH, BOD₅ × NH₄-N, BOD₅ × ortho-PO₄-P, Turbidity × SC, Turbidity × pH, Turbidity × Chl-a, Turbidity × NO₃-N, pH × NH₄-N predictors (independent variables and interactions) could explain 81.66% of R_c.

According to F-value, the most associated with R_c (response) and the significant predictor is the interaction of BOD₅ × SC. In terms of VIF values, DOY, NG (except Mudurnu reach), DOY × P_atm, DOY × BOD₅, DOY × SC, T_w × BOD₅, BOD₅ × SC, BOD₅ × pH, BOD₅ × NH₄-N, BOD₅ × ortho-PO₄-P, Turbidity × SC, Turbidity × pH, pH × NH₄-N predictors have multicollinearity issues (VIF > 40) (see Tab. 3). Multicollinearity problem is not desirable in terms of the model confidence. However, there is no exact assumption regarding the multicollinearity for MNLR models. DW value (DW = 1.117) of the model indicates that there is no autocorrelation issue.

The R_c model could predict the validation data in the ratio of 66.6% by comparing calculated (measured) and predicted R_c values (R ² = 66.6%, S = 25.98, p = 0.001, n=172: total number of explanatory variables of R_c MNLR model and measured R_c data for validation).$$\begin{array}{l}{R}_c=6642-11.43\cdot \mathrm{D}\mathrm{O}\mathrm{Y}-10.02\cdot P_{\mathrm{a}\mathrm{t}\mathrm{m}}+0\cdot \mathrm{A}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{t}+384.4\cdot BD-975\cdot AD\\ +302\cdot \mathrm{M}\mathrm{u}\mathrm{d}\mathrm{u}\mathrm{r}\mathrm{n}\mathrm{u}-0.01082\cdot \mathrm{D}\mathrm{O}\mathrm{Y}\times T_{\mathrm{a}\mathrm{i}\mathrm{r}}+0.01793\cdot \mathrm{D}\mathrm{O}\mathrm{Y}\times P_{\mathrm{a}\mathrm{t}\mathrm{m}}\\ +0.05699\cdot \mathrm{D}\mathrm{O}\mathrm{Y}\times \mathrm{B}\mathrm{O}{\mathrm{D}}_5-0.001648\cdot \mathrm{D}\mathrm{O}\mathrm{Y}\times \mathrm{S}\mathrm{C}-0.03896\cdot \mathrm{R}\mathrm{H}\times \mathrm{B}\mathrm{O}{\mathrm{D}}_5\\ +0.2618\cdot T_w\times \mathrm{B}\mathrm{O}{\mathrm{D}}_5-0.03089\cdot \mathrm{B}\mathrm{O}{\mathrm{D}}_5\times \mathrm{S}\mathrm{C}+6.248\cdot \mathrm{B}\mathrm{O}{\mathrm{D}}_5\times \mathrm{p}\mathrm{H}\\ -0.8503\cdot \mathrm{B}\mathrm{O}{\mathrm{D}}_5\times \mathrm{N}{\mathrm{H}}_4\text{−}\mathrm{N}-10.320\cdot \mathrm{B}\mathrm{O}{\mathrm{D}}_5\times \mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}-\mathrm{P}{\mathrm{O}}_4\text{−}\mathrm{P}\\ +0.00638\cdot \mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{y}\times SC-1.130\cdot \mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{y}\times \mathrm{p}\mathrm{H}\\ +0.01512\cdot \mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{y}\times \mathrm{C}\mathrm{h}\mathrm{l}\text{−}a+0.661\cdot \mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{y}\times \mathrm{N}{\mathrm{O}}_3\text{−}\mathrm{N}\\ +8.092\cdot \mathrm{p}\mathrm{H}\times \mathrm{N}{\mathrm{H}}_4\text{−}\mathrm{N}\end{array}$$(8)

The validation equation of R_c is as follows:$${\mathit{\text{R}}}_{\text{cpred}}=2.98+1.413R_{\text{cmeas}}$$(9)

The best-fit MNLR model of NEM was parameterized with R ²_adj of 84.6% and R ²_pred of 80.04%. DOY, NG, DOY × T_air, DOY × NH₄-N, T_air × TP, BOD₅ × SC, BOD₅ × pH, BOD₅ × NH₄-N, BOD₅ × ortho-PO₄-P, TN × Turbidity, TP × Turbidity, Turbidity × pH, SC × NH₄-N, NH₄-N × NO₃-N predictors could explain 80.04% of NEM. The most associated with NEM (response) and the significant predictor is the interaction of DOY × NH₄-N.

According to VIF values, DOY × NH₄-N, BOD₅ × SC, BOD₅ × pH, BOD₅ × NH₄-N, BOD₅ × ortho-PO₄-P, TN × Turbidity, TP × Turbidity, SC × NH₄-N predictors have multicollinearity issues (VIF > 40) (see Tab. 4). DW value (DW = 1.48) of the model indicates that there is no autocorrelation issue.

The NEM model could predict the validation data in the ratio of 74.9% by comparing calculated (measured) and predicted NEM values (R ² = 74.9%, S = 24.77, p ≤ 0.001, n = 184: total number of explanatory variables of NEM MNLR model and measured NEM data for validation).$$\begin{array}{l}\text{NEM}=-12\text{.7}-0\text{.0905}\cdot \text{DOY+0}\cdot \text{Abant}-6\text{.2}\cdot \text{BD}-1141\text{.5}\cdot \text{AD}-5\text{.7}\cdot \text{Mudurnu}\\ -0\text{.01349}\cdot \text{DOY}\times {\text{T}}_{\text{air}}\text{+}0\text{.2033}\cdot \text{DOY}\times {\text{NH}}_4\text{−}\text{N + 9.25}\cdot {\text{T}}_{\text{air}}\times \text{TP}\\ -0\text{.03072}\cdot {\text{BOD}}_5\times \text{SC +}6\text{.428}\cdot {\text{BOD}}_5\times \text{pH}-0\text{.3009}\cdot {\text{BOD}}_5\times {\text{NH}}_4\text{−N}\\ -2\text{.773}\cdot {\text{BOD}}_5\times {\text{ortho-PO}}_4\text{−P}\text{+ 0.8958}\cdot \text{TN}\times \text{Turbidity}\\ -8\text{.888}\cdot \text{TP}\times \text{Turbidity}-0\text{.3405}\cdot \mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{y}\times \text{pH}\\ -0\text{.03970}\cdot \text{SC}\times {\text{NH}}_4\text{−}\text{N + 10.36}\cdot {\text{NH}}_4\text{−}\text{N}\times {\text{NO}}_3\text{−}\text{N}\end{array}$$(10)

The validation equation of NEM is as follows:$$\mathrm{N}\mathrm{E}{\mathrm{M}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=\mathrm{-14.18}+0.9057-\mathrm{N}\mathrm{E}{\mathrm{M}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\text{.}$$(11)

Although high R²_adj and R²_pred values (87.14% and 81.66%, 84.6% and 80.04% respectively) of R_c and NEM models, low values (10.39%, 7.67%) were obtained for GPP. Instead of using GPP model, GPP estimations were obtained by differences among NEM and R_c predicted from the MNLR models (see Eq. (7)). Calculated (measured) and predicted GPP values were compared for validation with the following equation (R ² = 54.8%, S = 10.63, p = 0.152, n = 15).$$\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=\mathrm{-2.447}+0.675\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$$(12)

3.3 One-way ANOVA Tukey's multiple comparisons of GPP, NEM and R_c along the reach disturbance gradient (NG)

Results of the comparisons using all of 17-month-metabolism data between August 2015–December 2016 and mean metabolism rates were given in Table 5 and Figures 2 and 3, respectively. The results showed that there was no significant difference statistically among the group means of GPP in all reaches (GPP group means are A for all reaches). However, according to values of the means, it can be said that discharges of the WTP reduce GPP (GPP before discharges: 15.6 g O₂ m⁻² day⁻¹, GPP after discharges: 9.1 g O₂ m⁻² day⁻¹). There was no significant difference among the group means of NEM and R_c in reaches of Abant, Mudurnu, and BD, except AD reach (NEM and R_c group means are A for Abant, Mudurnu and BD reaches; B for AD reach). R_c increased after the discharges of the WTP (R_c before discharge: −30.6 g O₂ m⁻² day⁻¹, R_c after discharge: −130.9 g O₂ m⁻² day⁻¹).

Results of the comparisons according to the seasons and mean seasonal metabolism rates were given in Table 6 and Figures 4 and 5, respectively. Comparisons of the mean metabolism components in terms of the seasons enabled the revelation of the following differences (exceptions): In the winter, there was no significant difference statistically among the group means of NEM and R_c in all reaches including AD reach (NEM and R_c group means are A for all reaches), in contrast to the results obtained using all data (see Tables 5 and 6); Although there was not any significant difference among the group means of GPP in BD and AD reaches (GPP group means are A for all reaches), the average of GPP increased after the discharges (GPP before discharge: 2.385 g O₂ m⁻² day⁻¹, GPP after discharge: 8.72 g O₂ m⁻² day⁻¹) in the winter; In the summer, there was a significant difference among the group means of R_c in BD and Abant–Mudurnu reaches (see Tab. 6). In the autumn, there was a significant difference among the group means of GPP in BD and AD reaches, and GPP decreased after the discharges (see Tab. 6).

In general, the results showed that Bolu City WTP discharges affected the metabolism of Büyüksu Stream: R_c increased while GPP decreased (although there was no significant GPP difference between BD and AD, statistically) after the discharges.

Fig. 2 Mean metabolism rates with standard deviations of Abant and Mudurnu reaches according to all data.

Fig. 3 Mean metabolism rates with standard deviations of BD and AD reaches according to all data.

Fig. 4 Mean metabolism rates with standard deviations of Abant and Mudurnu reaches according to seasons.

Fig. 5 Mean metabolism rates with standard deviations of BD and AD reaches according to seasons.

4 Discussion

In this study, it was demonstrated whether the statistical models could be derived to predict the metabolism rates in the real environmental conditions (not laboratory conditions). The validation results of the models for NEM and R_c (except GPP) were obtained as high (74.9% and 66.6%, respectively) according to the real conditions. The metabolism rates under different values of the environmental variables in the four reaches can be roughly estimated by using these models.

The study demonstrated that linear relationships among the measured environmental variables and the metabolism rates varied spatiotemporally. It can be said that these environmental variables do not individually explain the metabolism rates since they correlated with the metabolism rates differently in each reach (negative or positive, low or high). Therefore, interactions of the environmental variables were included in the MNLR models. Furthermore, in order to represent this spatiotemporal heterogeneity, DOY and reaches of Abant, Mudurnu, BD and AD were used as explanatory variables (predictors) of the models. In this regard, second-order interactions among the different variables were considered in the MNLR models in addition to individual variables.

Low R ²_adj and R ²_pred values (10.39%, 7.67%) of GPP model indicated that the environmental variables used as predictors in R_c and NEM models were inadequate to explain GPP statistically. In other words, there is a lot of unexplained variation. Therefore, in addition to the environmental variables used in the study, other variables associated with GPP such as photosynthetically active radiation (PAR), shading percentage (%) (Burrell et al., 2014), etc. should be considered for modelling GPP. Similarly, other variables associated with R_c and NEM such as sediment oxygen demand (SOD) may be used as an individual predictor or may be interacted with the environmental variables in the models to obtain better results.

In general, mean NEM rates were negative in each reach. This indicated that biomass in the reaches decomposed, and generally heterotrophic environment (conditions) was dominant in the four reaches throughout the year. It was found that there was no significant difference between the group means of NEM and R_c in reaches of Abant, Mudurnu, and BD, except AD reach. This indicates that the discharges of the WTP effect NEM and R_c, and according to values of R_c means, the discharges enhance R_c. R_c might be increased since decomposition of high concentrations of organic matter, oxidation of high concentration of NH₄-N (nitrification), and heterotrophic microorganisms in discharges of the WTP. Similarly, in the related literature, several researchers such as Gücker et al. (2006), Chen (2013) and Chesworth (2016) found that WTP discharges increased R_c rates. It can be seen that the similarities and differences of mean R_c and GPP among this study and previous studies in Table 7. As it is seen, the mean R_c of this study (−130.9 g O₂ m⁻² day⁻¹) is greater than the other studies. Dense untreated wastewater discharges might cause this difference (mean BOD₅ = 79.6 mg L⁻¹, mean NH₄-N = 15.7 mg L⁻¹ in the water after discharges). Gücker et al. (2006) reported mean BOD₅ = 7–9 mg L⁻¹ and mean NH₄-N = 0.13–0.16 mg L⁻¹ in the downstreams of the discharges (a second greatest mean R_c in Tab. 7). In the summer, high T_air and T_w values might accelerate the growth and activity of the heterotrophic community in the water depending on the sites (In the summer, mean daily T_air and T_w: 19.0–21.3 °C, 20.2–15.4 °C, 21.1–20.0 °C for Abant, Mudurnu and BD-AD reaches, respectively). Thus, the means of R_c might differ significantly in BD and Abant − Mudurnu reaches (see Tab. 6).

It was found that there is a decrease (not significant statistically) in mean GPP after the WTP discharges. This decrease might be because the photosynthesis could not occur while the sunlight was absorbed by the untreated discharges which contain high concentrations of turbidity and suspended solids as well as treated discharges of the WTP. This reason is supported by Kicklighter (1987) who indicated that GPP was decreased because of the dark colour of the wastewater. In contrast, Gücker et al. (2006) found that discharges of treated wastewater containing nutrients (NH₄-N, NO₃-N, ortho-PO₄-P) stimulated GPP rates. But in this study, although there were high concentrations of nutrients in the discharges, the influence of the nutrients on the growth rate of photosynthesizing organisms might be restrained by the untreated wastewaters (due to light deficiency as noted above). Aristi et al. (2015) agreed with this assumption by demonstrating that the increase in the nutrient loads could not affect GPP in the reaches where sunlight was restricted. In the winter, mean GPP increased after the discharges. This might be since a negative effect of the WTP untreated discharges on the transmittance of sunlight was decreased by the dilution through high stream flow rates with precipitations in the winter. High ranges of the metabolism rates in BD and AD were measured due to extreme flow rates with the precipitations in the winter and spring (see Appendix Tab. 1, Fig. 5). Because, the metabolism calculations were also based on stream hydraulics. Moreover, unmeasured contaminants by runoff might disturb the metabolism rates. These high variabilities might result in that some statistical differences between the reaches were not significant.

5 Conclusion

With this study, the effects of the WTP discharges on metabolism rates were investigated. The results demostrated that (1) the effects on R_c were greater than GPP, (2) R_c was increased by the discharges, (3) GPP was reduced (except the winter) although there was no significant GPP difference between BD and AD, statistically. Moreover, the linear relationships among the environmental variables measured and the metabolism rates might vary spatiotemporally (in each reach). NEM and R_c rates as a function of measured environmental variables could be spatiotemporally modelled by including DOY and NG (R ² of 74.9% and 66.6% of model validation, respectively).

Stream hydraulics (stream flow rate, velocity, depth, width) were effective in estimations of the metabolism components. In this study, one-min-interval linear interpolations based on time were applied for Q, V, and D measured at the start and end of the deployment of the DO loggers instead of average (two measurements) for more accurate estimations of GPP, NEM, and R_c. However, it's recommended that simultaneously DO − T_w and Q–V–D measurements should be performed at each measurement interval as much as possible instead of average or linear interpolations. Also, experiences obtained from in-situ measurements in the study suggest that DO probes (data loggers) should not be left in the stream for more than two days because fouling could affect DO readings. In future studies, uncertainty analysis (such as Monte Carlo) can be carried out for MNLR models derived in this study.

WTPs are widespread around the world. Findings of the study prove that the WTP discharges increase the mean R_c. It causes decreases in DO concentration in the water. Thus, aquatic ecosystem and life can be endangered. Even though it was found that the discharges have no impact on the mean GPP as much as R_c in the study, GPP may be increased by the nutrients where turbidity is low and sunlight is sufficient. This may result in eutrophication. In this regard, metabolism rates can be used to show water quality of the streams, and they should be considered in the discharge standards of the WTPs for the receiving water bodies (river, lake, etc.).

Supplementary Material

Supplementary Tables S1 to S6. Access here

Acknowledgments

This study was financially supported by the Scientific Research Projects Unit of Bolu Abant Izzet Baysal University (Grant No: BAP–2015.09.02.851). We are grateful to H. Fidan and E. Özalp for their help with field campaigns.

References

All Tables

All Figures

	Fig. 1 Locations of the four reaches in the basin of Büyüksu Stream.
In the text

	Fig. 2 Mean metabolism rates with standard deviations of Abant and Mudurnu reaches according to all data.
In the text

	Fig. 3 Mean metabolism rates with standard deviations of BD and AD reaches according to all data.
In the text

	Fig. 4 Mean metabolism rates with standard deviations of Abant and Mudurnu reaches according to seasons.
In the text

	Fig. 5 Mean metabolism rates with standard deviations of BD and AD reaches according to seasons.
In the text