Statistical approaches to explore the linkages between physicochemical parameters and BQEs, and set river nutrient threshold concentrations in Hungary

BQE

Biological quality element

BOD₅

5-day biochemical oxygen demand

COD_ps

Chemical oxygen demand with potassium permanganate method

DOC

Dissolved organic carbon

EQR

Ecological quality ratio

OLS

Ordinary least squares

PB

Phytobenthos

PP

Phytoplankton

RMA

Ranged major axis

TN

Total nitrogen

TP

Total phosphorous

WFD

Water Framework Directive

Water Framework Directive (WFD) (EC 2000) determines the European water policy with its main objective: all waterbody should achieve ‘good status’ and the implementation of sustainable water use. Currently, in Europe, more than half of the water bodies are in moderate or worse ecological status, primarily due to nutrient surplus from diffuse and point sources (Poikane et al. 2019a). Despite the taken efforts, the status of surface waters has hardly improved in the first cycle of river basin planning (EEA 2018), so there is an urgent need to establish reliable nutrient criteria consistent with biological status.

Water bodies are categorized into five status classes: high, good, moderate, poor, and bad. WFD (EC 2000) requires the Member States to establish threshold values for biological and physicochemical (such as nutrients) metrics for each class. The main purpose of this paper is to show how statistical methods can help in the harmonization of biological class boundaries which were previously defined in international intercalibration (Kelly et al. 2009) and nutrient class boundaries. There are no conventional methods for the determination of the phosphorus and nitrogen limit values, which have been set individually by each Member State. As a result, there may have been significant differences between countries' methods (GDWM 2016).

Many studies (Rask et al. 2011; Nõges et al. 2016; Huo et al. 2018; Phillips et al. 2018a) show that there is a significant relationship between stressors (such as nutrients and hydromorphology) and the condition of biological quality element (BQEs). On the other hand, studies proved that the pressure–response relationships are very complex (Lyche-Solheim et al. 2013), therefore deriving nutrient criteria might require more refined methods.

Understanding how water bodies interact with the environment is very important in terms of planning remedial measures. Water ecosystems can be simplified into univariate statistical models which are not an accurate description of the reality. However, multivariate, nonlinear models and deterministic models bring us closer to reality, this is why several studies use them to describe water ecosystems and hydrology as a complex system (e.g. Kansoh et al. (2020) compute water budget components for lakes from Egyptian meteorological data; Ćosić-Flajsig et al. (2020) present a holistic approach to define environmental flow in a Croatian river).

To resolve the contradiction of univariate and multivariate stressor–biological response relationships in surface waters, various statistical analyses were performed on Hungarian river monitoring data in order to test their applicability on threshold evaluation. This paper also demonstrates that random forest – which in our knowledge previously was not used for this purpose – can be applied for predicting biological class from physicochemical variables without setting their limits individually.

Exploring the linkages between physicochemical parameters and BQEs and setting river nutrient thresholds consistent with ecological status have not received enough attention in Hungary; thus, nutrient targets have not been set with statistical methods yet. Therefore, the aims of this study are to:

• 1.

  Derive nutrient thresholds consistent with good ecological status with statistical methods, where possible. Criteria assignment is made to phytoplankton (PP) and phytobenthos (PB), which are the most sensitive organisms for the nutrient load. This paper only considers the good–moderate threshold, which is most relevant in terms of the WFD, given that if a waterbody does not reach good status, actions need to be taken to improve it.

• 2.

  Explore relationships between several stressors and BQEs that are sensitive to these pressures according to literature.

• 3.

  Compare the results with other studies and discuss the advantages and disadvantages of the statistical methods used for threshold evaluation.

• 4.

  Show an alternative method that can be used for biological classification without setting limits for each parameter applied in the ecological classification.

The statistical analysis was carried out using a guide (Phillips et al. 2018b) and the related ‘Toolkit’. The purpose of the guide is to help the EU Member States to find nutrient concentrations that support good ecological status. As part of the ‘Toolkit’, we used the appended R code that contains several statistical tools. We also used RStudio 3.6.2 software to create boxplots and to build univariate and multiple regression, and random forest models (see the research methodology in Figure 1). Statistical analyses were performed for every Hungarian watercourse type. The types are differentiated according to altitude, size of the catchment area, geology, sediment roughness, and bottom slope (GDWM 2016). As the nutrient class boundaries that we were looking for are type-specific, we analysed each watercourse type separately and we also merged similar types where it was possible in order to increase the sample size.

