Bayesian Machine Learning and Functional Data Analysis as a Two-Fold Approach for the Study of Acid Mine Drainage Events

Abstract

Acid mine drainage events have a negative influence on the water quality of fluvial systems affected by coal mining activities. This research focuses on the analysis of these events, revealing hidden correlations among potential factors that contribute to the occurrence of atypical measures and ultimately proposing the basis of an analytical tool capable of automatically capturing the overall behavior of the fluvial system. For this purpose, the hydrological and water quality data collected by an automated station located in a coal mining region in the NW of Spain (Fabero) were analyzed with advanced mathematical methods: statistical Bayesian machine learning (BML) and functional data analysis (FDA). The Bayesian analysis describes a structure fully dedicated to explaining the behavior of the fluvial system and the characterization of the pH, delving into its statistical association with the rest of the variables in the model. FDA allows the definition of several time-dependent correlations between the functional outliers of different variables, namely, the inverse relationship between pH, rainfall, and flow. The results demonstrate that an analytical tool structured around a Bayesian model and functional analysis automatically captures different patterns of the pH in the fluvial system and identifies the underlying anomalies.

1. Introduction

Acid mine drainage (AMD) is defined as the generation and transportation of highly acidic waters, which contain important quantities of heavy metals. The causes and detrimental effects that AMD has on water, soils, and biota are well known and have been reported in numerous scientific studies [1,2,3]. The low pH levels reached in the areas affected by this phenomenon impair the proper functioning of the vital processes in plants, humans, and aquatic biota [1,4,5,6,7]. Furthermore, the usually acidic pH of these effluents favors the leaching and transportation of heavy metals present in the geological formations impacted and also regulates the precipitation processes of these elements [3,8,9,10]. In the case of coal mines, the degree of acidity reached is highly conditioned by the percentage of sulfides present in the coal layers or in the till formations of the mineral layers [3,4,6,11].

AMD may initially be considered a localized problem restricted to areas near mine sites and mine waste dumps, but its influence can increase if acid mine water discharges into major watercourses [5,6,11,12]. Analysis and knowledge of the behavior of river systems affected by acid drainage and runoff are essential for the proper management of the environmental impact. In Spain, water quality monitoring and flood control programs, on a regional scale, are based on two automatic monitoring networks—SAIH (Automatic Hydrological Information Systems) [13] and SAICA [14] (Automatic Water Quality Information System) networks. The SAIH network collects data every day on several hydrological variables—precipitation, flow, and the surface water level—and the SAICA network records physicochemical parameters, which are indicators of water quality—water temperature, dissolved oxygen, electrical conductivity, pH, and turbidity. Both networks provide real-time information and valuable historical data series, whose correct processing and mathematical analysis can provide knowledge about the behavior of the river system and the detection of anomalous events. The use of water sensors for the improvement of water quality has been documented before [15,16].

The objective of this research is to evaluate the potential of two advanced mathematical techniques—Bayesian machine learning and functional data analysis—for the analysis and exploitation of knowledge contained in these valuable data series. These models will be implemented to achieve three specific goals: (i) the identification of behavioral patterns of the river systems affected by acid events, (ii) the inference of future behavioral responses, and (iii) the automatic detection of reliable anomalies. Bayesian machine learning models are statistical methods based on uncertainty quantification focused, in this case, on behavioral analysis and inference in complex multivariable systems [17,18,19] such as freshwater systems. On the other hand, functional data analysis [20] is a branch of statistics implemented, in this case, for treating continuous data series and detecting anomalies in different variables [21,22,23,24,25,26,27].

In particular, this research analyzes the case of the Cúa River during the last decade, which runs through the coal mining region of Fabero, located in the NW part of Spain. Mining activity in this area started in the 1910s and ceased definitively in 2018. However, water quality in the Cúa River is still affected by runoff and acid drainage, making it an environmental problem of special interest for the basin management organization. Therefore, the Cúa River is monitored by the two control networks available in this basin. In this work, the results of the Bayesian analysis describe a structure fully dedicated to explaining the behavior of the fluvial system and the characterization of the pH, delving into its statistical association with the rest of the variables in the model. Additionally, the resulting model solves inference problems based on Bayes’ theorem and the calculation of conditional probabilities, which enables predicting future anomalous behaviors of the fluvial system based on the available information from the SAIH and SAICA networks. The main findings of the functional analysis define several time-dependent correlations between the functional outliers of different variables, namely, the inverse relationship between pH, rainfall, and flow.

