*Article* **A River Temperature Model to Assist Managers in Identifying Thermal Pollution Causes and Solutions**

#### **Reza Abdi \* and Theodore Endreny**

Department of Environmental Resources Engineering, SUNY ESF, Syracuse, NY 13210, USA; te@esf.edu **\*** Correspondence: reabdi@syr.edu; Tel.: +1-315-470-6565

Received: 27 April 2019; Accepted: 15 May 2019; Published: 22 May 2019

**Abstract:** Thermal pollution of rivers degrades water quality and ecosystem health, and cities can protect rivers by decreasing warmer impervious surface stormwater inflows and increasing cooler subsurface inflows and shading from riparian vegetation. This study develops the mechanistic i-Tree Cool River Model and tests if it can be used to identify likely causes and mitigation of thermal pollution. The model represents the impacts of external loads including solar radiation in the absence of riparian shade, multiple lateral storm sewer inflows, tributaries draining reservoirs, groundwater flow, and hyporheic exchange flow in dry weather steady flows and wet weather unsteady flows. The i-Tree Cool River Model estimates the shading effects of the riparian vegetation and other features as a function of heights and distances as well as solar geometry. The model was tested along 1500 m of a New York mountain river with a riparian forest and urban areas during 30 h with two summer storm events in 2007. The simulations were sensitive to the inflows of storm sewers, subsurface inflows, as well as riparian shading, and upstream boundary temperature inflows for steady and unsteady conditions. The model simulated hourly river temperature with an R<sup>2</sup> of 0.98; when shading was removed from the simulation the R<sup>2</sup> decreased 0.88, indicating the importance of riparian shading in river thermal modeling. When stormwater inflows were removed from the simulation, the R2 decreased from 0.98 to 0.92, and when subsurface inflows were removed, the R2 decreased to 0.94. The simulation of thermal loading is important to manage against pollution of rivers.

**Keywords:** River thermal pollution; Mechanistic model; Urban hydrology; Riparian shading; Heat balance

#### **1. Introduction**

Excessive river temperatures are detrimental to water quality and ecosystem health [1]. River warming can result from increased inflows of warm point and non-point source discharges, decreased inflows of cool sub-surface waters, removal of riparian shade, increased air temperatures, and changes in channel substrate and depth that increase absorption, conduction, and convection in heat transfer [2–4]. River thermal pollution is often associated with discharges of coolant water used by industry, but it is also associated with land-use change, including urbanization, river impoundment, channel management, and regulation [5,6]. River temperature is a critical water quality parameter for riverine systems, that affects the saturation of dissolved oxygen [7,8], kinetic reactions and resulting pollutant concentrations [9,10], and fish distribution, metabolism, growth, reproduction, and mortality [11,12]. Urbanization can elevate river temperatures through changes in riparian land cover which affects shade on the water surface, through river morphology which affects water depth, surface area, and velocity, and through flow connectivity with groundwater, stormwater, and other point and non-point source inflows [13–18]. When precipitation strikes hot impervious surfaces of urban areas, this generates warmer stormwater relative to the temperature of river water [19–21].

River temperature management and mitigation of thermal pollution are best planned with simulation models that enable scenario evaluation, to explore relationships between river temperature response and drivers that vary in space and time [22,23]. Caissie [24] provides a review of research into the spatial and temporal drivers of river temperature, and the evolution of modeling approaches, including statistical models and cause-effect deterministic models, which we call mechanistic models. An illustration of field observed spatial and temporal variation is provided by Webb et al. [25] who monitored 11 reaches in south-west England, through July 1992 to February 1993, noting a correspondence between variations in hourly river temperature and variation in discharge and drainage area (180 km2 to 0.4 km2), channel surface area (17 m2 to 439 m2), depth (0.1 m to 0.5 m), slope (0.002 to 0.05), orientation (north-south to east-west), and riparian cover (pasture to dense woodland). They found drainage basin land cover, especially riparian vegetation, overwhelmed other drivers of temperature. In two separate studies of more than16 rivers in the Washington, DC area, it was determined that runoff from impervious land uses entering rivers through urban storm sewers was the major thermal stressor, causing rapid (< 3 h) surges in temperature greater than 3 ◦C [1,26].

A comprehensive mechanistic model developed for river managers is Heat Source [27], which allows the user to specify local climatology, hydrology, morphology, and land use into Microsoft Excel and ESRI Arc View software, to simulate spatial and temporal variation in river temperature as a function of shortwave and longwave radiation, sensible and latent heat, riverbed conduction, and inflows from tributaries, groundwater, and hyporheic exchange. Data and algorithm limitations can constrain utilization of input-intensive mechanistic models, but these limitations can sometimes be overcome with innovations. Yearsley [28] developed a semi-Lagrangian scheme to advect river heat within a large river network when channel morphology data needed by Heat Source [27] and similar models (e.g., HSPF, CE-QUAL-W2, and QUAL2K) were not measured. To better capture abrupt changes in velocity and riparian shading, Crispell [29] created a retention time alternative to the advection-dispersion routing algorithm used in Heat Source [27], which maintains numerical stability at very fine spatial but coarse temporal discretization. For cases when observed boundary condition data are not available, Sun et al. [18] modified the DHSVM–RBM mechanistic model of river temperature to use Mohseni et al. [30] non-linear regression between weekly air temperature and river temperature to generate the upstream river temperature time series needed as a boundary condition.

The portability and accessibility of river temperature models are significant limitations for users interested in river management, pollution mitigation and restoration scenarios. The Heat Source model provides a balance between scenario simulation options and model parsimoniousness that made it our choice as the base code for developing a free, open-source, lower-complexity river temperature model useful in river pollution mitigation and restoration. The complexity of HSPF, CE-QUAL-W2, and QUAL2K is high, each containing many non-temperature routines, and CE-QUAL-W2 representing a 2-dimensional (2D) domain. These three models do not simulate ecological processes important in scenario analysis, including hyporheic exchange and temporal variation in the riparian shade; HSPF does have a pre-processor to provide temporal variation in riparian shade, which we use in our code development [14]. Glose et al. [31] noted that the major limitation of Heat Source [27] is lack of automation in making multiple simulations for parameter calibration and sensitivity analysis, given it is written in the Visual Basic for Applications language within Microsoft Excel. Glose et al. [31] addressed this limitation by using Matlab, a well-supported scientific programming language, to create the steady state model HFLUX, which represents many of the mechanistic processes in Heat Source. The HFLUX model [31] does not include the shade factor estimation and unsteady state routing algorithms of Heat Source [27], which are important in cases of temporal variation in shading and storm flow dynamics. Unfortunately, neither Heat Source [27] nor HFLUX [31] can be compiled into an executable, and therefore cannot be deployed outside of the VBA or Matlab environment.

This study created the i-Tree Cool River Model to address limitations of the Heat Source [27] and HFLUX [31] models and advance mechanistic model simulation of river management, pollution mitigation, and restoration scenarios in a parsimonious manner. The i-Tree Cool River Model is designed to allow for flexible shading factor algorithms, steady and unsteady flow, as well as other heat and mass transfer processes. The i-Tree Cool River Model is an open-source tool written in C++, and its package contains the routines and an executable file for running the code, which can be downloaded from http://www.itreetools.org/research\_suite/coolriver. The model executable is called at the command line along with a configuration extensible markup language (XML) file, which includes the required initial information. The i-Tree Cool River Model C++ algorithms can be edited and recompiled with Visual Studio 2017 Community Edition or later, which is freeware. The simulation output includes the simulated river temperature and the heat fluxes.

The objectives of this paper are to present the theory of the i-Tree Cool River Model, to apply the model in a case study with unsteady stormwater inflows, and to evaluate the importance of the heat and mass transfer processes. To that end, following the model development, the manuscript provides a model testing to address the application of the model. The science questions are: When analyzing sources of thermal pollution, and possible mitigation scenarios, what is the relative contribution of (a) storm sewer inflows, (b) subsurface inflows of groundwater and hyporheic exchange, (c) riparian shading and weather, on the accuracy of simulated river temperature?

#### **2. Methods**

#### *2.1. Heat Flux Formulation*

The i-Tree Cool River Model simulates an advection-dispersion equation with inflows and heat fluxes following Martin et al. [32]

$$\frac{\partial T\_w}{\partial t} = -lI \frac{\partial T\_w}{\partial \mathbf{x}} + D\_L \frac{\partial^2 T\_w}{\partial \mathbf{x}^2} + R\_{l\bar{l}} + R\_{\bar{l}} \tag{1}$$

where *Tw* is the cross-sectional averaged river temperature (◦C), *t* is time (s), *U* is the reach average flow velocity (m/s), *x* is river distance (m), *DL* is the dispersion coefficient (m2/s), *Rh* is the heat flux reaction term, also known as heat transfer [3,27], and *Ri* is the reaction term of the external inflows. When *Ri* is combined with the advection and dispersion terms in Equation (1), they are collectively referred to as mass transfer [3,27]. The *Rh* and *Ri* are defined as

$$R\_h = \frac{\Phi\_{net}}{\rho \mathbb{C}\_p y} \tag{2}$$

$$R\_{\bar{i}} = \frac{Q\_W T\_W + Q\_{GW} T\_{GW} + Q\_{Hyp} T\_{Hyp} + Q\_{SS} T\_{SS}}{Q\_{\bar{i}} + Q\_{GW} + Q\_{Hyp} + Q\_{SS}} - T\_{W,t-1} \tag{3}$$