With the aim of making a comparison, two independent databases were used for the analysis. One of them was the dataset of the River Basin Management Plan (GDWM 2016), which contains long-term average values of water body aggregated monitoring data from the period 2007–2012 for all 889 Hungarian watercourses. In parallel, the national monitoring open database (NEIS 2021) provided sampling-site-level data from the period 2007–2017 for 741 water bodies and 899 sample sites. All datasets contain measured concentrations (total phosphorous (TP), total nitrogen (TN), dissolved organic carbon (DOC), HCO₃, 5-day biochemical oxygen demand (BOD₅), chemical oxygen demand with potassium permanganate method (COD_ps), dissolved oxygen concentration, and water temperature), and PP and PB EQR (ecological quality ratio) values. Although the research focuses on PP and PB, to analysing the linkages between hydromorphologycal features and the state of BQEs, we used the aggregated monitoring data of all the five BQEs and data of the hydromorphologycal features of all Hungarian rivers from the dataset of the River Basin Management Plan (GDWM 2016). We chose the mentioned parameters for the analysis because of their assumed strong relationship with the BQEs and because of the suggestion of the guide (Phillips et al. 2018b).

With the sample-site-level database (NEIS 2021), analyses were conducted by using both year-round and vegetation period nutrient concentrations. Based on the type-specific information, natural, artificial, and heavily modified water categories were evaluated separately.

The ordinary least squares (OLS) and type II ranged major axis (RMA) regression were used. Both are univariate regression methods, the difference between them is how they estimate uncertainty in the dependent and independent variables (Legendre & Legendre 1998).

The good–moderate threshold concentration can be predicted from the regression equation using the good–moderate EQR threshold (Phillips et al. 2018b), which is 0.6 in Hungary. Also, regression methods can show how strong is the relationship between the variables involved in the analysis (Figure 2).

This is a categorical method that estimates the threshold value by minimizing the mismatch between biological status and nutrients and uses a series of nutrient concentrations. The method provides EQR into two categories (good or better and moderate or worse) (Phillips et al. 2018b).

Nutrient thresholds were generated by plotting the percentage of water bodies with mismatched classification in dependence of nutrients for potential nutrient criteria values. One pile of the lines represents where water bodies are at good or better status according to the EQR value, but moderate or worse status based on the nutrient concentration, and the other pile of lines represents the inverse of this (Figure 3). The cross-over point, where the two piles of lines meet, is where the misclassification is minimized and equal, hence this can be a potential nutrient threshold (Phillips et al. 2018b).

The advantage of this method is that it is easier to interpret and less sensitive to outliers than regression methods. The use of this method is advantageous in the case of poor nutrient–biology relationship (low r²) and if the range of data is not wide enough, i.e., we do not have data from all biological classes (Poikane et al. 2019b).

Boxplots were created to demonstrate the corresponding nutrient range in each biological class. On boxplots, boxes show the range of measured concentrations in each class from the first quartile to the third quartile with a marked median and vertical lines between the first quartile and the minimum and also from the third quartile to the maximum (Dawson 2011).

First, we illustrated all the five biological classes in respect of TN and TP concentrations for each watercourse type. To see the position of the biological classes with other stressors, boxplots were also made for the following variables with PP and PB biological classes: BOD₅, COD_ps, and additional hydromorphological parameters (such as hydrological regime, river continuity, and morphological conditions). We also made boxplots with PP and PB biological classes on which the position of the two BQEs to each other can be seen.

Besides nutrients, we analysed several stressors that have been found earlier as significant stressors with univariate and multivariate linear regression. Additional data of measured DOC, BOD₅, COD_ps, and information about hydromorphological features were involved in univariate analysis, while multiple regression was made with the combinations of TN and TP concentrations, PB, and PP EQR, and hydromorphological status of the waterbody.

To understand the behaviour and nutrient dependence of PP and PB, long-term EQR, TP, and TN time-series data were plotted. The graphical visualization of the data can show attributions that statistical methods conceal.

Random forest analysis was made with the randomForest package (Liaw & Wiener 2002). Random forest is a classification method that aggregates numerous decision trees by the bootstrapping method to give a prediction for the response categories. The method develops a predetermined number of trees. This was defined to 500 in each analysis. Each tree and each new split are made from random data selected by the bagging method (Breiman 2001). From the variable selection measures in tree development, we used the Gini index that indicates the frequency of the selection for a split for each variable and their overall discriminative value for the classification problem (Breiman et al. 1984). First, we ranked several variables (pH, alkalinity, N emissions from point sources, P emissions from point sources, zinc, DOC, BOD₅, HCO₃, COD_Cr, calcium, chloride, chromium, manganese, magnesium, inorganic nitrogen, organic nitrogen, sodium, NH₄-N, nickel, silicon, sulphide, sulphate, iron, nitrate, nitrite and phosphate concentrations, hardness, conductivity, TP and TN concentrations, temperature, total dissolved solids, slope, velocity, transparency, morphological status, hydrological status, status of longitudinal and transversal continuity, hydromorphological status, land use (percentages of urban, agricultural, and forest areas), category (natural, artificial, and heavily modified) water level, discharge, orthophosphate concentration, and oxygen saturation) that we considered having a big impact on the biota than we choose the five variables that have the biggest importance according to the Gini index for each watercourse type. Random forest prediction for the biological classes was made with the five chosen variables.