The novelty of this work is in the implementation and validation of Bayesian networks for: (i) the study of hidden behavioral patterns of the fluvial system dependent on the changes in the pH levels, (ii) the identification of the potential hydrological and physicochemical parameters which regulate these patterns, and (iii) the inference or prediction of future behavior of the fluvial system based on the quantitative variation of its control parameters. Likewise, the FDA algorithms are used for the automatic identification of anomalous events in SAIH and SAICA time series and for their capacity to explain the behavior of the system based on the analysis of anomalous events in their control parameters are validated.

This paper is structured as follows: the first section shows the theoretical framework of the mathematical techniques used and a detailed description of the case study and the data to be processed. Then, the main results of this work are presented in two distinct methodological modules: (i) those resulting from the use of Bayesian models for understanding the behavior of water and inference analysis and (ii) those resulting from the analysis of outliers using functional data analysis. Finally, a discussion of results is carried out to analyze the behavior of the fluvial system based on the results obtained. Lastly, the main conclusions drawn from both study approaches are presented.

2. Materials and Methods

2.1. Fluvial System Description and Data Acquisition

The data analyzed were recorded by a measuring station, which belongs to the SAIH and SAICA networks [13,14], and is located in the Miño-Sil River Basin. Both monitoring systems operate autonomously, collecting data on different variables relative to the control of the fluvial systems in Spain. In particular, the SAIH station measures water flow (m³/s), rainfall (mm), and temperature ($°\mathrm{C}$), which can affect the biological properties of the aquatic system, as well as the behavior of other substances dissolved in the water.

While the SAICA station records various physicochemical parameters of water quality, including: (i) conductivity (μS/cm), which measures the capacity of an aqueous solution to carry an electric current based on the presence of ions and their concentrations [28], (ii) dissolved oxygen (mg/L), which is essential for the breathing of aquatic organisms and enters the water through diffusion from the atmosphere and as a byproduct of photosynthesis [28], (iii) pH of the water to determine its acidity, using a logarithmic scale ranging from 0 to 14, with values below 7 indicating acidity, and values above 7 indicating basicity [28], and (iv) the amount of suspended particles in the water, known as turbidity (NTU), measured in Nephelometric Turbidity Units to determine the level of cloudiness or muddiness of the water [29].

More specifically, the station was positioned in the Fabero area, in the Alto Bierzo region within the Province of Leon, Spain. Its climate is temperate with a yearly mean temperature of 11.5 $°\mathrm{C}$, ranging from a mean maximum of 25 $°\mathrm{C}$ in July and a mean minimum in January of 0 $°\mathrm{C}$. Lastly, the pluviometry in this area is moderate, with a total of 850 mm/year. December and January are the rainiest months, while September and August register the lowest pluviometry values.

The economy of the Fabero region has been heavily influenced by the coal mining industry. The business started in 1917 because of the presence of anthracite coal and did not cease operations until November 2018. The Cúa River and its tributaries are the main hydrological resources of the area, and they are considered of special environmental interest by the water basin organization. The water quality in the area is adversely affected by acid mine drainage resulting from coal mining activities. The dumps containing pyrite undergo natural physical, chemical, and biological alterations, including spontaneous combustion and acid drainage production, releasing heavy metals into the environment and negatively impacting the landscape [30]. Acid drainage resulting from precipitation events reacts with rainwater before flowing to the nearest water body and traveling to the Cúa River, causing a gradual decline in pH and deterioration in water quality, particularly when the pH of most aquatic systems falls below 6, and especially below 5 [1].

Figure 1 shows a map of the Fabero area, including the location of the largest and last coal mine in the region, the position of the water control station, and the Cúa and Rioseco rivers.

The water quality station, which recorded the data for this research, is located downstream of the main drainage of the coal mine. The data analyzed spans from 1 January 2011 to 31 August 2022, giving a total of 11 years and 8 months recorded with daily frequency. This results in a total of 4342 data points.

2.2. Mathematical Background

The different mathematical methods adopted in this research are introduced in this section. Statistical process control is used as a first approach to identify anomalous events in the data and understand the general behavior of each variable studied; functional data analysis is implemented to achieve reliable detection of outliers and the discovery of relations between them. Lastly, Bayesian machine learning aims to find correlations between different variables and draw inference cases to help model and understand the trends and changes that control the Cúa River's fluvial system.

2.2.1. Statistical Process Control