where Φ*net* is the net exchange of thermal energy (W/m2), ρ is the water density (kg/m3), *Cp* is the specific heat capacity of water (J/kg ◦C), *y* is the average water column depth (m), *Q* is discharge (m3/s), *T* is water temperature (◦C), and subscripts *<sup>W</sup>* is the river flow, *GW* is groundwater flow, *Hyp* is hyporheic exchange, and *SS* is stormwater inflow, *t* − 1 is prior time step. River velocities, dispersion, and inflows are calculated using standard methods, described in Supplementary Materials Section S1. The subsurface inflows distinguish between hyporheic and groundwater inflows due to their different environmental processes and use a separate mathematical formulation for each term. For surface inflows, users can include tributaries in place of storm sewers, and assign an unlimited number of surface inflows for each cross section.

The net exchange of thermal energy is defined as in Boyd et al. [27] as

$$
\Phi\_{\text{nrt}} = \Phi\_{\text{longuure}} + \Phi\_{\text{shortuure}} + \Phi\_{\text{latent}} + \Phi\_{\text{scusible}} + \Phi\_{\text{sediment}} \tag{4}
$$

where the Φ is the heat flux (W/m2), and subscripts *net* is the net heat flux at the water surface, *longwave* is the longwave radiation flux at the water surface, *shortwave* is the shortwave radiation at the water surface, *latent* is the latent heat flux from evaporation, *sensible* is the sensible heat flux representing the convective thermal flux from the water surface, and *sediment* is the bed sediment heat flux representing conduction forcing at the water column interface.

The longwave radiation flux in Equation (4) is composed of positive downward fluxes from the atmosphere and land cover over the water surface, and a negative upward flux from the waterbody to the air, following the approach of Boyd et al. [27]

$$
\Phi\_{\text{longmuure}} = \Phi\_{\text{longmuure}}^{\text{atmospheric}} + \Phi\_{\text{longmuure}}^{\text{land}\text{covev}} + \Phi\_{\text{longmuure}}^{\text{huck}} \tag{5}
$$

where Φ*atmospheric longwave* is the atmospheric flux (W/m2), <sup>Φ</sup>*land*cov*er longwave* is the land cover flux (W/m2), and <sup>Φ</sup>*back longwave* is the back-to-air flux (W/m2). Atmospheric longwave radiation is a function of air temperature and exposure from the river surface to the atmosphere, called the view-to-sky factor (*f*), calculated using Boyd et al. [27]

$$
\Phi\_{\text{longurve}}^{\text{atmosplex}} = 0.96 \varepsilon\_{\text{atm}} \sigma (T\_{\text{air}} + 273.2)^4 \min(f\_1, f\_2, f\_3) \tag{6}
$$

where *Tair* is air temperature (◦C), the ε*atm* is the emissivity of the atmosphere (0 to 1), σ is the Stefan-Boltzmann constant (5.6696 <sup>×</sup> 10−8, W/m2K4), and min(*f* 1, *f* 2, *f* 3) is the minimum of the three view-to-sky factors (0 to 1), where *f* <sup>1</sup> represents building effects, *f* <sup>2</sup> represents vegetation effects, and *f* <sup>3</sup> represents topographic effects (Figure 1). The emissivity of the atmosphere ε*atm* is calculated using [33]

$$\kappa\_{atm} = 1.72 (\frac{0.1 \text{g}\_{\text{el}}}{T\_{air} + 273.2})^{\frac{1}{T}} (1 + 0.22 \text{C}\_{\text{L}}^{-2}) \tag{7}$$

where *ea* is the actual vapor pressure (mbar), and *CL* is the cloudiness, which ranges from 0 for a clear sky to 1 for full cloud cover [34].

**Figure 1.** Shading of the river surface. A cross-sectional view, in which BSA is the building shading angle, VSA is the vegetation shading angle, and TSA is the topographic shading angle. hbuilding, htree, and hbank are building, vegetation, and bank heights respectively. Dbuilding, Dcanopy, and Dbank are building to the bank, canopy to the bank, and bank.

The view-to-sky factors value of 1 indicates a full unobstructed sky view [27,35,36]. The general sky-view-factor formula for *fi* is computed for each cross-section based on Chen et al. [14]

$$f\_i = 1 - \frac{2}{\pi} S A\_i \tag{8}$$

where *i* indicates the object at that cross-section, where 1 = building, 2 = vegetation, or 3 = topography; and *SA* is the shade angle (radians), computed as and *hc* is the combined height of the objects above the water (e.g., if a tree is set on a hill, *hc* = *htree* + *hbank*), and max (*Di*) is the maximum distance from all objects at that cross-section to the edge of the water.

$$SA\_i = \tan^{-1} \left( \frac{h\_c}{\max(D\_{i, 1-3})} \right) \tag{9}$$

The land cover longwave radiation in Equation (5) also uses the view-to-sky factors. The land cover radiation represents the land cover, e.g., vegetation such as trees' influence on water temperature, and the model by default sets land cover temperature equal to atmospheric temperature, following the approach of Boyd et al. [27]

$$\Phi\_{\text{longmur}}^{\text{landcover}} = 0.96(1 - \min(f\_1, f\_2, f\_3))0.96\sigma(T\_{\text{air}} + 273.2)^4 \tag{10}$$

The waterbody to air radiation term in Equation (5) is a function of water temperature, representing heat flux emitted from the water surface, following the approach of Boyd et al. [27]

$$
\Phi\_{\text{longware}}^{\text{back}} = -0.96 \sigma (T\_w + 273.2)^4 \tag{11}
$$

where the *TW* is the river temperature (◦C).

See the Supplementary Materials, Section S2 for the methods used to find the remaining right-hand side terms in Equation (4), which are short wave radiation, latent heat flux, sensible heat flux, and bed sediment heat flux. Table S1 lists the 10 input files required by i-Tree Cool River, and names and describes the parameters in each of the files.

#### *2.2. Study Area and Model Inputs*

The i-Tree Cool River Model's accuracy in representing thermal loading was tested in unsteady state (i.e., wet weather) using unpublished data from 11 to 12 June 2007 for a 1500 m reach of Sawmill Creek, in Tannersville, New York (42.1955 N, 74.1339 W, WGS84). Sawmill Creek is a second-order mountainous river with varying watershed land use, starting in forests and transitioning to urban land. At the end of the Sawmill Creek study reach, the time of travel was approximately 30 min and the upstream watershed area is 8.16 km2, which includes a nested urban watershed of 1.8 km2 draining to the river in storm sewers (Figure 2a,b). Sawmill Creek flow at the upstream boundary was estimated using stage-discharge relations, monitoring stage with pressure transducers (manufactured by Global Water Instruments) at the upstream and downstream stations (0 m and 1500 m respectively) and in storm sewers. Stage was converted to discharge using the Manning equation, with stage converted to channel area and hydraulic radius using geometry relations, and the Manning roughness coefficients estimated from pebble counts Wolman [37] at each cross section by Crispell et al. [29]. We installed compound weir plates in the sewers and used the weir plate manufacturer equations to convert stage to discharge for the storm sewer inflows. Observed storm flows were corroborated with simulated flows by the i-Tree Hydro model (Yang, et al. [38]), calibrated to match the estimated baseflow.

Rainfall occurred twice during 12 June 2007 the first time with a 2 h duration totaling 3.3 mm and the second time witha3h duration totaling 8.4 mm. The storm sewer inflows were active during dry and wet weather, in dry weather due to illicit connections draining buildings, and in wet weather due to storm runoff. Crispell et al. [29] monitored river temperatures and storm sewer drainage temperatures at 12 river stations in the Sawmill Creek reach and the two inflow locations at 10-minute intervals using ibutton temperature data loggers, which were used to set boundary conditions. The temperature monitoring ibuttons in Sawmill Creek reach were strategically placed and considered representative of the reach temperature, capturing the influence of stormwater inflows after they had distance to mix with the channel water [29]. In the upstream section of the reach, between the cross section at station

0 m and a station at 600 m, the observed average river temperature increased by 0.03 ◦C per 100 m (3.3%), and between the cross section at station 900 m and a station 1500 m, temperature increased by 0.008 ◦C per 100 m (1%). In the middle section of the reach, between the cross section at station 600 m and station 900 m, the temperature increased by 0.1 ◦C per 100 m, three times the rate of the upstream reach, an increase attributed to the warming effect of the Tannersville's storm sewer inflow (Table S2).

The i-Tree Cool River was also tested in steady state mode, e.g., no rainfall events, for a 475 m reach of Meadowbrook Creek (43.0306 N, 76.0680 W, WGS84), a first order and urbanized river in the city of Syracuse, New York (Figure 3a,b). Flow at the upstream boundary, cross section survey data, and river temperature at 30 monitoring locations at 5-minute intervals were provided by Glose et al. [31], who used these data to develop HFlux. We simulated the 5-day period of 13–19 June 2012.

**Figure 2.** (**a**) New York State with the Sawmill Creek study area denoted by the star. (**b**) Monitoring stations and reach distances along Sawmill Creek.

