We can define a database of n fish farms as ℱ = { F 1 , F 2 , … , F n } . Each fish farm F ∈ ℱ contains several cages   F = { C 1 , … , C m F } , where m_F is the number of cages in the farm. Finally, each cage C ∈ F can be represented as a time series   C = { c 1 , c 2 , … ,   c t C } , where c i ∈ C is the observation of cage C at time i and t C is the number of observations in the cage. Each observation c ∈ C contains information about the oxygen concentration in the cage (in   m g L − 1 ). Effectively, each cage consists of a univariate time series. In the original dataset, other information was available, such as temperature or cage depth. However, in the experiments, variables other than oxygen concentration were not found to improve the hypoxia detection accuracy of the model. Therefore, these were excluded from the model.

The case study comprises data collected from 14 fish farms located in different places across eastern Canada. The data ranges over five years, from 2017 to 2023. The sampling period is detailed in Table 1 . The observations of each time series are captured with a high but, due to operator stocking plans, irregular sampling frequency. The period between consecutive observations can range from 1 to 3 minutes. We aim to forecast hourly values of oxygen concentration. Therefore, we aggregated the data in each cage to an hourly granularity. This is accomplished by taking the median of the collected values of each hour. If no observation was captured in a given hour this represented a missing observation. As an example, Figure 1 illustrates a sample of the time series of oxygen concentration in a given cage.

Figure 2 shows the number of data points over each month by fish farm, from October 2017 to December 2022. Until 2019, the data came from only two farms (SI and LP). The data for the year 2022 comprises the most data points, while there is a period with very little data during 2018.

We carried out an extensive exploratory data analysis to draw insights about the oxygen concentration dynamics. We illustrate the main findings in this section. Figure 3 shows the oxygen concentration over time at each farm. For visualization purposes, we aggregated the data in each farm by taking the median across all available cages. Overall, this plot indicates the presence of yearly seasonality influencing oxygen concentration.

There are some periods of missing data, especially in the years 2018 and 2019. The mechanisms behind the missingness are either missing completely at random (e.g. random sensor malfunction in a given observation) or missing at random (e.g. sensor shutdown during some period for maintenance), according to the definitions described by Schafer and Graham (2002).

We explored the distribution of oxygen concentration using violin plots in Figure 4 . For this analysis, we studied each farm as a whole and disregarded information about individual cages. Overall, the distribution varies across farms. Besides the different locations of each farm, one possible reason for this variability was the sampling period at each farm is different, covering various lengths of time and capturing data during different seasons which is known to influence oxygen fluctuations (Burke et al., 2021).

This work is devoted to the prediction of hypoxia events, which can be defined as an observation with an oxygen concentration below 7 mgL^{− 1} (Oppedal et al., 2011). In this context, we analyzed the distribution of the oxygen condition (normal or hypoxia) across fish farms and over different months. Figure 5 shows the relative frequency (%) of each condition across all observations at each fish farm. Hypoxia is generally a rare occurrence across all farms; however, events are more common at some farms such as DH and SI compared to others. Figure 6 also shows that hypoxia is more common during the summer and fall seasons. This is expected as these events are correlated with higher water temperatures, which reduces the solubility of oxygen in seawater (Burke et al., 2021).

The general goal behind forecasting is to predict the value of the upcoming observations of oxygen concentration in a given cage, c t C + 1 , … , c t C + h , given the historical data, where h denotes the forecasting horizon.

We formalize the problem of time series forecasting based on an auto-regressive strategy. Accordingly, observations of a time series are modeled based on their recent lags. More precisely, the value of c_i is modeled based on the most recent q timesteps: X i = { c i − 1 , c i − 2 , … , c i − p } , where c i is the observation we want to predict and X i ∈ X ⊂ ℝ p represents the corresponding input lagged features (Bontempi et al., 2013). This approach leads to a multiple regression problem where the temporal dependency is modeled by using past observations as explanatory variables. Using h = 1 for illustration purposes, the goal is to transform each time series into a matrix structure as follows:

Each row in the matrix above is a sample (X,c). The last column of the matrix denotes the target variable (future observations) and the remaining columns represent the lagged explanatory variables (X). D i denotes the reconstructed time series in a matrix structure for a given cage C i . For each forecasting horizon   h , the goal is to build a multiple regression model f that can be written as c i + h = f ( X i , Z i ) , where Z i represents additional covariates known at the i-th instance. These covariates can be, for example, explanatory time series related to the phenomenon. For example, the season of the year.