Statistical process control is defined as the application of a set of statistical methods aiming to monitor and control the quality of a production process. Studying the variation of the data with these methods allows the detection of anomalies within the process. Control charts are the main tool in statistical process control, and they are used to represent and analyze a given process. They were initially created by Walter A. Shewhart [31] in the 1930s.

Significant variations may be easily identified by control charts. However, since they only analyze the most recent samples of groups, control charts fail to detect smaller changes or trends over a considerable time range. In order to address this problem, several new sets of rules were defined [32,33]. In this research, the enhanced version by Lloyd. S. Nelson [34] of the WECO rules [35] was implemented to increase the quality and precision of the results.

2.2.2. Bayesian Networks

Bayesian networks are directed acyclic graphical models (DAGs) that, for a given problem domain, represent a set of variables, called nodes, and their conditional dependencies, called arcs [36]. Their inherently graphical architecture leads to robust models ideal for easily exploring and interpreting complex problems. For the purposes of this study, the BayesiaLab v.10.2 AI software [37] was used to design a supervised Bayesian network based on the collected data in order to optimize the depiction and prediction of water acidification behavior.

To measure the variable impact on the target variable $\left(X\right)$, the relative mutual information (RMI) was calculated to quantify the percentage of knowledge obtained about $X$ by observing $Y$, thus being able to point out the predictors that give rise to a higher information gain [38]:

$$I_R\left(X,Y\right)=\frac{I(X,Y)}{H\left(X\right)}=\frac{H\left(X\right)-H\left(X\right|Y)}{-\sum _{xϵX}p\left(x_i\right){log}_2\left(p\left(x_i\right)\right)}$$

(1)

where $I(X,Y)$ represents the mutual information in absolute terms, and the entropy $\left(H\right)$ is an advanced quantification of uncertainty introduced in information theory as the amount of noise or disorder contained in a random variable. Entropy is maximal when all possible states of a distribution are equally likely. On the other hand, if there is no uncertainty, the network has zero entropy [39]. Finally, $x_i$ represents the states of node $X$, and $p\left(x_i\right)$ denotes the probability function.

To complete the characterization of the target node, the contributions of each variable to the predicted mean value of the target node were calculated. This metric was derived from counterfactuals based on the analysis of a specific neutral state representing the neutral condition of the system.

Finally, the Contingency Table Fit (CTF) metric was used to measure the quality of the representation of the Joint Probability Distribution (JPD) via the target variable. A common alert threshold below which the target value should be relearned is 70%. From a mathematical viewpoint, CTF may be expressed as [40]:

$${CTF}_N=100\frac{H_D\left(X\right)-H_N\left(X\right)}{H_D\left(X\right)-H_C\left(X\right)}$$

(2)

where $H_N\left(X\right)$ represents the entropy value of the constructed Bayesian network, and $H_C\left(X\right)$ and $H_D\left(X\right)$ are the entropy values of the fully connected and disconnected network, respectively.

2.2.3. Functional Data Analysis

The analysis of functional data is applied to observations of a continuous random process that is observed at discrete points [41]. In the case of an initial sample composed of different data points, the data must be transformed to generate a functional sample. This process consists of defining the best-fitting functions or curves to the initial data points that discretely represent the experimental process. This transformation is known as smoothing.

Considering a set of samples $x\left(t_j\right)$ in a set of $n_p$ points $t_j\in \mathbb{R}$, where $t_j$ represents every time instant, all samples can be considered discrete observations of the function $x\left(t\right)\in \chi \subset \mathrm{F}$, with $F$ being a functional space. The function $x\left(t\right)$ is estimated by taking into account that $F=span\left\{{\varphi }_1,\dots ,{\varphi }_{n_b}\right\}$ is a functional space consisting in a set of basis functions $\left\{{\varphi }_k\right\}$, with $k=1,2,\dots ,n$, where $n$ is the required number of basis functions to define the functional space $F$. Several types of basis exist in statistics, but given the fact that the data used in this research has periodic nature, the Fourier basis becomes the most suitable [42,43]. Once the different parameters have been defined, the smoothing problem may be expressed as [26]:

$$\underset{x\in F}{\min }{\sum }_{j=1}^{n_p}{\left\{z_j-x\left(t_j\right)\right\}}^2+\lambda \Gamma .$$

(3)