The i-Tree Cool River Model simulated Sawmill Creek using input data from multiple sources. Specification of hourly weather data, including air temperature *Tair*, dew point temperature *Tdew*, shortwave radiation *Sin*, fraction cloudiness *C*, and wind speed *Uwind* were obtained from the National Solar Radiation Data Base (NSRDB) station ID#1227776, located 23 km from the study site. The NSRDB provides satellite estimated surface radiation at 30 min intervals. The shading factor, *SF*, in Equation (S2) and view-to-sky coefficients, *f* in Equation (S4) were estimated at 1 m intervals along Sawmill Creek using the TTools algorithm from observations of riparian vegetation and aerial images of the study area [29]. The river base width and bank slope were obtained from field surveys at each cross section [29], which defined the irregular pattern of river widening and narrowing. The simulated Sawmill Creek reach was delineated into 15 segments considering the locations where the temperature was observed, with segment lengths no greater than 100 m, which resulted in a simulation timestep

of 0.5 seconds to satisfy the i-Tree Cool River Model stability criteria. The simulations represented the observed conditions, as well as alternative scenarios to determine model sensitivity to shading, subsurface inflow, and the calculated upstream boundary condition, which are sometimes difficult to obtain. Our calculated upstream boundary condition was derived with Mohseni et al. [30] Equation (12), a non-linear regression between air temperature and river temperature.

$$T\_W = \mu + \frac{\alpha - \mu}{1 + e^{\gamma(\beta - T\_{air})}} \tag{12}$$

In the equation, the coefficient α is the estimated maximum stream temperature, γ is a measure of the steepest slope of the function, and β represents the air temperature at the inflection point.

**Figure 3.** (**a**) New York State with theMeadowbrook Creek study area denoted by the star. (**b**)Monitoring stations and reach distances along Meadowbrook Creek.

Our observed upstream boundary condition was obtained from ibutton thermistor measurements. We analyzed the simulated and observed river temperatures for each of the cross sections averaged with respect to time to obtain a 30-hour average at each of the 12 cross sections. Simulations were written hourly for each cross section, which can be written at any timesteps for each meter of the river. We ran this simulation using Equations (5) to (11) for longwave radiation, Equation (S2) for shortwave radiation, Equation (S11) for latent heat, Equation (S12) for sensible heat, and Equation (S15) for the sediment heat.

#### **3. Results**

#### *3.1. Model Evaluation*

A scatterplot of the 30-hours of simulated and observed river temperatures for each of the 12 cross sections along the 1500 m of Sawmill Creek reach provides insights on the relative goodness of fit for each cross section and associated drivers of i-Tree Cool River Model accuracy (Figure 4). At cross section 1, along the upstream boundary, as expected, the observed boundary condition, resulted in a model fit with an *R*<sup>2</sup> of 1.0. Downstream, the fit degraded. Initial conditions of 14.7 ◦C for all cross sections caused the largest deviations between simulated and observed temperatures for most scatterplots. The model underestimated observed temperatures at cross sections 4 to 12 by approximately 1 ◦C, while upstream cross section temperatures were cooler and closer to the initial condition. The falling limb of storm event hydrograph corresponded with deviations in simulations at 22:00 hours of day 1 and 02:00 h of day 2 for cross sections 9 to 12, overestimating temperatures by approximately 0.3 ◦C and 0.5 ◦C, respectively.

**Figure 4.** Observed versus simulated river temperatures for the 12 cross sections (XS) and stations from 0 m to 1500 m in Sawmill Creek. The coefficient of determination, *R*<sup>2</sup> for each cross section is shown in the plot.

The i-Tree Cool River model simulated hourly water temperatures were not significantly different than the observed, for reach averaged data, based on the p-values calculated using a paired-samples t-test and the α = 0.05 (See Table S3 for more details). The 30-hour average simulated and observed river temperature, along the entire 1500 m Sawmill Creek reach, increased by 0.4 ◦C at a slope of approximately 0.02 ◦C per 100 m, but with longitudinal variation in that slope. For the 1500 m reach, the model had a root-mean-square error (RMSE) of 0.03 ◦C and a coefficient of determination (*R*2) of 0.98 with a p-value of 0.87 which was greater than the α of 0.05. The model simulated the relatively rapid increase in water temperature recorded by the sensors, between cross section 600 m and 900 m, corresponding to the reach with storm sewer inflow. This relatively rapid increase in temperature leveled at station 900 m, which is the first station downstream of the last storm sewer outfall. Relatively warm water in the Tannersville storm sewer entering Sawmill Creek between cross sections at station 600 m and at station 900 m was a major driver of the i-Tree Cool River Model forecasting a rise in river temperature during both wet and dry weather conditions (Figure S1a). During the wet weather, a total of 7 hours, the rate that simulated river temperature increased from station 600 m to station 900 m at a rate of 0.32 ◦C per 100 m (Figure S1b), much steeper than during dry weather. During dry weather, the simulated river temperature from station 600 m to station 900 m increased at a rate of 0.04 ◦C per 100 m (Figure S1c), approximately 12% of the wet weather slope. The i-Tree Cool River algorithms for shading, net groundwater discharge, hyporheic exchange, and upstream boundary condition temperature influenced the simulation of longitudinal river temperature and the model goodness of fit. In all of the scenarios, the calculated paired t-test p-values were smaller than the α = 0.05, rejecting the null hypothesis, *H*<sup>0</sup> of a significant difference between the means of the simulated and observed reach averaged river temperatures (See Table S4 for more details). When the shading algorithm was disabled, i.e., no shade was simulated, the i-Tree Cool River Model overestimated the river temperature for all cross sections by 0.34 ◦C, at a rate of 0.02 ◦C per 100 m, for the 30-hour period, 11 to 12 June 2007 (Figure S2), and the model RMSE increased to 0.36 ◦C and the *R*<sup>2</sup> decreased to 0.88.

Diurnal sinusoidal patterns of simulated and observed river warming and cooling were driven by the heat balance but disrupted by abrupt pulses of inflow due to warm runoff during the two storm events on 11 and 12 June 2007. The Sawmill Creek mean temperature, the average of measurements at the 12 cross sections, diurnally peaked at 15.8 ◦C by 15:00 h June 11, 2007 (Figure 5a), two hours after

the peak in shortwave radiation (Figure 5b). By 20:00 hours, shortwave radiation has declined to 0, net radiation became negative, and river temperature has decreased from the peak of 15.8 ◦C to 15.1 ◦C. A storm event at 21:00 h 11 June 2007, and then again at 01:00 h of 12 June 2007, generates inflow of warmer water from the upstream and storm sewer, creating a temporary increase in temperature, which disrupts the sinusoidal pattern in cooling toward the diurnal minimum temperature at 05:00 h of 12 June 2007. Dry weather extends through the remainder of the simulation, and at 06:00 h of 12 June 2007, the increasing shortwave and thus net radiation reestablish heat flux as the main driver of the increasing river temperature, which peaks at 16.4 ◦C at 14:00 h. There was no significant difference between the simulated and observed time averaged river temperature datasets based on the paired-samples t-test and α = 0.05 (See Table S5 for more details). The model simulations of the abrupt pulses in river temperature during the wet weather, extending from hour 20 of day 1 to hour 3 of day 2, had a Nash-Sutcliffe Efficiency (NSE) coefficient of 0.9 and a p-value of 0.80 (Figure 5a). The magnitude of the simulated river temperature changes due to the inflow of stormwater from the Tannersville storm sewer system was 0.3 ◦C during the first storm and 0.4 ◦C during the second storm. The model simulations had their poorest fit with observed river temperatures duringa6h period on 12 June 2007, between 03:00 and 09:00 h, centered at sunrise, when it overestimated the river temperature by an average of 0.13 ◦C.

**Figure 5.** (**a**) Reach average air temperature and observed and simulated river temperatures in Sawmill Creek, between 12:00 h of 11 June 2007 to 17:00 h of 12 June 2007; (**b**) Simulated reach averaged heat fluxes and precipitation into Sawmill Creek between 12:00 h of 11 June 2007 to 17:00 h of 12 June 2007.

The river temperature simulated by the i-Tree Cool River Model was a function of spatially and temporally varying contributions of groundwater, hyporheic exchange, and storm sewer inflow. Analysis of these components to thermal loading can assist in developing pollution mitigation or river restoration scenarios. In cross sections without storm sewer inflow and in the absence of rain events, the river temperature was predominantly determined by river flow from the upstream reach, and groundwater only contributed approximately 1%, while hyporheic exchange contributed approximately 10% (Figure 6a; 350 m, and 1100 m). In cross sections with storm sewers, even in the

absence of rain events, when the storm sewers discharged flow from illicit connections, they contributed approximately 25% of the flow influencing the cross section river temperature (Figure 6a; 800 m), which warmed the river water (see Table S2 reflecting warmer average temperature of the storm sewer in the dry weather). During wet weather, there was inflow from the storm sewer due to the flows from impervious areas, and this inflow contributed approximately 50% of the flow thereby influencing the river temperature (Figure 6b, 800 m) and provided the thermal load to the river system. Integrated along the 1500 m of Sawmill Creek reach, the contribution of groundwater summed to 15% of the total river volume, while the hyporheic exchange, which flows in and out within each sub-reach associated with a cross section, averaged approximately 10% of inflow in each sub-reach.

**Figure 6.** Simulated contribution to river temperature of the river flow (Str.), storm sewer (SS.), hyporheic exchange (Hyp.), and groundwater flow (GW.) in two timesteps including (**a**) before storm and (**b**) during the storm at three representative cross sections from the upper reaches (between cross sections at station 0 m and 600 m), middle reach (between cross sections at station 600 m and 900 m), and downstream reach (between cross sections at station 900 m and 1500 m).