The application of regression and categorical methods demonstrated that the relationship between nutrients and BQEs is not strong. Nevertheless, according to the instructions of the EU guidance (Phillips et al. 2018b), we were able to perform statistical analyses and suggest nutrient class boundaries in several cases. With the dataset containing long-term averages, nutrient class boundaries were determined in 29 cases, and with sample-site-level annual averages, data class boundaries were suggested in 39 cases. Furthermore, by using data from the vegetation period only, class boundaries could be suggested in 63 cases. The threshold values varied in the case of highland rivers 1.8–6.2 mg TN/l, 180–400 μg TP/l; in the case of lowland rivers 1.4–5.0 mg TN/l, and 100–350 μg TP/l. Generally, we could suggest thresholds for those water types where a large amount of data were available (e.g. small and medium highland rivers with a mild bottom slope (359 waterbody, 1,380 annual average data, and 1,257 vegetation period data); small and medium lowland rivers with smooth sediment (376 waterbodies, 1,637 annual average data, and 1,510 vegetation period data)). However, in these cases, the r² values (p < 0.05) were lower than water types with fewer monitoring data. Despite the lower r² values in the case of small and medium rivers, the analyses of PP and nutrients yielded consistent and realistic thresholds (Table 1), those of TN are similar to what other studies described (Phillips et al. 2018a; Poikane et al. 2019b) and similar to current Hungarian nutrient thresholds (3–5 mg/L for TN and 150–300 μg/L for TP) (GDWM 2016). In contrast, the suggested thresholds for TP tended to be much higher than what the studies (Phillips et al. 2018a; Poikane et al. 2019b) described.

As the regression and mismatch methods did not lead to reasonable results in many cases, we created boxplots to see the nutrient (TP and TN) ranges in each biological class. The plots showed significant overlap between biological classes. Boxplots made with phosphorous tended to be ‘better’, the boxes overlapped less than in the case of nitrogen (Figure 4).

We also made boxplots with PP and PB biological classes from which we found out that apparently there is no relationship between the status of the two BQEs.

Boxplots and univariate linear regression were also made for assessing the relationship between PP and PB biological classes and BOD₅ and COD_ps. On the plots, boxes were better separated from each other, without significant overlapping. Statistical indicators were only acceptable (r² ≥ 0.36, p ≤ 0.05, Phillips et al. 2018b) in the case of small streams and big rivers.

DOC showed a stronger linear relationship with EQR than the other analysed variables (Figure 5). The relationship was the strongest in small streams and big rivers showed (0.2 ≤ r² ≤ 0.4, p ≤ 0.05).

Hydromorphological features were analysed by boxplots and regression analysis with all the five indicator organisms of the WFD (PP, PB, benthic invertebrates, macrophytes, and fish), since higher-order organisms reflect the hydromorphological status of the water flows better. Hydromophology was analysed separately as morphology, hydrology, longitudinal and transversal continuity, and overall, as hydromorphology. In general, the relationship between the involved variables was very weak (in most cases, r² ≤ 0.1, p ≥ 0.1) besides boxplots overlapped and did not figure in the expected shape.

The joint use of the two nutrient concentrations improved the relationship compared to the results of the univariate regressions (r² = 0.1–0.4). When hydromorphology was involved in the multiple regression model apart from nutrients, only r² ≈ 0.1 was obtained.

With the application of the random forest method, we were able to perform predictions for the biological class from five chosen variables selected according to their Gini index. For PP class prediction, BOD₅, water temperature, TP, TN, and HCO₃ concentration were used, while for PB prediction, we used BOD₅, water temperature, TP, TN, and dissolved oxygen concentration (Figure 6). With the mentioned variables, we were able to predict biological class with 35–81% error in the case of PP and with 18–47% error for PB (Table 2).

Analysing the long-term data series demonstrated that PP EQR values fluctuate significantly, even though there are not many outliers in nutrient concentrations outside of seasonality. It appears that there is no relationship between PP and nutrient concentrations in rivers. PB behaves differently, EQR values are very stable in time, unresponsive to changes in nutrient concentration. These phenomena appear in all watercourse types, as well as both in the case of nitrogen and phosphorus. Even when examining the same watercourses, it can be observed that PB and PP behave differently.