In Equation (3), $z_j=x\left(t_j\right)+{\epsilon }_j$ is the result of evaluating $x$ at the point $t_j$, ${\epsilon }_j$ is the zero-mean random noise, and $\Gamma$ is an operator which penalizes the complexity of the solution. The purpose of this penalty is to ensure a good fit to the data in the sense that ${\left\{z_j-x\left(t_j\right)\right\}}^2$ is small, but also some aspect of the data captured by $\Gamma$ is kept under control. Lastly, the $\lambda$ parameter adjusts the intensity of the regularization. After the discrete data is converted to its functional counterpart, outliers, defined as portions of the data that present different parts and/or values, can be detected based on the concept of functional depth. This statistical measure was first introduced for multivariate analysis to quantify the centrality of a given point compared to a set of observations. Consequently, the points represented in a Euclidean space that is close to the center of the distribution will have a greater depth value compared to those located in the periphery [44]. This definition has been extended to the functional domain, where the depth measure analyzes the centrality of a curve $x_i$ compared to a set of curves $x_1,\dots ,x_n$. Therefore, the concept of depth makes it possible to work with observations, defined in a given time interval, in the form of curves, instead of having to summarize the information contained in these curves into a single value, such as the mean [25]. In this case, the Modified Band Depth (MBD) [45] has been selected as it has demonstrated a better performance in the analysis of environmental data with this approach [46]. Functional directional outlyingness is considered to increase the accuracy in the detection of outliers. This method separates functional outlyingness into two main components: shape outlyingness and magnitude outlyingness, allowing the study of the centrality of the curves and their variability. A magnitude outlier is defined as an observation that is shifted from the mass of the data, while a shape outlier is an observation that differs in shape from the rest of the functions in the set (even if it lies completely inside the bulk of the data) [47]. In this work, the detection of anomalies is performed on the magnitude and shape value pairs according to the outlier detection algorithm introduced by Dai and Genton [48].

3. Results

3.1. Variability of the Data and Outlier Detection with SPC

The complete dataset of all variables was initially analyzed with the $\stackrel{-}{x}$ control chart, which displays the mean of all values in weekly subgroups. Each dot in the plot corresponds to a specific week of the database. Among all variables available, the pH is regarded as the main indicator of the acid mine drainage problem in the Cúa River, while the flow acts as a central measure of the environmental dynamics in the fluvial system. The result of the 8 Nelson rules on the pH $\stackrel{-}{x}$ control chart is shown in Figure 2. This information helps visualize the data and its variability, besides serving as a first approach to the anomaly detection issue. Nevertheless, the overfitting of the Nelson rules decreases the validity of the results for this purpose.

The pH data, as displayed in the $\stackrel{-}{x}$ control chart of Figure 3, does not present a clear behavioral cycle, and most values are within 1 standard deviation from the mean. Furthermore, Nelson rule number 1 detects a total of 8 weeks with a pH mean value below 3 standard deviations from the mean and 2 weeks with a mean value over 7.86, which is 3 standard deviations above the mean. Nevertheless, the rest of the rules detect an unrealistic number of alarms, as it represents more than 50% of the pH data.

Another example of this analysis is presented in Figure 3. In this case, the variable water flow is analyzed with the $\stackrel{-}{x}$ control chart. Contrary to the pH data, the water flow displays a clear pattern, marked by a fast increase from December until reaching peak values in the winter months of the Northern Hemisphere and a steep decline in spring to the low and steady water flow levels of the remainder of the year. Additionally, Nelson rule number 1 performs a correct detection of the most outstanding anomalies in accordance with the functional results, but rules number 2 and 3 define as outliers a disproportionate number of weeks, mostly, all those values below the mean which correspond to normal flow in the months of summer, fall, and the end of spring.

Consequently, the results of the statistical process control do not present a degree of reliability high enough for accurate anomaly detection in these data. However, they shed light on its variability, the seasonal patterns of the Cúa River, and the different standard deviation intervals defined, which will be used to classify individual values in different categories in the Bayesian analysis.

3.2. Variables Influence Analysis in pH Distribution