In addition to the analysis of unsteady simulations in Sawmill Creek, the i-Tree Cool River Model performance was analyzed for the steady state condition in Meadowbrook Creek for 13–19 June 2012. The i-Tree Cool River Model simulated the time averaged river temperatures at 30 cross sections with an RMSE of 0.2 ◦C. We combined these 30 cross sections into reach averaged river temperature data to examine the diurnal pattern driven by the heat balance (Figure 7a). There was no significant difference between the simulated and observed reach averaged river temperatures based on a paired-samples t-test and α = 0.05 (See Table S6 for more details). The model simulations of the temperature for steady state in the Meadowbrook Creek study reach for 13–19 June 2012 had a NSE coefficient of 0.9 and a p-value of 0.72. The model captures how Meadowbrook Creek cools by 0.25 ◦C as it flows along the 475 m reach (Figure 7b), driven by the constant inflow of cooler groundwater.

**Figure 7.** (**a**) Reach average observed and simulated river temperatures in Meadowbrook Creek, for 13–19 June 2012; (**b**) Time averaged observed and simulated river temperatures for 475 m reach of the Meadowbrook Creek.

#### *3.2. Sensitivity Analysis*

A sensitivity analysis of the i-Tree Cool River Model examined the fluctuation in simulated temperature with changes in input data to identify the most sensitive parameters, which is useful when considering impacts of environmental change. The sensitivity analysis was performed for steady and unsteady simulations. We used global sensitivity analysis to identify the most important parameters and coordinated this analysis with that for the Meadowbrook Creek reach in summertime, by Glose et al. [31], noting both models are based on Heat Source [27]. Glose et al. [31] used an observed boundary conditions and our Equations (1), (2), (4)–(6), (10), (11), (S2), (S11) and (S12), and varied discharge by±10%, groundwater temperature ±15%, varied shading factor and view-to-sky factor by ±20%. They identified groundwater as the most sensitive parameter, with a ±0.2 ◦C change on average stream temperature. We replicated this sensitivity analysis to the 1500 m Sawmill Creek reach, confirming these sensitivities. We then extended the analysis in the 1500 m Sawmill Creek reach to consider varying parameters of storm sewer temperature, sediment temperature, and recorded boundary conditions temperature by ±15% (Figure S3) and varying parameters of substrate hydraulic conductivity (SHC), cloudiness factor (Cl), and groundwater discharge (GW) by ±20% (Figure S4). When storm sewer temperature was varied by ±15%, the reach-averaged temperature changed by 1.65% (0.27 ◦C). When sediment temperature was varied by ±15%, the reach-averaged temperature changed by 0.3%. When upstream boundary conditions temperature was varied by ±15%, the reach-averaged temperature changed by 9.5%. When substrate hydraulic conductivity was varied by ±20%, the reach-averaged temperature changed by 0.15%. When cloudiness factor was varied by ±20%, the reach-averaged temperature changes by 0.02%. When groundwater discharge factor was varied by ±20%, the reach-averaged temperature changes by 0.03%. Based on this analysis, the most sensitive model parameters, ranked in order of importance, are upstream boundary conditions, storm sewer temperature, sediment temperature,

substrate hydraulic conductivity, groundwater discharge, and cloudiness (additional sensitivity analysis is presented in supplementary materials Figures S3–S6).

#### **4. Discussion**

The i-Tree Cool River Model simulated the warming effects of the many potential sources of thermal pollution, including radiation fluxes and urban runoff, both dry weather illicit connections and wet weather stormwater. The model also shows how groundwater and hyporheic exchange inflows can provide a cooling effect, providing a comprehensive approach to assessing and perhaps mitigating thermal loading. To determine which factors are most effective in such management, this discussion provides some perspective on the effects of each warming and cooling effect.

The impact of urban runoff on the average temperature was rapid, within 1 hour of the onset of precipitation, and caused a temperature increase of 0.3 ◦C for the first storm, and 0.4 ◦C for the second storm of 12 June 2007. The rapid and large change in river temperature can be attributed to the short duration event, which Herb et al. [1] suggest when rainfalls only last 2 to 3 h will make the largest impact on raising stormwater and water temperature. The urban storm sewer area was approximately 21% of the watershed drainage area and had 35% impervious cover, which contributed to a relatively large volume of flashy, warm, stormwater response. Relative differences between air and water temperatures contribute to the warming, as noted by Herb et al. [1]; for the Sawmill Creek reach the average air temperature was 21.2 ◦C and average dew point temperature was 18.0 ◦C, both warmer than the average river temperature which was 15.2 ◦C. Even though the two rainfall events occurred at night during 12 June 2007, when solar radiation was not present, the prior day averaged 20% cloud coverage, allowing 80% of summer shortwave radiation to reach the small albedo impervious surface.

Simulation of the effects of the nighttime stormwater thermal load in riverine receiving waters on Sawmill Creek contributes additional data and tools to the investigation and management of the urban heat island. A common signature of the urban heat island is elevated nighttime air temperatures in urban areas relative to rural areas, due to physical differences affecting solar heating, such as albedo and thermal capacity, and anthropogenic heat sources [39]. Daytime insulation is a common driver of thermal loading of receiving waters [40], but for 12 June in Tannersville, the daytime solar heating of impervious area did not cause thermal loading of the river until the nighttime wet weather event. The nighttime precipitation landed on warm impervious surfaces, retaining much of their daytime elevated surface temperatures due to high capacitance, and this surface warmth was conducted into the stormwater entering the relatively cool, rural origin, receiving water. The effect of urban heat islands on rivers was studied by Somers et al. [41], who noted a 1.6 ◦C higher warm season temperature in urban rivers than forested rivers, and 8 ◦C greater spatial variation in urban rivers than in rural river temperatures along a 1 km transect. During a daytime storm event affecting all rivers, the temperatures in urban rivers rose as much as 4 ◦C, compared with a negligible rise in temperature in the forested rivers [41]. Nelson et al. [26] forecasted the thermal impact of individual storm events and found storm-induced river temperature surged by 3.5 ◦C for the warm season in urban watersheds near Washington, DC, USA; with drainage areas averaging 8 km2.

Proper simulation with the i-Tree Cool River Model of unsteady flows and their thermal pollution of receiving waters requires consideration of model goals and limitations. Typical model goals are either model inter-comparison for contrasting scenarios, such as varying impervious or tree cover, or model simulation for hindcasting or forecasting. In cases of model inter-comparison, the model has fewer limitations and the model physics will allow users to consider changes in river temperature for changes in study site conditions; model simulation has accuracy constrained by the accuracy of inputs as well as a model epistemic error [31,42]. This project attempted to improve accuracy of model simulations in Sawmill Creek by obtaining accurate data of the storm volumes and temperatures entering at the upstream and storm sewer locations along the boundary, using ibutton sensors, which are widely used for river temperature monitoring [31,43,44]. In cases where point or diffuse sources enter at multiple, unspecified locations along the river channel, such as with groundwater seeps, ibuttons may be inefficient and a better monitoring approach may involve using distributed temperature sensing system [8,35] for high-frequency time series, or from forward-looking infrared radar [45] temporally coarser data. An alternative to monitoring storm sewer inflow temperatures is estimating those temperatures using models of impervious runoff [25,29], and upstream boundary conditions can be estimated using air temperature records with the Mohseni et al. [30] non-linear regression.

Groundwater and hyporheic exchange were significant factors of temperature regulation during wet and dry weather. The section of Sawmill Creek simulated by the i-Tree Cool River Model had groundwater flow rates of approximately 0.0024 m3/s per 100 m, or 1% of flow, and riverbed longitudinal slopes that generated 10% contributions of hyporheic exchange (Figure 6). This combined subsurface flow, through the model inflow routines, contributed a cooling effect for the simulated summer period in 11–12 June 2007. Removing these inflows from the simulation caused the model to achieve RMSE of 0.18 ◦C and *R*<sup>2</sup> of 0.94, less important than the cooling through shading and the heat flux routines, when removed generated a RMSE of 0.36 ◦C and *R*<sup>2</sup> of 0.88. In winter, when river temperature is typically below subsurface water temperatures, this inflow would likely contribute a warming effect, as observed in other rivers by Risley et al. [46] and Kurylyk et al. [47]. From a survey of other studies, the relative contribution of groundwater and hyporheic exchange inflow with river water varies by site conditions and time. Poole et al. [48] working in mountain rivers, with bed slopes above 2%, also found the surface water received a larger volumetric inflow from hyporheic exchange than from groundwater, while Glose et al. [31] working in valleys with approximately 1% slopes did not identify significant hyporheic exchange and set groundwater as the only subsurface source of inflow.