Presumably due to the small sample size and low reliability of the Hungarian BQE data (GDWM 2016), as well as high uncertainties in the methods, the relationship between BQEs and nutrients proved to be weak, weaker than other studies have found. Although based on the analysis with both sampling site and watercourse average monitoring data, we were able to establish nutrient thresholds in several cases.

By analysing each hydromorphological category (natural, artificial, and heavily modified categories) separately, the thresholds estimated by the statistical models became more unified and the value of r² grew higher, but still did not reach the target value of 0.36 defined by the guide (Phillips et al. 2018b). The value of r² was typically lower and the TP limits were higher than in the case of other countries (Phillips et al. 2018a, 2018b; Poikane et al. 2019b). Comparing the results derived from the year-round and vegetation period nutrient and EQR data, it can be said that although we obtained results in several cases using vegetation period averages, neither the r² values nor the uniformity of the estimated values suggests clearly that the use of vegetation period averages is more advantageous.

Boxplots showed no significant differences between nutrient concentrations in adjacent classes, boxes that symbolize different classes, often overlapped, and there were many outliers. Most of the cases boxes did not show the expected formation (‘bad’ biological class at the highest nutrient concentration, ‘high’ biological class at the lowest nutrient concentration).

The methods used for the analysis often resulted in a strong relationship only in the case of the water types with a small sample size. An also remarkable outcome is that nitrogen almost in every case led to a weaker relationship with BQEs than phosphorous. DOC, BOD₅, and COD_ps showed an even stronger relationship with BQEs than nutrients, while hydromorphology led to a very weak one.

The efforts made for a better understanding of the stressor–response relationship emphasized that physicochemical elements other than nutrients may have a stronger effect on the biological status of the rivers. Multiple regressions also show that these elements may weaken each other's effect. Among the examined variables, only DOC has a relatively strong but still statistically weak relationship with BQEs.

Random forest predicts the biological status with quite significant error in the case of some water types, but in other cases, predictions are appropriate. Trying to predict from different variables might decrease the errors.

The analyses highlighted that basic statistical methods did not provide reliable results; therefore, they are not suggested for determining nutrient criteria in Hungary for determining nutrient criteria in Hungary, though others (Phillips et al. 2018b; Poikane et al. 2019b) have found them capable for this purpose with European water monitoring data. One reason for the unsuccessful analysis lies in the small sample size and the low reliability of the Hungarian data set. There are only 889 river water bodies in Hungary, which is divided into 10 types and 15 subtypes. Some types include less than 10 watercourses, while others include more than 200. Statistical methods cannot be applied to types with only a few river bodies. Another reason might be that the suggested thresholds allow very high phosphorous concentrations in the Hungarian rivers, which value one order of magnitude above the limiting concentration of algal growth. In such conditions, there is no direct relation between BQEs and the nutrient concentrations; thus, the difference between the biological status in adjacent classes is not significant.

Also, given the complexity of water ecosystems, the shortage of basic statistical methods (such as OLS and RMA regression and the mismatch) is that they do not consider other factors that might influence the status of the biological elements; therefore, we suggest using them as decision support methods supplemented with more complex methods and expertise. Consequently, the guide (Phillips et al. 2018b) is also suggested only as an assistant throughout criteria establishment as the methods it presents most likely give a numerical result, but the value it gives might be an unrealistic nutrient threshold. A slight relationship that can be described with low r² (in the case of our study mostly r² < 0.3) is not enough convincing to provide reliable results for water quality status evaluation to be considered as a basis of acts and legislation.

The research also pointed out some mistaken implementation strategies of the WFD, namely that not each BQE in each water type reacts differently to pressures so considering every BQE throughout criteria development might be a mistaken strategy.

The paper demonstrated a new approach based on the random forest method that is a promising tool that can be used as a decision support method for biological classification, as it can be used as a nonlinear, multiparametric predictive model. It needs some refinement to decrease the error of classification, which might be done by using different variables, for example, residence time and shading. The benefit of this algorithm is that it can be used for biological classification without setting limits for each parameter applied in the ecological classification and without doing the really expensive and time-consuming biological sampling.

The research reported in this paper is part of project no. BME-NVA-02, implemented with the support provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021 funding scheme.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All relevant data are available from an online repository: NEIS (National Environmental Information System), 2021. Available from: http://web.okir.hu/en/ (16.01.2021) and from River Basin Management Plan (GDWM 2016).