In this section, a supervised Bayesian network representing the general problem under study was constructed by analyzing the statistical association stemming from the pH status and each variable in the model. For that, considering the entire database, the standard deviation intervals defined from the $\stackrel{-}{x}$ control charts of each variable in the database were used as discretization threshold values for the Bayesian analysis (

Table 1. Threshold values selected from the $\stackrel{-}{x}$ control chart analysis.

Conductivity (µS/cm)

$\le 106.89 (x \le \stackrel{-}{x}$ − 1σ) | $\le 204.04 (\stackrel{-}{x}$ − 1σ,$\stackrel{-}{x}$] | $\le 301.2 (\stackrel{-}{x}$,$\stackrel{-}{x}$ + 1σ] | $> 301.2 (x>$$\stackrel{-}{x}$+ 1σ)

Dissolved oxygen (mg/L)

$\le 8.81 (x \le \stackrel{-}{x}$ − 2σ] | $\le 9.79 (\stackrel{-}{x}$− 2σ,$\stackrel{-}{x}$ − 1σ] | $\le 11.7 (\stackrel{-}{x}$ − 1σ,$\stackrel{-}{x}$+ 1σ] | $> 11.7 (x>$$\stackrel{-}{x}$+ 1σ)

pH (u. pH)

≤5.5 [49]$| \le 6.5 (5.5$,$\stackrel{-}{x}$ − 2σ] | $\le 7.1 (\stackrel{-}{x}$ − 2σ,$\stackrel{-}{x}$] | $> 7.1 (x>$$\stackrel{-}{x}$)

Water Temperature (°C)

$\le 6.73 (x \le \stackrel{-}{x}$ − 1σ] | $\le 10.07 (\stackrel{-}{x}$− 1σ,$\stackrel{-}{x}$] | $\le 13.4 (\stackrel{-}{x}$,$\stackrel{-}{x}$ + 1σ] | $> 13.4 (x>\stackrel{-}{x}$ + 1σ)

Turbidity (NTU)

$\le 4.28 (x \le \stackrel{-}{x})$$| \le 11.35 (\stackrel{-}{x}$$, \stackrel{-}{x}$+ 1σ) | $\le 18.42 (\stackrel{-}{x}$+ 1σ$, \stackrel{-}{x}$+ 2σ$] | > 18.42 (x > \stackrel{-}{x}$ + 2σ)

Rainfall (mm)

$\le 2.29 (x \le \stackrel{-}{x}$) |$\le 5.68 (\stackrel{-}{x}$,$\stackrel{-}{x}$+ 1σ] | $\le 9.07 (\stackrel{-}{x}$+ 1σ,$\stackrel{-}{x}$+ 2σ] | $> 9.07 (x>$$\stackrel{-}{x}$+ 2σ)

Temperature (°C)

$\le 5 (x \le \stackrel{-}{x}$ − 1σ] | ≤ 11.15 °$\mathrm{C} (\stackrel{-}{x}$ − 1σ,$\stackrel{-}{x}$] | ≤ 17.3 °$\mathrm{C} (\stackrel{-}{x}$,$\stackrel{-}{x}$ + 1σ] | > 17.3 °$\mathrm{C} (x>$$\stackrel{-}{x}$+ 1σ)

Water flow (m3/s)

$\le 5.54 (x \le \stackrel{-}{x}$) | $\le 11.38 (\stackrel{-}{x}$,$\stackrel{-}{x}$+ 1σ] | $\le 17.22 (\stackrel{-}{x}$+ 1σ,$\stackrel{-}{x}$+ 2σ] | $> 17.22 (x>$$\stackrel{-}{x}$+ 2σ)

).

From the resulting Bayesian network, a relationship analysis was carried out, between the target node (pH) and its child nodes, with the goal of identifying the most representative parent–child connections. These relationships are shown in Figure 4, where water temperature (35.65%) and flow (31.29%) stand out for their high impact, followed by turbidity (10.68%), air temperature (8.93%), and conductivity (8.09%) variables. To clarify the interpretation, the RMI was also computed exclusively between the arcs of the variables and the target node. The blue number shows the RMI concerning the secondary node, while the red number refers to the main node. Thus, knowing the pH state reduces the uncertainty of the conductivity characteristics by 18.9196%. In the opposite direction, river flow is the factor that provides the most significant reduction in the uncertainty of the pH state (16.5696%). Similar results are found for turbidity and conductivity, with 12.7214% and 12.1323%, respectively. From the point of view of the predictive importance of the variables, dissolved oxygen is the factor that provides the lowest amount of information in the case analyzed.

To go deeper into the problem of acid mine drainage in the Cúa River,

Table 2. Local impact analysis with pH target states. Results of the Relative Binary Mutual Information (RBMI) calculated.

	≤5.5 1	(5.5, 6.5]	(6.5, 7.1]	(x > 7.1)