Riparian shading from tree canopy, hillslope, and buildings provided the only land-based reduction in shortwave radiation and the view to the sky for the river, which influenced the longwave radiation. We used model inter-comparison simulations to contrast a scenario with and without shading and determined shading cooled river temperatures by an average of 0.34 ◦C during the 30 hours period. The landscape contribution to shading varied longitudinally along the reach, and at cross-section 9 was primarily from building shade, while upstream at cross sections 1 to 5 was primarily forest; hillslope topography provided minimal shading at this site. Shading is a concern in river thermal loading, and shallow and slow moving water is more vulnerable to such warming, and others have modeled this effect. Sun et al. [18] simulated 6 years along 6 separate reaches ranging in length from 85 to 1185 m of Mercer Creek in Washington State, and determined that tree and hillslope shading reduced the annual maximum temperatures by 4 ◦C. Roth et al. [49] simulated three cloud-free summer days in August 2007 along a 1260 m section of the Boiron de Morges River in southwest Switzerland, and determined riparian shading, by decreasing shortwave radiation, decreased daily average water temperatures by 0.7 ◦C. Guoyuan et al. [50] demonstrated predictions of shade from riparian vegetation (e.g., the Chen et al. [14] method used in i-Tree Cool River) were sensitive to the interaction of river azimuth and latitude, with E-W rivers in low latitudes benefiting least from riparian shade. Lee et al. [51] recommended for effective reduction in shortwave radiation, riparian areas utilize shading angles of 70◦ (1.22 radian) and view-to-sky factors smaller than or equal to 0.22. In the 1500 m Sawmill Creek reach, less than 10% of the view-to-sky factors were smaller or equal to 0.22 and shading angles averaged 50◦, and therefore additional thermal management opportunities are present.

The i-Tree Cool River Model was designed to assist river managers assess mechanistic causes of thermal pollution using free, open source, relatively simple algorithms in order to negotiate the balance between complexity and accuracy. While the model requires several input files, many of these can be obtained from publicly available data, site surveys, or estimation approaches; for model inter-comparison studies the accuracy of input data become less critical than in forecast simulations. The number of input files required by the model is comparable to other mechanistic models simulating river temperature, such as HFlux, HSPF, and QUAL2K, which require approximately 25 to 40 parameters, spatial data of river geometry and riparian features, and time series data describing the weather and discharge. Obtaining these inputs is a potential limitation of the i-Tree Cool River Model, and methods to obtain or estimate these input files are discussed above.

#### **5. Conclusions**

In this study, we developed the one-dimensional mechanistic i-Tree Cool River Model to simulate river temperature considering a combination of advection, dispersion, heat flux, and inflow processes. The i-Tree Cool River Model has the ability to analyze the impacts of external loads including multiple lateral storm sewer inflows, groundwater flow, and hyporheic exchange flow in steady and unsteady flows. The i-Tree Cool River Model estimates the shading effects of the riparian vegetation and other features as a function of heights and distances as well as solar geometry. The model performance was tested in steady and unsteady modes for the Meadowbrook reach in Syracuse, New York and Sawmill Creek in Tannersville, New York, respectively. The i-Tree Cool River Model performed satisfactorily in both simulations. The model can be used to conduct thermal pollution analysis of urban areas and investigate land cover and hydrology-based mitigation methods. The simulated river temperature of the i-Tree Cool River Model can be used for other environmental models, such as urban development models, atmospheric models, climate change models, and hydrology models.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/5/1060/s1, Table S1: List of the input files required for the simulation process of the i-Tree Cool River Model, Table S2: Observed water temperatures for three reaches of Sawmill Creek and for the Tannersville storm sewer, during 11 and 12 June 2007, as the average for all time steps during the dry or wet weather conditions, Table S3: Statistical analysis (paired t-test) of the reach averaged observed and simulated river temperature in Sawmill Creek for the (a) original condition including both wet and dry weather (b) wet weather, and (c) dry weather, Table S4: Statistical analysis (paired t-test) of the reach averaged observed and simulated river temperatures using the scenarios for the (a) no shading effect, (b) no groundwater and hyporheic exchange inflows, and (c) calculated boundary condition, Table S5: Statistical analysis (paired t-test) of the observed and simulated river temperatures in Sawmill Creek, between 12:00 h of 11 Jun 2007 to 17:00 h of 12 June 2007, Table S6: Statistical analysis (paired t-test) of the observed and simulated river temperatures in Meadowbrook Creek, for 13–19 June 2012, Figure S1: Time averaged observed and simulated river temperature in Sawmill Creek for the (a) original condition including both wet and dry weather (b) wet weather, and (c) dry weather, Figure S2: Time averaged observed and simulated river temperatures using the scenarios for the (a) no shading effect, (b) no groundwater and hyporheic exchange inflows, and (c) calculated boundary condition, Figure S3: Simulated time averaged river temperatures along the 1500 m Sawmill Creek reach for the original condition (Base) and for conditions with ±15% changes in (a) storm sewer temperature (TSS), (b) sediment temperature, and (c) boundary conditions temperature, Figure S4: Simulated time averaged river temperature along the 1500 m Sawmill Creek reach for the original condition (Base) and for conditions with ±20% changes in (a) substrate hydraulic conductivity (SHC), (b) cloudiness factor (Cl), and (c) groundwater discharge (GW), Figure S5: Fluctuations of shading factors and daily average shortwave radiation along the 1500 m Sawmill Creek reach. The shading factors denoted by a triangle are measured at each of the 12 monitoring stations, and the minimum and maximum shading factors were selected from the 5 m interval set of shading factors measured between each station, Figure S6: Temperature differences between the observed and simulated river temperature when using Mohsni et al. [30], Δ*Tcalc* versus recorded, Δ*Trecorded* boundary conditions, for (a) nighttime and (b) daytime.

**Author Contributions:** Conceptualization, R.A. and T.E.; methodology, T.E.; software, R.A.; validation, R.A. and T.E.; formal analysis, R.A. and T.E.; investigation, R.A. and T.E.; resources, R.A. and T.E.; data curation, R.A. and T.E.; writing—original draft preparation, R.A.; writing—review and editing, R.A. and T.E.; visualization, R.A.; supervision, T.E.; project administration, T.E.; funding acquisition, T.E.

**Funding:** This research was supported by the USDA Forest Service Northern Research Station i-Tree grant No. 15-JV-11242308-114.

**Acknowledgments:** We would like to thank Omid Mohseni, Laura Lautz, Ning Sun, and Jill Crispell for explaining their models.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **An Assessment of Self-Purification in Streams**

#### **Valentinas Šaulys, Oksana Survile \* and Rasa Stankeviˇ ˙ ciene˙**

Department of Environmental Protection and Water Engineering, Faculty of Environmental Engineering, Vilnius Gediminas technical university; Sauletekio al. 11, LT-10223 Vilnius, Lithuania; ˙ valentinas.saulys@vgtu.lt (V.S.); rasa.stankeviciene@vgtu.lt (R.S.)

**\*** Correspondence: oksana.survile@vgtu.lt; Tel.: +37-069-980-153

Received: 14 November 2019; Accepted: 21 December 2019; Published: 25 December 2019

**Abstract:** The territory of Lithuania is characterized by a prevailing moisture excess, therefore in order to timely remove excess water from arable lands, the drainage systems have long been installed. In order to drain excess water people used to dig trenches, to regulate (deepen or straighten) natural streams. The length of regulated streams has reached 46,000 km and they are deteriorated ecosystems. Investigations showed that the self-purification of streams from nitrates and phosphates is more effective in natural stretches than in stretches regulated for drainage purposes. Decrease in the average concentration of nitrates in natural and regulated stretches are 8.8 ± 5.0 and 3.0 ± 2.9 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup>−1, respectively. The average coefficient of nitrate self-purification, at a confidence level of 95% in natural stream stretches is 0.50 <sup>±</sup> 0.22, and in regulated is <sup>−</sup>0.15 <sup>±</sup> 0.21 km<sup>−</sup>1, and this difference is essential. The change in the average concentration of phosphates in natural and regulated stretches is almost the same, 0.2 <sup>±</sup> 0.1 and 0.2 <sup>±</sup> 0.2 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup>−1, respectively. The average coefficient of phosphate self-purification, at a confidence level of 95%, in natural stream stretches is 0.28 ± 0.12, in regulated <sup>−</sup>0.14 <sup>±</sup> 0.12 km−1, and this difference is not essential. In terms of the need for the renovation of drainage systems it is suggested that soft naturalization measures are first applied in the streams of Western (Samogitian) Highlands, Coastal Lowlands, and South-Eastern Highlands to improve their self-purification processes.

**Keywords:** water quality in streams; self-purification; nitrates; phosphates

#### **1. Introduction**

Small-sized rivers, their headwaters, and wetlands have a significant effect on the adjustment of river flow, retention of outwash and pollutants, and preservation of biological diversity. Unfortunately, intensive urbanization affects the channels of small rivers and adjacent ecosystems and strongly influences their ecological status [1]. Intensive use of European rivers for human purposes in the last hundred years has changed natural water flows, their physical and chemical properties, channel morphology, and species composition of local flora and fauna. The former natural, curving river channels were straightened and deepened, artificial slopes were formed having nothing to do with nature. Seeking to adjust flow rate, and due to recreational purposes, people built dams, and sometimes water energy was used in hydropower engineering [2–4].

At the end of the 20th century, with a constant growth of world population the cereals production almost doubled, the use of nitrogen (N) and phosphorus (P) fertilizers, related to the increasing food demand, has grown by 6.9 and 3.5 times. From environmental point of view, a large amount of N and P fertilizers was used improperly in the lands of agricultural designation and has essentially changed biogeochemical cycles of these two main nutrients [5]. The expanded network of small rivers and trenches was filled with a large amount of pollutants, including nutrients, pesticides, heavy metals, and soil particles washed out from adjacent fields. Due to the physical and chemical impact on river