Traditionally, forecasting models are trained using the historical data of the time series of interest, often referred to as local models (Januschowski et al., 2020). In our case, we need to forecast time series from several cages that belong to multiple fish farms. Leveraging the historical observations of all available time series can be valuable to build a forecasting model. The dynamics of the time series are often related, and a model may be able to learn useful patterns from some time series that have not revealed themselves in others. Models that are trained on multiple time series are often referred to as global forecasting models (Godahewa et al., 2021).

The popularity of global forecasting methods has surged following their successful use in the M4 forecasting competition, which featured 100000 time series from different application domains. The contest was won by a global forecasting method developed by Smyl (2020) that combines an exponential smoothing method with a recurrent neural network. Since then, several other works have been developed for training forecasting models with multiple time series (Godahewa et al., 2021).

The key motivation for our use of a global forecasting approach is giving the learning algorithm access to additional data. Machine learning algorithms tend to perform better with larger training sets, especially those with many parameters (Cerqueira et al., 2022). We will build a global model to forecast oxygen concentration based on the available information across all cages from all fish farms.

Conveying the uncertainty around forecasts is critical for better decision-making. Uncertainty in numerical forecasts is usually quantified by predicting intervals (Khosravi et al., 2011), or probabilities (Gneiting and Katzfuss, 2014).

For binary events, we can simply output the probability of that event occurring. An instance of a binary event is exceedance, which, as mentioned before, refers to when a time series exceeds a predefined threshold.

The probability of a hypoxia event h_i denotes the probability that the oxygen concentration c_i falls below 7 mgL^{− 1} in any given instant i. In this type of problem, h_i is usually modeled by resorting to a binary target variable b_i , which can be formalized as shown in Equation 1 :

This setting leads to a classification problem where the goal is to build a classifier g of the form b_i = g(X_i ). Logistic regression is commonly used for this purpose (Taylor and Yu, 2016). However, in our problem, a classification model would be unable to output the numeric oxygen forecasts which can be important for on-farm decision-making. Accordingly, in Section 4, we present an alternative approach to standard classification.

Understanding and forecasting water quality and oxygen dynamics is a relevant task and many works have been devoted to it (Zhang et al., 2022; Hai et al., 2023; Zhao et al., 2023). Liu et al. (2021) developed an ensemble based on Bayesian model averaging using a hybrid approach that includes artificial neural networks. Therein, the authors outline several state-of-the-art approaches to forecasting dissolved oxygen in the context of aquaculture.

Several works address the problem of detecting hypoxic events using data-driven methods. Arepalli and Naik (2024) tackle this problem using a deep learning method based on self-attention and LSTM layers. Besides dissolved oxygen, they also use other input variables such as temperature or turbidity. They report 99.8% accuracy, though this metric is usually a poor choice for problems involving imbalanced distributions (Bronco et al., 2016). Similar to us, Politicos et al. (2021) modeled hypoxia using tree-based models and reported an F1 score ranging between 0.84 and 0.93. The F1 score is a classification metric that combines precision and recall measures. However, while we deal with North Atlantic ocean data, the data used in their study was collected from an enclosed lagoon in the Greek Mediterranean.

We aim to predict hypoxic events and forecast oxygen concentration using a single model, whereas previous research addresses only one of these problems. As mentioned before, coupling forecasts with hypoxia probabilities can be useful to foster more informed decision-making by farmers.

The development of the model works in three main steps:

In the next subsections, we explain each step in turn.

In this section, we detail the performance estimation procedure and evaluation metrics used.

For illustration purposes, the model was tested on a sample cage in the LP fish farm and its predictions are shown in Figures 10 , 11 . Figure 10 illustrates the forecasts for the next 24 hours, along with the actual and previous values of the series. In this particular example, the model can match the actual values closely.

Figure 11 shows the estimated probability of a hypoxic event in an example sample period. In it, a hypoxia event occurs in the final part of the sample (when the red line equals 1), and the model can detect it.

This section presents a quantitative evaluation of the model described in Section 4.

We start by comparing the developed method that leverages forecasts into exceedance probability estimates using the CDF with a probabilistic binary classification model according to the log loss. The log loss evaluates how close the probability estimates are to the actual values (0 or 1). Having accurate probability estimates is important so that farmers can assess and react to a situation in a reliable way. The classification model is trained using a lightgbm using the same procedure we used for training the forecasting model (c.f. Section 4.1.2). The only change is that the classifier is optimized for binary classification rather than regression.