Conductivity	5.271%	26.418%	10.594%	20.361%
Flow	17.072%	27.955%	11.091%	19.603%
Month	15.260%	21.067%	7.699%	11.168%
Turbidity	16.980%	23.783%	3.183%	8.707%
Water Tª	7.268%	9.921%	3.614%	5.945%
Temperature	7.907%	10.968%	3.494%	5.793%
Rainfall	10.635%	9.411%	0.672%	2.411%
Dissolved oxygen	4.291%	4.961%	0.919%	1.629%

Notes: ¹ A pH level of 5.5 is the limit defined by Spanish legislation [49] for the change of class condition from Good/Moderate to Moderate/Deficient for the river type: R-T31 small Cantabrian–Atlantic siliceous axes.

presents the calculation of the amount of information brought by each variable to the knowledge of the different ranges of pH values. Previous high-impact studies have shown that freshwater quality begins to gradually deteriorate as we move away from the normal value ($\stackrel{-}{x}$). This deterioration in water acceptability is most noticeable when pH values fall below 6 and, in particular, below 5 [1]. Likewise, for maximum productivity, Simate and Ndluvo emphasize that the pH value range be maintained between 6.5 and 8.5 [1].

That said, this local analysis between each of the 4 states of the target variable and its predictors allowed us to identify that flow in pH ranges below 5.5 [49] is the factor with the highest correlation. Given the importance of freshwater quality when the pH value drops below 6.5, the percentage values of mutual information obtained for conductivity, flow, turbidity, and month variables in the pH range between 5.5 and 6.5 should be highlighted. However, if the pH level is higher than 7.1, the most important predictor variables would be conductivity (20.361%) and flow rate (19.603%), respectively.

In terms of model validation, CTF is a very useful metric to measure the quality of the JPD representation. A CTF higher than 70% implies an adequate quality of the induced factors, and specifically the value obtained for the supervised Bayesian network created is 74.581%.

3.3. Functional Analysis Approach

This section presents the results of the functional analysis of each variable for outlier detection purposes. Most notably, several relationships are drawn between two or three given variables of the database related to the date of their respective outliers. In accordance with the results of the mutual information analysis between each variable and the pH, the functional results begin with the analysis of the water flow variable, whose graphical results are displayed in Figure 5. The left side of this figure contains a Cartesian representation of the magnitude–shape pair of values of each function. As can be seen, those functions with an anomalous magnitude and/or shape outlyingness value (red points) lay outside the main cluster of points, which contains the pair of values magnitude-shape of the nonoutlying functions (blue points). The right side shows the functional plot of the weekly rainfall values. Those functions that are considered anomalous are colored in red, while those considered standard are represented in blue. Water flow has a mean value of 5.35 m³/s, with a maximum of 67.37 m³/s on 20 December 2019, and a minimum of 0.46 on 9 October 2017, with most outliers taking place between December and March of every year.

Furthermore, Figure 6 shows the results of the functional analysis on the pH data in the same fashion as the previous image. This variable has an average value of 7.07 with a minimum of 4.3 on 5 February 2017, and a maximum of 8.0 on several days in June and July 2020, with most acid anomalies happening also between December and March.

This analysis was performed on the data of each variable, and the results are summarized in

Table 3. Results of the functional analysis on all variables. The first column contains the name of each variable. Columns 2 to 4 include information on the minimum, average, and maximum values for their respective variables. Lastly, column 5 contains the number of weeks identified as outliers in each variable out of the 620 analyzed.

Variable	Min.	Avg.	Max.	Outliers
Rainfall	0 Several days	2.61	93.2 10 December 2017	84
pH	4.3 5 February 2017	7.07	8.0 Several days	84
Flow	0.46 9 October 2017	5.35	67.37 20 December 2019	83
Conductivity	21.8 13 December 2020	204.10	642.2 14 October 2011	81
Temperature	−4.0 8 January 2021	11.26	26.0 17 June 2017	83
Dissolved oxygen	7.58 11 November 2016	10.76	13.91 8 November 2016	82
Turbidity	0 Several days	4.27	199.90 6 May 2012	81

. This table contains basic statistics of all variables analyzed and the number of weeks identified as functional outliers.

After the functional outliers are detected, it can be inferred that there could be some degree of correlation between the variables, as the number of outliers is around the same value in all cases. The existence of anomalous events on the same dates and in several variables defines a multivariate anomaly, and its detection validates the functional analysis for the identification of such events in the data provided by the SAIH and SAICA networks. These results are complemented with a time dependency analysis of the detected outliers and the correlations between the different variables. The rest of this section presents the relationships found between two or more variables based on the anomaly detection results of the functional analysis. The main outlier correlations found are presented graphically in Figure 7.