channel and its water flow, part of local flora and fauna species became extinct, and part of them were replaced by invasive or introduced species [1].

Baczyk A. et al. [6] maintains that this negative anthropogenic effect has been widely described in many world-wide countries and even 96% of the articles analyzed reported a unilateral negative anthropogenic impact on water ecosystems, especially on invertebrates and fishes.

The territory of Lithuania is characterized by a prevailing moisture excess; the average annual amount of precipitation reaches 750 mm, evaporation comes to 512 mm [7,8], and therefore, the excess water is frequently accumulating in the soils of low-gradient and heavier gradation soils (clays, loam) of plains and depressions of hilly relief, due to low infiltration. To timely remove excess water from arable lands the drainage systems has long been installed, i.e., a complex of functionally related hydraulic structures located at the territory to be drained and intended for adjusting the soil moisture regime and creating favorable growing conditions for vegetation.

At the beginning of the 19th and 20th century, the largest attention when draining lands was paid to the regulation of rivers and streams and to digging of discharge trenches. Until 1958, lands were still drained by trenches, but later people started using drainage, since this was a world-recognized and more effective draining method. For the streams to be regulated, water intakes are usually curved, waterlogged, overgrown with water plants, having the already formed natural riparian lanes. Water depth and flow rate in these channels are usually low, water lies close to the surface, and floods take a long period.

Human activities, aimed to increase used agricultural and wooded areas, most of all changed small streams and the upper reaches of rivers: they were regulated to remove excess water collected by drainage, straightened to increase the flow rate and water capacity, or deepened to accommodate for at least 1.50 m depth of the river channel required for installing drainage system.

Thus, the regulated stream does not differ from a trench that was dug. The regulated streams and trenches served as a draining network.

During almost 100 years of land drainage in Lithuania, the streams were regulated, trenches were dug, and drainage systems were installed. Within this period, the area of drained land reached 3.02 million ha (47% of the total area of the country), of which 2.62 million ha were drained with the help of drainage, and 52,454 km of main drainage trenches were regulated and dug.

Land drainage has changed landscape structure, especially by a majority of newly dug trenches (Figure 1). The lowest density of trenches was found in plain regions where the most fertile clayey soils were drained. In hilly relief conditions the density of trenches remained higher.

**Figure 1.** Dynamics in regulated stream channels and main drainage trenches in Lithuania.

The length of drainage trenches had an effect on the density of the total Lithuanian hydrographic network. What concerns the total length of trenches in the territory of Lithuania it was obviously increasing. In 1930, the density of trenches was 0.1 km/km2, in 1945—0.34, in 1960—0.69, in 1975—0.76, and in 1999 even 0.96 km/km2. Such density of ditches was determined by the deepening, widening and straightening of natural stream channels, i.e., after they became main drainage trenches, and later—regulated streams. In such a way the length of regulated streams reached 46,000 km [8,9].

If the concept of regulated stream can be related to the adaptation of stream to serve draining function (to timely discharge certain amount of excess water necessary to drain the area), from the ecological point of view this is a deregulation of the natural stream and should be related to the destruction of the stream's ecosystem.

When carrying out regulation of streams the fact of deepening and straightening of natural stream channels used to be emphasized, though at the same time the riparian zones of streams were destroyed. There are not many studies focused on the retention processes of biogenic materials in protective riparian zones. Results of nitrogen retention in protective riparian zones are different [10–13], however water running through the ecosystem of riparian zone is able to retain up to 74% ± 4% of nitrogen. Data shows that with the increasing width of protective riparian zone the retention efficiency also increases [14], however this relationship is statistically stronger when the zone between the cultivated land and the river border >50 m.

Some authors argue that the efficiency of protective zones in retaining phosphorus is highly dependent on parameters such as vegetation cover type and density, topography, soil, climatic conditions, however, it is of short-term, seasonal in character, as long as plant vegetation continues [11]. The largest amount of total phosphorus is washed out in early spring before the start of vegetation. At that time, the efficiency of protective zone is poor, especially in the case where riversides are overgrown mostly with grass. The outwash is slightly better retained in protective zones overgrown with trees and bushes, though due to decomposition of organic material the amount of mobile phosphorus in these zones is increased [15]. Despite this, many agricultural consultants, planners and practitioners recommend the grass-overgrown protective riparian zones as a measure to settle phosphorus compounds in the areas of intensive agriculture.

There have been various recommendations for defining the main parameter of protective zones, i.e., width [16–18]. It was proven that a long-term efficiency can be reached where the zones are wider, 30 m-wide. In agricultural landscape, as suggested by foreign authors, the absolute minimum width should be 10 m. When designing protective riparian zones, it is necessary to take into account not only width requirements but also no less important parameters, such as vegetation species and distribution, slope gradients, soil types, topography, rainfalls [19,20].

In Lithuania, the minimum required depth of stream channel for the installation of drainage systems was about 1.5 m. In the upstream, where the basin area is not large, the stream channel has excessive depth and water collected from the area to be drained does not fill the whole cross-section area with water. Thus, there is a possibility to the improve ecological situation by growing bushes and small trees on the upper part of channel slope, and to manage the lower part of channel in a way that channel capacity corresponds to the requirements for draining function.

The authors suggest an alternative way of maintaining regulated streams (trenches): to allow the overgrowing of slopes with woody vegetation, the crowns of which would shadow the channel and strengthen the slopes [8,21,22]. Depending on the tendencies of species distribution in landscape and in cross-section of the trench, the restoration of dendrology in a desirable direction can be promoted artificially by planting special species or correcting their varietal composition in slopes, also by forming protective zones. This makes it possible when operating trenches to develop their natural functions and to also preserve drainage functions. The process of overgrowing of drainage network with trees can be assessed with ambiguity: reduces hydraulic capacity but makes a positive effect on landscape structure, reduces deflation probability, accumulation of outwash, and pollution of water bodies. Once the right of land ownership was given back to people, the land use intensity started to change. In less-favored areas for developing agriculture the drainage systems is used not intensively. Still, there are locations where drainage systems are outdated, their condition is bad. More than 45% of drainage systems are older than 40 years [23].

Over the last few decades, the number of projects of river naturalization has highly increased in most developed countries in Europe and all over the world [24–27]. Implementation of river naturalization projects is usually very expensive and the real benefit is not always evident. The analysis, carried out by German specialists [28], showed that the accomplished naturalization processes made no improvements to the population of benthic invertebrates. It was also reported that restoration of this population depends on whether there are any specimens of this species in adjacent waters. This research shows that the processes are not yet fully investigated and that slightly different results may be obtained under different natural conditions.

A very important indicator, showing an overall ecological status of the river, is self-purification [29]. It occurs due to water attenuation with surface and ground waters or due to certain hydrological, biological and chemical processes, such as sedimentation, coagulation, evaporation, sedimentation of colloidal and their further consolidation on the bottom of water body or, finally, due to pollutants assimilation with the living organisms. The level of self-purification in each water body depends on certain factors, such as temperature, water level, river flow rate, hydrological regime, tidal regime, amount of inorganic compounds in water, sediment characteristics, amount of pollutants, phytoplankton, benthic invertebrate fauna, fish fauna, and algae species, and their distribution [30,31].

The aim of the research is to determine the distribution of nitrate and phosphate concentrations in the water of natural and regulated for drainage streams and the influence of regulation on the self-purification efficiency of the streams.

#### **2. Materials and Methods**

#### *2.1. Location, Sampling Points*

The research was carried out in Lithuania, in the Baltic Sea region, in the Nemunas and Venta river basins. The research included regulated and unregulated stretches of the following streams: Terpine˙ (T), Žalesa (Z), Kuosine (K), M ˙ ekla (M), Durbinis (D), and Uogis (U). Distribution of the relevant ˙ stretches in the territory of Lithuania is given in Figure 2 and Table 1. The relevant stretches were coded. The first letter indicates the relevant stream according to the first letter of its name. The second letter indicates whether the stream is regulated or natural (the sampling place in the natural stream channel is marked by N, in the regulated—by R). The beginning of the stretch is marked by s; the end—by e.

**Figure 2.** Distribution of the relevant stretches in the territory of Lithuania.


**Table 1.** Characteristics of used streams.

#### *2.2. Time of Research, Nitrogen and Phosphorus Concentrations*

Water samples for water quality analysis were taken in accordance with the sampling standard [32] taking into consideration all the water sampling aspects. Water samples were taken once per month in the period from August 2013 to February 2019. The investigation in different stretches lasted 12–24 months. Concentrations of nitrates (NO− <sup>3</sup> ) and phosphates (PO3<sup>−</sup> <sup>4</sup> ) were determined. The amount of nitrates in water was measured by the photometer HANNA HI 83,205 using the cadmium reduction method and HI 93728-01 reagents. Concentrations of nitrates and phosphates in the samples taken were determined in the Laboratory of Hydraulics of Vilnius Gediminas Technical University.

According to the average annual value of each index a water body is attached to one of the five classes of ecological status. The obtained average annual values were compared to the values presented in the ecological status classification according to the physicochemical quality elements (Table 2) [33–35].

**Table 2.** Ecological status classification according to the physicochemical quality elements [33–35].


In the assessment many parameters (physicochemical, hydromorphological, biological) are used to determine the ecological status class. In the given article we have only rated according to two of them, NO− <sup>3</sup> and PO3<sup>−</sup> 4 .