The results are shown in Figure 12 . Overall, the developed model (denoted as CDF) shows better performance across almost all fish farms. We conducted a Wilcoxon signed rank test, recommended by Demšar (2006), which confirmed that the difference in performance is significant.

Figure 13 shows the results for each farm in terms of event recall and false alarm rate, along with the number of hypoxia events that occurred during the testing period. The results vary considerably across fish farms. In some cases, the results are not representative due to the low number of events.

(e.g., RB farm). For the farms with a higher number of events (e.g., LP, RI, and SI), the event recall ranges from 60% to 85%. The false alarm rate also varies across farms. For example, in the SI farm, the score is about 0.05. This means that we can expect around 5 false alarm events (sequences of false alarms) every 100 days. High scores in terms of ER are associated with either a low number of events (which means the score is not representative) or a high FAR score. The end-users can control this trade-off using the decision threshold, which in these experiments is set to 0.5. Notwithstanding, this parameter can be optimized using cross-validation.

As mentioned before, the model is designed to make predictions with a lead time of 12 hours. However, sometimes the model might detect events with some delay. Figure 14 shows the average anticipation time (in hours) for each fish farm. The results indicate that the model detects events with an anticipation time that ranges from about 5 to 8 hours.

Finally, we also evaluated the numeric forecasts produced by the model. Figure 15 shows the distribution of absolute error (i.e., the absolute difference between the forecast and actual values) across the forecasting horizon. The median absolute error varies across farms. For example, this value is about 0.25 for the AD farm. Absolute errors above 0.5 are usually outliers (depending on the farm) according to the boxplot.

Pre-processing the data through differencing was found to be an important step in obtaining the reported forecasting performance, however, stabilizing the variance and modeling seasonal effects did not have a significant impact. We remark that, during the development of our model, we attempted to add extra meteorological and oceanographic variables. These included wind speed and direction, air temperature, solar radiation, or upwelling index. However, these did not improve forecasting performance. Another important aspect is that we average the values collected within each hour by taking the median. We adopted this approach to avoid potential anomalies due to sensor errors. However, we could also potentially have missed some brief hypoxia events that occur within each hour.

The learning algorithm applied in this work is a lightgbm that was optimized using random search. During the development stage, we attempted other algorithms, such as linear models and different types of neural networks. Overall, the lightgbm presented the best forecasting accuracy according to the metrics listed.

An important aspect of our work is that we evaluate our model based on realistic assumptions regarding its potential application. Usually, hypoxia detection models are evaluated using standard classification metrics such as precision, recall, and F1-score (e.g (Arepalli and Naik, 2024)). However, these metrics are not suitable for this problem due to the time-dependency among observations (c.f. Section 4.2.2).

The results suggest that the CDF-based approach for computing exceedance probability estimates outperforms a binary classification model. Beyond performance differences, it is important to note that the classification model is unable to produce forecasts, which are desirable in our application.

We could build two models, one regression model for forecasting and one classification model for exceedance probability estimation. However, two different models are susceptible to inconsistencies. For example, instances where the forecasts suggest a high chance of hypoxia but the classifier outputs a small probability for this event. In effect, being able to serve the two types of predictions using a single value provides consistency in the output. This can be important for a reliable use of the predictions by end-users.

There are several ways the model can be improved. We assume a constant variance when estimating the event probability using the CDF. Specifically, at each time step, we evaluate the CDF using the standard deviation of oxygen concentration of the corresponding cage computed with the whole training set. However, this statistic may vary over time, for example, due to seasonal effects. Thus, an adaptive approach (e.g., a rolling standard deviation) may lead to better probabilistic estimates.

Another point of improvement may be a better data compilation process. We apply a global approach by merging the training observations from all fish farms in a single dataset. This approach has several advantages such as the increased training sample size. This feature can be important if the user wants to deploy the model in fish farms with fewer data points available. However, introducing data from other farms may not bring value in some cases due to the idiosyncrasies of the farm. We can improve the global-local trade-off by, for example, clustering fish farms before modeling (Bandara et al., 2020).

Finally, the proposed model is based on auto-regression. We can also carry out an extended feature engineering process to improve the representation of the dataset. For example, previous works in the literature (e.g (Liu et al., 2021)) report that feature engineering based on empirical mode decomposition improves oxygen concentration forecasting accuracy.