Flow and pH contain 26 coincidental outlying weeks, giving a match percentage of 32% between those 2 variables. They present an inverse relationship since when the flow increases, the pH levels drop. Weeks from numbers 17 to 20 seem to break away from the defined trend. However, these weeks, which belong to November and December 2019, present mean flow levels between 14.31 m³/s and 24.71 m³/s. Consequently, their values are considerably above the 3rd quartile of the flow data, which is 6.53 m³/s, evidencing their outlying behavior. This information points towards the fact that the pH level can decrease, but when it reaches a certain point, the water flow loses its influence. In other words, there are certain scenarios when an increase in the flow level will not necessarily imply a further drop in the pH.

Rainfall and pH have a total of 23 outlying weeks in common, which translates into a 28% match. It can be observed that their relationship is inverse, meaning that when it rains, the pH levels drop. Nevertheless, there is an exception to this rule: outlying week number 22 (10 December 2020 $\to$ 16 December 2020), when the rainfall values decrease in synchronization with the pH levels. This behavior is explained by analyzing the precipitation levels in the previous week, which are around a high value of 7.1 mm each day. Additionally, in week number 22, there is a peak of 17 mm and 26 mm on Monday and Tuesday, but Wednesday, Thursday, and Friday are close to 0 mm, lowering the mean value of that week. Thus, it can be inferred that when a precipitation event takes place, the pH levels do not drop immediately, and when the rain ceases, the pH level does not go back up to its standard numbers after at least 3 to 7 days.

Turbidity and pH present an inverse relationship, which matches with the behavior of the fluvial system, meaning that a higher level of dissolved particles implies acid waters with a lower pH value. The correlation between these 2 variables translates into a matching percentage of 27% with 21 outlying weeks shared between turbidity and pH. Conductivity and pH have a clean direct relationship: when the conductivity decreases, so does the pH and vice versa. These 2 variables have a matching percentage of 40% with a total of 32 coincidental anomalous weeks. Temperature and pH have 6 outlying weeks in common, resulting in a matching percentage of 9%. Nevertheless, their proper direct relationship can be observed in Figure 7e.

Temperature and dissolved oxygen include a higher number of common outliers for a total of 33, resulting in a match percentage of 41% and a well-defined inverse relationship since when the temperature rises, especially in the summer months, the dissolved oxygen declines. Lastly, three variables were found to have a direct relationship: rainfall, flow, and turbidity. They present a matching percentage of 40% due to their 32 coincidental outliers. This relationship can be seen in Figure 7g.

Paying close attention to the pH variable, which is paramount in this analysis, the total number of outliers explained by the relationship between the pH and any of the other 6 variables in the database is 55, equivalent to 65% of the pH anomalies. This implies that there is a 35%, or 29 outliers, of the 84 total pH anomalies unexplained by their correlations with the other variables. Out of these 29 outliers, 24 are weeks with a mean value above the average. Therefore, these anomalies, of which 61% take place in the driest months of the year between May and October, could be due to a natural process whose information is not contained in the current variables available or due to uncontrolled dumps of polluting elements in the Cúa River. Regarding the 5 outliers remaining unexplained, these present a mean value lower than the mean pH of 7.07, but after careful analysis of the rainfall in the weeks prior and after those anomalies, it can be concluded that there is some lag between the rainfall anomaly and the associated acid pH outlier. In these circumstances, every time an important precipitation begins, the acid pH anomaly does not take place until approximately 1.5 days on average have passed. Similarly, when the rainfall stops, it still takes typically 4.5 days for the pH levels to rise.

4. Discussion

The hydrological and environmental data of the Cúa River was analyzed with three different statistical methods. The results of the control charts, in particular, the $\stackrel{-}{x}$ chart provides valuable information on the variability of data, the seasonality of the different natural processes taking place in the river, and the different standard deviation intervals defined by the control limits. On the other hand, the detection of anomalies through this method does not offer consistent results. Control charts consider the temporal correlation between data, but even when the use of subgroups is helpful for graphical representation, it also accounts for an important loss of information. Additionally, the seasonality of the environmental data damages the performance of these methods as the anomaly detection rules become erratic. It is this set of conditions that leads to the detection of an unrealistic number of outliers, marked by the definition of frequent values as anomalies.