For the assessment of river self-purification from biogenic substances, the following simplified river purification formula was used [36]:

$$\alpha = \ln \left( \mathbb{C}\_0 \mathbb{C}\_L^{-1} \right) \text{L}^{-1} \text{ \AA} \tag{1}$$

where: *C*0—concentration of chemical material at the beginning of the relevant river stretch mg L<sup>−</sup>1; *CL*—concentration of chemical material at the end of the relevant river stretch mg L−1; *L*—length of river stretch km; α-river purification coefficient km<sup>−</sup>1.

The statistical method of one-way analysis of variance (ANOVA), at confidence level of 95%, was used to determine whether the mean values of the self-purification coefficients for natural and regulated streams differed significantly. SPSS Statistics software was used for statistical processing.

#### *2.3. Distribution of Sediment and Used Agricultural Areas*

Distribution of sediment and used agricultural areas in the relevant river basins was studied according to the information provided in the Spatial Information Portal of Lithuania [37]. The portal provides data about typological units and gradation of the analyzed soil contour of the used agricultural area, also about land irrigation status and water logging, variety of soil cover, climatic conditions, soil stoniness, agrochemical properties, base index, and soil productivity index of the analyzed soil contour. However, the determination of the agents influencing the self-purification intensity of the stretches of the studied streams is very complex, due to the abundance and diversity of the factors and will be probably addressed in the future.

#### **3. Results**

According to the Spatial Information Portal of Lithuania, the main sediment in the relevant river basins is sandy loam and light loam which are prevailing in all basins. The river basins also contain light and medium loam, adhesive sand, light and average clay, and peat. Distribution of the used agricultural areas within the river basins showed that the largest part of the river basins is taken up by permanent grasslands, arable lands, forests, urbanized territory, and pastures (Table 3).


**Table 3.** Soil gradation, used agricultural areas and soil productivity index of the river basins.

The research and analysis of results showed that nitrate concentrations according to their amount and period are distributed rather differently. Many results showed that water quality in the relevant streams exceeded the limit value of good ecological status (nitrate concentrations 10.18 mg NO− <sup>3</sup> <sup>L</sup><sup>−</sup>1) and could be attributed to moderate and poor ecological status. The worst results were represented by Terpine, M˙ ekla, and Kuosin ˙ e streams. The best ecological status was found in Žalesa stream [ ˙ 38].

The worst water quality was observed in winter months since at this time of the year maximum concentrations of nitrogen oxides were determined which corresponded to the poor or bad ecological status: in Terpine stream, at point TRe (41.5 mg ˙ NO− <sup>3</sup> <sup>L</sup>−1, in December); in Mekla stream, at point ˙ MNs (81.4 mg NO− <sup>3</sup> <sup>L</sup><sup>−</sup>1, in March); MRe (76.4 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup><sup>−</sup>1, in February); in Kuosine stream, at point ˙ KRs (24.5 mg NO− <sup>3</sup> <sup>L</sup><sup>−</sup>1, in January); in Durbinis stream, at point DRe (26.1 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup><sup>−</sup>1, in December). During vegetation period nitrate concentrations were lower. This was especially characteristic of the Durbinis, Mekla (Figure ˙ 3), and Uogis streams (Figure 4). However, in the Kuosine (Figure ˙ 5) and Terpine streams the relatively large nitrate concentrations were measured even in vegetation period. ˙ This could be caused by a more abundant fertilization of arable lands in summer and autumn. The average annual consumption of fertilizers in Lithuanian agriculture in recent years (2016–2018) is 162,000 tons of nitrogen (N) and 52,000 t of phosphorus (P2O5). Winter crops were fertilized with 100–150 kg (N) ha−<sup>1</sup> and 50–75 kg (P2O5) ha<sup>−</sup>1, sugar beet 100–130 kg (N) ha−1, and 80–90 kg (P2O5)

ha−<sup>1</sup> [39]. As a result, a higher amount of nitrate pollutants is washed out into the rivers. In summer time, due to a shallow water the river water is less diluted.

**Figure 3.** Nitrate concentrations (NO− <sup>3</sup> ) in natural (N) and regulated (R) stretches of Mekla stream. ˙

**Figure 4.** Nitrate concentrations (NO− <sup>3</sup> ) in natural (N) and regulated (R) stretches of Uogis stream.

**Figure 5.** Nitrate concentrations (NO− <sup>3</sup> ) in natural (N) and regulated (R) stretches of Kuosine stream. ˙

The measurement results showed that phosphate concentrations according to their amount and period, same like nitrates, are distributed rather differently and water quality in streams frequently exceeded the limit value for good ecological status (phosphate concentrations 0.28 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup><sup>−</sup>1) and corresponded to the poor or bad ecological status (Figure 6). Very high phosphate concentrations were measured in the Durbinis stream where no seasonality was observed since phosphate concentrations were very high during all months. The maximum concentration measured was 2.5 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup><sup>−</sup>1. Based on this concentration the river water corresponded to the bad ecological status (>1.23 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup>−1). In Uogis stream (Figure 7), according to phosphate concentrations in water in a period of even four moths (March, May, June, and December) at the beginning of regulated stretch (point Urs) the ecological status of the stream was bad, however at the end of natural stretch (point Une) in a period of six months (March, April, June, July, November, and February) the stream had a prevailingly good ecological status, and in the remaining six months (May, August, September, October, December, and January) high ecological status. According to phosphate concentrations the Kuosine stream (Figure ˙ 8) represented high ecological status of water quality in January, February, June, and October as well as poor ecological status in December, March, and July in KRs stretches

**Figure 6.** Phosphate concentrations (PO3<sup>−</sup> <sup>4</sup> ) in natural (N) and regulated (R) stretches of Mekla stream. ˙

**Figure 7.** Phosphate concentrations (PO3<sup>−</sup> <sup>4</sup> ) in natural (N) and regulated (R) stretches of Uogis stream.

**Figure 8.** Phosphate concentrations (PO3<sup>−</sup> <sup>4</sup> ) in natural (N) and regulated (R) stretches of Kuosine stream. ˙

Based on the measured concentrations of nitrates and phosphates the average values at the end of stretches and self-purification coefficients were calculated in order to find out how the river is able to purify itself from pollutants. The average concentration of nitrates in natural stretches was twice lower than that in regulated stretches, i.e., 9.3 ± 2.7 and 18.5 ± 6.1 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup>−1, respectively. The average concentration of phosphates in natural stretches was also lower than that in regulated stretches, i.e., 0.4 <sup>±</sup> 0.2 and 0.7 <sup>±</sup> 0.2 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup><sup>−</sup>1, respectively.

The average change of concentrations, determined at the beginning and at the end of river stretches, and the calculated self-purification coefficients showed that streams purify better in natural stretches, further research into organic matter could, of course, further substantiate the purification efficiency. Decrease in nitrate concentrations in natural stretches was 8.8 ± 5.0, in regulated−3.0 ± 2.9 mg NO− <sup>3</sup> <sup>L</sup>−1, decrease in phosphate concentrations was 0.2 <sup>±</sup> 0.1 and 0.2 <sup>±</sup> 0.2 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup>−1, respectively (Table 4). The average nitrate self-purification coefficient of all streams was 0.50 ± 0.22 km−<sup>1</sup> in unregulated stretches, and 0.15 <sup>±</sup> 0.21 km−<sup>1</sup> in regulated stretches. The average phosphate self-purification coefficient in natural stretches was 0.28 <sup>±</sup> 0.12 km<sup>−</sup>1, in regulated <sup>−</sup>0.14 <sup>±</sup> 0.12 km<sup>−</sup>1.


**Table 4.** Average decrease in nitrate and phosphate concentrations in natural and regulated stream stretches and self-purification coefficients (R-regulated stretch, N-natural stretch).

<sup>1</sup> Average data of all relevant stream stretches.

What concerns separate streams, the best phosphate self-purification coefficient was determined in natural stretches of Uogis (0.74 <sup>±</sup> 0.34 km<sup>−</sup>1) and Žalesa (0.32 <sup>±</sup> 0.61 km<sup>−</sup>1). The nitrate self-purification coefficient in the natural stretch of Žalesa even reached 1.01 <sup>±</sup> 0.38 km<sup>−</sup>1, and in the natural stretch of Kuosine˙ <sup>−</sup>0.99 <sup>±</sup> 0.41 km<sup>−</sup>1. The negative values of self-purification coefficient were calculated in the regulated stretches of Terpine stream (Table ˙ 4).

A high value of self-purification coefficient in the natural stretch of Žalesa stream could be influenced by the adjacent permanent grasslands and forests. Decrease in pollutant concentrations can also be determined by the dilution of polluted water with surface and underground waters. It should be mentioned that the right tributary flows into the Žalesa in the stretch ZNs–ZNe. The inflowing stream water could dilute the nitrate-polluted water of the Žalesa and thus contribute to the decrease in nitrate concentrations. The natural stretch of the Kuosine stream is surrounded by forests. ˙ The self-purification was influenced not only by a natural river stretch but also by the adjacent used agricultural land, since it was likely that the access of nitrates into water was limited.

Nitrate self-purification coefficients are obviously different when comparing natural and regulated stretches of Terpine stream (Table ˙ 4; Figure 9).

**Figure 9.** Nitrate (*NO*− <sup>3</sup> ) self-purification coefficient in Terpine stream. ˙

The whole TNs-TNe stretch is a natural part of the Terpine stream abundant in meanders thus the ˙ river flow rate is lower and detention of nitrates is higher. Nitrate (NO− <sup>3</sup> ) self-purification coefficient α in a vegetation period (stretch TNs-TNe) is positive α = 0.10; 0.12 km−1. In winter season due to rotting vegetation and frozen land the increase in the amount of nitrates could be noticed in the stretch of natural channel <sup>α</sup> <sup>=</sup> <sup>−</sup>0.02; <sup>−</sup>0.07 km<sup>−</sup>1.

In the regulated part of stream, the amount of nitrates is frequently increasing in a flowing direction. In spite of that, in the samples of TRs the least amount of nitrates was determined in the whole period of research. This phenomenon could be explained by the fact that there is a dam before the stretch TRs–TRe and the existing pond operates as a settler, therefore retention of nitrates takes place all year round. The average self-purification coefficient of the stretch α is equal to 0.10 km<sup>−</sup>1. The slopes of the stretch TRs–TRe are very high and steep, a sharply descending relief increases river flow rate thus causing intensive bank erosion and slope landslips, whereas vegetation prevails on the upper part of the slope and in a protective riparian zone [40]. The average self-purification coefficient of the stretch <sup>α</sup> <sup>=</sup> <sup>−</sup>0.09 km<sup>−</sup>1.

The one-factorial dispersion analysis showed that the average self-purification coefficients, at a confidence level of 95%, vary essentially in the Žalesa, Terpine, and Kuosin ˙ e streams. The average ˙ self-purification coefficients, with a confidence level of 95%, do not essentially vary in the Durbinis, Mekla, and Uogis streams. Phosphate self-purification coe ˙ fficient in regulated and natural river stretches, at a confidence level of 95%, is not essential in the Kuosine, Durbinis, M ˙ ekla, Žalesa, and ˙ Terpine streams, though it varies essentially in Uogis stream. The results obtained show that nitrate ˙ self-purification is more effective than phosphate self-purification. The average nitrate concentration in natural stretches is lower than that in regulated stretches, 9.3 ± 2.7 and 18.5 ± 6.1 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup>−1, respectively. Decrease in the average nitrate concentration in natural and regulated stretches is 8.8 ± 5.0 and 3.0 ± 2.9 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup>−1, respectively. The average nitrate self-purification coefficient, at a confidence level of 95%, in natural stretches is 0.50 <sup>±</sup> 0.22 km−1, in regulated stretches <sup>−</sup>0.15 <sup>±</sup> 0.21 km<sup>−</sup>1, and this difference is essential.

The average phosphate concentration in natural stretches is also lower than that in regulated stretches, 0.4 <sup>±</sup> 0.2 and 0.7 <sup>±</sup> 0.2 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup>−1, respectively. The change in the average phosphate concentration in natural and regulated stretches is almost the same, 0.2 <sup>±</sup> 0.1 and 0.2 <sup>±</sup> 0.2 mg PO3<sup>−</sup> 4 L<sup>−</sup>1, respectively. The average phosphate self-purification coefficient, at a confidence level of 95%, in natural stretches is 0.28 <sup>±</sup> 0.12 km−1, in regulated stretches <sup>−</sup>0.14 <sup>±</sup> 0.12 km−1, and this difference is not essential.

The distribution of nitrate and phosphate concentrations in rivers is also influenced by the underground water, which we consider relatively small and similar, unfortunately it is not analyzed in detail in our research. However, in the future, in order to get a more detailed analysis of nitrate and phosphate changes, a more detailed assessment of groundwater (underground water fluxes, quantify the nutrient fluxes) should be carried out.

#### **4. Discussion**

The research results show that nitrate and phosphate self-purification in the regulated stretches of separate streams takes place slower or not at all since the self-purification coefficients were negative. Negative self-purification rates indicate that larger amounts of nitrate and phosphate pollutants are sometimes discharged into the river, although the inflow from the surrounding areas is similar but still variable, so regulated streams, sometimes natural also fail to treat itself. Seeking to prevent pollutants from getting into the rivers the protective riparian zones are defined along the riverbanks [41]. Industrial activities, carried out in protective riparian zones, has a direct, usually negative impact on water quality, therefore it is very important that they are observed. In order to improve water quality in rivers and streams of Lithuania it is suggested to naturalize them to the extent possible. After getting into the rivers and streams of Lithuania, nitrates and phosphates are better self-purified in natural than in regulated stretches. In order to restore a disturbed hydrological and environmental balance of regulated streams it is important to retain their drainage function without increasing maintenance costs. Therefore, seeking to improve self-purification of streams the so-called soft naturalization measures are suggested: to allow woody vegetation grow on river slopes, to form natural obstacles for water flow, and other elements in flood plains that are characteristic to wetlands [2,42]. The research evaluates the water quality of streams only in terms of nitrate and phosphate concentrations in water, but many

other physicochemical, hydro morphological, and biological parameters, for instance organic fraction of nutrients that also have a significant impact on a river s ecological status should be taken into consideration to fully assess the ecological status of the river and differences in the ecological status of natural and regulated streams.

Nitrate and phosphate concentrations are greatly influenced by man-or animal-made dams. The formed ponds act as precipitators, resulting that significantly less nitrate concentrations flow out from the pond [43]. Research in the Nevežis River Basin [ ˙ 22] shows that beaver ponds (comparing annual concentrations of inputs and outflows) retain relatively more phosphorus than nitrogen and more soluble mineral N and P compounds: on average 28% nitrate and ammoniacal nitrogen and 43% of orthophosphate phosphorus. But organic forms are much less retained. The retention of total N and total P is therefore 7% and 27%, respectively.

When planning renovation of drainage systems, the authors Šaulys and Barvidiene [ ˙ 44] suggest to take into consideration the Lithuanian hydrological regime, and natural and industrial conditions. Though the territory of Lithuania is relatively small, due to a variety of factors, which form and redistribute run-off, the feeding type and hydrological regime of surface water bodies are rather different. Lithuanian surface water bodies are divided into three hydrological regions based on relief, precipitation and soil [45,46]. Moreover, having taken into consideration the increase in soil productivity index due to drainage, abandoned agricultural land, possible breakdown of drainage systems, and dynamics in crop production, it was determined that the highest need for drainage systems renovation is in the Middle Lithuania region, the average in the Coastal Lowlands and South-Eastern Highlands, and the lowest in the Western (Samogitian) Highlands.

In terms of the need for renovation of drainage systems it is suggested that soft naturalization measures are first applied in the regulated stretches of the Durbinis and Uogis streams, and in other regulated streams of Western (Samogitian) Highlands. A more intensive application of naturalization measures is suggested for the Žalesa, Kuosine, and Terpin ˙ e streams, and other regulated streams of ˙ South-Eastern Highlands, since the mentioned streams are crossing territories where the prevailing soil is sand and sandy loam, moreover, the soil productivity index is relatively low.

The works focused on mechanical naturalization of regulated streams in Lithuania have already started [47], though mechanical naturalization of regulated streams financially is rather expensive, the benefit for water quality and ecological diversity has not yet been comprehensively investigated and scientifically justified [24,48,49]. Therefore, no intensive mechanical naturalization works are proposed so far for the regulated streams of Lithuania. In the future, the ecological status of regulated rivers should be studied and evaluated in as many parameters as possible, at the same time contributing to the broader public awareness of the rationalization of the application of naturalization processes.

#### **5. Conclusions**

It was determined that the self-purification of streams from nitrates and phosphates is more effective in natural stretches than in stretches regulated for drainage purposes.

Decrease in the average nitrate () concentration in natural and regulated stream stretches was 8.8 ± 5.0 and 3.0 ± 2.9 mg NO<sup>−</sup> <sup>3</sup> <sup>L</sup><sup>−</sup>1, respectively. The average nitrate self-purification coefficient, at a confidence level of 95%, in natural stretches was 0.50 <sup>±</sup> 0.22, in regulated stretches <sup>−</sup>0.15 <sup>±</sup> 0.21 km<sup>−</sup>1, and this difference is essential.

The change in the average phosphate concentration in natural and regulated stream stretches is almost the same, 0.2 <sup>±</sup> 0.1 and 0.2 <sup>±</sup> 0.2 mg PO3<sup>−</sup> <sup>4</sup> <sup>L</sup>−1, respectively. The average phosphate self-purification coefficient, at a confidence level of 95%, in natural stretches is 0.28 ± 0.12, in regulated stretches <sup>−</sup>0.14 <sup>±</sup> 0.12 km<sup>−</sup>1, and this difference is not essential.

During the warm vegetation period, concentrations of nitrates in water were (could be) in most cases lower than during the cold period of the year. Increased concentrations of nitrates could be caused by fertilization of arable lands. Increased nitrate concentrations could also be affected by the adjacent communities of garden-plots, this is a possible unevenness of inflow from surrounding areas.

In terms of the need for renovation of drainage systems it is suggested that soft naturalization measures are first applied in the streams of Western (Samogitian) Highlands, Coastal Lowlands, and South-Eastern Highlands to improve their self-purification processes and where the need for renovation of drainage systems is low.

**Author Contributions:** V.Š. conceived the study; O.S. and R.S. procured and analyzed the data, prepared the illustrations and wrote the manuscript under the supervision of V.Š. All authors contributed with ideas to the analysis, discussed the results and commented on the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