Concerning Bayesian analysis, the use of artificial intelligence algorithms in combination with other mathematical approaches is expected to have a positive impact on public water protection policies in the coming years. These models have a high capacity and flexibility for the analysis of complex problems and can handle the heterogeneity of hydrological circumstances. As an example of the potential of this combined approach, Figure 8 shows the Bayesian results of an inference analysis on the target node pH, the physical variable flow, and the chemical node turbidity. Overall, the results reveal the direct relationship of pH with flow rate, with a nearly 9-fold increase in the probabilities of pH < 5.5 occurring if the water flow rate was high (x $>$ $\stackrel{-}{x}$ + 2σ). Furthermore, as it is consistent, a potential change in water turbidity regimes was detected when the flow rate increases above 17.22 m³/s.

The functional analysis approach studies the whole process and extracts continuous information from the dataset, which translates into a much smaller loss of information and a reliable study of the different trends in the data. Additionally, the smoothing process, which transforms the discrete data into functions, achieves a precise approximation of the original data points enhancing the elimination of erroneous values, such as failures in the data-gathering sensors.

The results of the functional approach automatically identify a series of outliers in each variable. A time-dependency analysis was implemented to assess the reliability of the detected outliers. This process reveals that, on average, several variables present 31% outliers on the same dates. Additionally, studying the values of those anomalies in each correlated variable reveals hidden correlations between them. Most interesting is the case of water flow, rainfall, and turbidity, which are inversely correlated with the pH and present. On the other hand, conductivity shows a direct relationship with the pH, which seems to be a particular behavior of this fluvial system. Consequently, these results validate the functional approach for the automatic detection of outliers and the modeling of a fluvial system based on the time series data collected by the SAIH and SAICA water control networks.

5. Conclusions

The occurrence of acid mine drainage events in the Cúa River (NW part of Spain) was analyzed through statistical process control, Bayesian machine learning, and functional data analysis. The data studied were recorded by a water control station in the Cúa River with daily records of conductivity, dissolved oxygen, pH, rainfall, temperature, turbidity, and water flow between January 2011 and August 2022.

Statistical process control acts as a first approach to the problem. This method takes into account the time correlations within the data due to the use of rational subgroups, but this implies the loss of information. Adding the lack of normality, and high variability of the data, the results end up being plagued with false alarms. Nevertheless, the $\stackrel{-}{x}$ chart provides valuable information regarding the variability of the parameters that control water quality in the Cúa River, as well as the different standard deviation intervals for each variable, which are leveraged by the Bayesian analysis.

The Bayesian model can achieve a higher degree of understanding of the behavior of the fluvial system. A valuable measure for evaluating the effectiveness of the JPD representation among the monitored variables, known as CTF, was calculated. The resulting CTF of 74.581% means that the model provides an accurate representation of the behavior of the aquatic system. In order to know in which potential ranges of values the pH will be found, the calculation of the RMBI between pH states and their predictor variables revealed the most important variables to watch out for. The results showed that for pH ranges below the lower limit set in the RD 817/2015, the potential predictors were flow, month, and turbidity of the water.

Lastly, the functional approach is not dependent on the statistical distribution of the data and studies it as a complete time unit, which leads to better trend analysis, a reduction in the information lost in the process, and more refined outlier detection. The validation of the functional results was carried out through a multivariate temporal analysis of the detected anomalies. Several correlations were found by taking a closer look at the time dependency between the outliers of the different variables. In particular, the variable rainfall presents a direct relationship with the increase of the water flow and the turbidity of the fluvial system, as well as the decrease in the pH and conductivity levels. Additionally, the temporal analysis showed that in some cases, the pH levels do not fall right at the beginning of a precipitation event; rather, they decrease typically 1.5 later and remain low for 4.5 days on average despite the precipitation event having been completed.

It is worth noting that the methods used adapt effectively to the data collected continuously by the hydrological and water quality monitoring networks installed in the study area. This makes it possible to analyze and comprehend the specific behavior of a river over an interval of time, as well as draw inference scenarios and automatically detect reliable multivariate outliers between the different variables studied.

The reliability of the methods implemented is dependent on the availability of time series data without significant portions of missing values. As with most statistical methods, the absence of continuity can translate into a negative impact on the accuracy of the models; therefore, the obtention of faithful results is heavily dependent on the adequate maintenance of the monitoring networks. The national authorities in charge of the protection of water quality in this river basin are planning the development of specific actions with the intent to treat water in the event of acidic conditions. Thus, this study may be an interesting and powerful tool to identify anomalous events and to establish and determine the requirements of future treatments.