Probabilistic models for harmful algae: application to the Norwegian coast

Abstract

We have developed probabilistic models to estimate the likelihood of harmful algae presence and outbreaks along the Norwegian coast, which can help optimization of the national monitoring program and the planning of mitigation actions. We employ support vector machines to calibrate probabilistic models for estimating the presence and harmful abundance (HA) of eight toxic algae found along the Norwegian coast, including Alexandrium spp., Alexandrium tamarense, Dinophysis acuta, Dinophysis acuminata, Dinophysis norvegica, Pseudo-nitzschia spp., Protoceratium reticulatum, and Azadinium spinosum. The inputs are sea surface temperature, photosynthetically active radiation, mixed layer depth, and sea surface salinity. The probabilistic models are trained with data from 2006 to 2013 and tested with data from 2014 to 2019. The presence models demonstrate good statistical performance across all taxa, with R (observed presence frequency vs. predicted probability) ranging from 0.69 to 0.98 and root mean squared error ranging from 0.84% to 7.84%. Predicting the probability of HA is more challenging, and the HA models only reach skill with four taxa (Alexandrium spp., A. tamarense, D. acuta, and A. spinosum). There are large differences in seasonal and geographical variability and sensitivity to the model input of different taxa, which are presented and discussed. The models estimate geographical regions and periods with relatively higher risk of toxic species presence and HA, and might optimize the harmful algae monitoring. The method can be extended to other regions as it relies only on remote sensing and model data as input and running national programs of toxic algae monitoring.

Impact Statement

Some algae produce toxins that contaminate shellfish and poison humans upon consumption, posing significant risks to human health. They also impact the aquaculture business since shellfish farmers cannot sell their products once the shellfish are contaminated by toxins. Assessing toxic algae risks can allow farmers to make better-informed decisions and reduce economic loss. We calibrate machine learning models that assess the probability of the presence and hazardous levels of eight toxic algae in Norway based on environmental factors, such as temperature, light, salinity, and mixed layer depth. The probabilistic models can be fed with predictions of environmental factors, and thus also predict the probability of toxic algae. Public agencies can use this information to execute mitigation actions against toxic algae.

1. Introduction

Toxic algae pose a significant threat to human health due to their ability to contaminate shellfish with toxins, resulting in poisoning symptoms or, in extreme cases, death upon consumption. Prevention of poisoning outbreaks relies on monitoring programs to assess the abundance of harmful algae and toxins in shellfish farms. When harmful abundances (HAs) are detected, shellfish sales are forbidden and public advice is provided. Nonetheless, both monitoring programs and preventing the consumption of contaminated shellfish incur significant costs, and there still remains a small possibility of outbreaks that can harm the public (Castberg et al., 2004; Hoagland et al., 2002; Hoagland and Scatasta, 2006; Karlson et al., 2021; Martino et al., 2020; Pettersson and Pozdnyakov, 2013). Assessing geographical regions and time periods with an elevated probability of toxic species detection can offer several advantages, such as optimizing monitoring programs by redistributing resources and efforts, enhancing the protection of public health, enabling early harvesting prior to toxin outbreaks, improving business planning and investment decisions, and fostering increased consumer confidence (Jin et al., 2020). Probabilistic models can serve this purpose and inform likelihood changes based on external factors that influence the growth of toxic algae, such as sea surface temperature (SST), mixed layer depth (MLD), photosynthetic active radiation (PAR), and sea surface salinity (SSS) (Anderson et al., 2011, 2012; Bates et al., 2018; García-Portela et al., 2018; Jauffrais et al., 2013; Kim et al., 2008; Klemm et al., 2022; Paz et al., 2006; Reguera et al., 2012).

Few advancements have been made in the development of probabilistic models for toxic algae. For example, at the USA West Coast, logistic generalized linear models (GLMs) were utilized to model the probability of Pseudo-nitzschia spp. exceeding 10,000 CellsL⁻¹, incorporating observations of ocean color, temperature, and salinity (Anderson et al., 2011). In French Mediterranean lagoons, decision tree rules based on temperature and salinity thresholds were employed to model the probabilities of hazardous abundances of Alexandrium tamarense and Dinophysis spp. (Bouquet et al., 2022). In Irish coastal waters, gradient boosting models (GBMs) were applied to model the probability of the presence/absence of A. tamarense using inputs such as temperature, salinity, and a water stratification index (Klemm et al., 2022). Although these studies have made significant contributions, none have addressed a wide range of toxic algae taxa, nor have they been calibrated for the Norwegian coastal shelf.

Along the Norwegian coast, the Norwegian Food Safety Authority (NFSA) regularly monitors several toxic algae taxa, including Dinophysis acuminata Claparède & Lachmann 1859, Dinophysis acuta Ehrenberg 1839, Dinophysis norvegica Claparède & Lachmann 1859, Alexandrium Halim, 1960, A. tamarense (Lebour) Balech 1995, Pseudo-nitzschia H. Peragallo, 1900, Azadinium spinosum Elbrächter & Tillmann 2009, and Protoceratium reticulatum (Claparède & Lachmann) Bütschli 1885. These taxa are monitored because of their association with Diarrhetic Shellfish Toxins, Paralytic Shellfish Toxins, Amnesic Shellfish Toxins, Azaspiracid Shellfish Toxins, and yessotoxins. The monitoring program has been in operation since 2006 and carries out extensive weekly monitoring of algae abundance and monthly monitoring of toxins in all shellfish farms in operation along the Norwegian coast. This long-term data series makes possible a robust calibration and validation of probabilistic models encompassing a wide range of toxic algae taxa. Moreover, the Norwegian coastline stretches across a high range of latitudes (58°N–71°N) and exhibits substantially different environmental settings, making it an optimal natural laboratory for studying harmful algae responses to a wide range of environmental conditions (Wells et al., 2020).

This study calibrates and validates models for estimating the probability of presence and HA of key toxic algae species, including D. acuminata, D. acuta, D. norvegica, Alexandrium spp., A. tamarense, Pseudo-nitzschia spp., P. reticulatum, and A. spinosum, along the Norwegian coast. The probabilistic models are based on support vector machines (SVMs), which among the many existing machine learning methods is well suited for harmful algae bloom (HAB) modeling as it requires a small amount of training data to find an optimal solution (Silva et al., 2023; Cruz et al., 2021). Since algae in situ data are scarce, this trait is optimal for modeling presence and HA. Some studies have already employed SVM for HAB modeling and showed superior performance against other methods such as artificial neural networks and random forest (Li et al., 2014; Ribeiro and Torgo, 2008). As model inputs, we use SST and PAR estimated from satellite observations and MLD and SSS from operational ocean reanalysis. These inputs are chosen for two reasons: (i) they have enough time span and spatial coverage matching the algae observations in the shellfish farms and (ii) algae growth is mostly driven by these factors. Algae show distinct temperature-related traits that cause them to grow or diminish in different temperature ranges (Basti et al., 2018; Fehling et al., 2004; Guerrini et al., 2007; Nagai et al., 2004; Rial et al., 2023; Röder et al., 2012; Thomas et al., 2012); PAR corresponds to the light available for photosynthesis and therefore strongly influences algae growth (Bill et al., 2012; García-Portela et al., 2018; Jauffrais et al., 2013); salinity variations affect algae by inducing osmotic stress, creating ion stress through the unavoidable absorption or loss of ions, and altering the cellular ionic ratios due to selective mechanisms (Jauffrais et al., 2013; Kirst, 1990; Klemm et al., 2022; Nagai et al., 2004; Rial et al., 2023; Weber et al., 2021); shallower MLD—a common proxy to well-stratified waters—is commonly associated with HABs (Klemm et al., 2022; Reguera et al., 2012). Note that other important variables, such as nutrients, are not included as no product available matches the farms in a long time series. By using those inputs, the SVM skill in modeling each toxic algae probability is evaluated. The models are validated using an operational setting, trained with data spanning from 2006 to 2013, and tested with data covering the period from 2014 to 2019. The influence of the input predictors is assessed for each taxon’s probabilistic model. Seasonal probability and annual risk maps for all targeted algae species are presented for the Norwegian coastline.

2. Material and methods

2.1. Study region

The Norwegian coast is surrounded by the Skagerrak Strait, the North Sea, the Norwegian Sea, and the Barents Sea (Figure 1). Two significant current systems dominate circulation: the Norwegian Atlantic Current (NwAC) and the Norwegian Coastal Current (NCC). The NwAC is an extension of the North Atlantic Current as it flows between the Faroe Islands and Scotland and continues northward along with the Norwegian Continental Shelf break up to the Barents Sea (Eldevik et al., 2009; Furevik et al., 2002). The NCC flows from the south of Norway and along the coast north to the Barents Sea. The NCC is substantially fresher than the NwAC (composed of Atlantic Water), as it transports fresh waters from the land inflow, the Baltic, and the North seas.

Figure 1. Study region. The farm locations are represented by dots, and the circles encompass the area over which the satellite and models are averaged (44 km). Areas used in Figures 9 and 10 are highlighted in red.

The Norwegian coastal waters extend from the sub-Arctic and to the Arctic regions and comprise disparate environmental conditions. In northern Norway, the polar night is from November 18 to January 23, and the sun is above the horizon all day between the May 20 and July 24 (Giesen et al., 2014). SST can vary from 5 °C in winter to 20 °C in summer in the North Sea and from 1 °C to 15 °C in the Barents Sea opening (Chen et al., 2021; Jakowczyk and Stramska, 2014). In the Skagerrak Strait, fresher waters are brought from the Baltic Sea and river runoff, varying from 10 PSU in summer to 30 PSU in winter (Frigstad et al., 2020; Hordoir et al., 2013). In northern Norway, salinity is above 34 PSU most of the time and can decrease to 20 PSU in episodic events of freshwater input during the summer melting season (Frigstad et al., 2020). The MLD can be deeper than 50 m in the winter and shallower than 30 m in the summer due to the input of fresh waters and surface heating (Peralta-Ferriz and Woodgate, 2015).

2.2. In situ data collection

Abundance (CellsL⁻¹) of Alexandrium spp., Alexandrium tamarense group, D. acuta, D. acuminata, D. norvegica, Pseudo-nitzschia spp., P. reticulatum, and A. spinosum, was provided by the monitoring program of algae toxins in mussels and dietetic advice to the public from NFSA. Algae samples are collected at several aquaculture mussel sites weekly (every Monday). The monitoring program is run routinely, and the data used in this study cover from 2006 to 2019. The program sampling method consists of collecting water samples by lowering a tube from the surface to a 3 m depth. A subsample of 25 mL is taken from the water sample and preserved with acidic Lugol’s iodine before being transported to the laboratory for analysis. The subsample (25 mL) is filtered on a membrane filter, and the genus and species present are identified and counted on the whole filter under a light microscope at 200× magnification.

Only 35 shellfish farms in the coastal area, shown in Figure 1, are included in the study because the remote sensing and model data used are less reliable or unavailable (due to spatial resolution) in the inner fjords (see Section 2.3). The amount of data available depends on each farm’s operation time. For 8 locations, we only have 1 year of data, while for 11 sites, we have more than 9 years of data. In situ data are not collected during the winter months—it is out of the productivity season. The total data for training and evaluating the probabilistic models comprise 5919 samples for each taxon. An example of algae abundance time series in Arendal 2019 is shown in Figure 2a-h to illustrate the data used as the input for modeling.

Figure 2. Arendal (in the south of Norway) time series in 2019 as an example of the input data for training and testing the models. The time series for (a) Alexandrium spp., (b) Alexandrium tamarense, (c) D. acuta, (d) D. acuminata, (e) D. norvegica, (f) Pseudo-nitzschia spp., (g) P. reticulatum, (h) A. spinosum, (i) SST, (j) PAR, (k) MLD, and (k) SSS.

2.3. Satellite and model reanalysis data

We use SST (°C) from the ESA SST CCI and C3S global SST reprocessed product level 4, retrieved from the Copernicus Marine Environment Monitoring Service (CMEMS). The product is created by running the Operational Sea Surface Temperature and Sea Ice Analysis system (Good et al., 2020), which combines satellite (AATSR, ATSR, SLSTR, and AVHRR) and in situ observations to produce gap-free maps of daily average SST at 0.05° of spatial resolution (Merchant et al., 2019). In the Nordic seas, the SST uncertainty is below 0.4 °C (Good et al., 2020).

The MLD (m) and SSS (PSU) are provided by the CMEMS Arctic MFC TOPAZ modeling system (Sakov et al., 2012; Xie et al., 2017). The TOPAZ system is a coupled ocean–sea ice model and data assimilation system for the North Atlantic and Arctic Oceans. The ocean model couples a Hybrid Coordinate Ocean Model (Bleck, 2002) with an elasto-viscous-plastic sea ice model (Hunke and Dukowicz, 1997). TOPAZ weekly assimilates available ocean and sea ice data with the ensemble Kalman filter (Evensen, 2003). The MLD is calculated using a density criterion with a 0.01 kg m⁻³, as in Petrenko et al. (2013) and Ferreira et al. (2015). The SSS is extracted from depths of 0–3 m. The TOPAZ product performs well concerning ocean variables near the surface. In the first 200 m, salinity root mean squared difference (RMSD) and bias are below 0.3 PSU and between −0.05 and 0.05 PSU. The temperature RMSD and bias are below 1 °C and between −0.5 and 0.5 °C, respectively (Xie et al., 2017; Lien et al., 2016). Note that MLD uncertainties are not provided, but since it is computed using density, temperature and salinity indicates the quality of MLD estimations. Nevertheless, the MLD error is expected to be about 10 m on the Norwegian coast (L. Bertino, personal communication).

PAR ( $ {\mathrm{Em}}^{-2}{\mathrm{d}}^{-1} $ ) irradiance onto the ocean surface from 2006 to 2019 is retrieved from the GlobColour project, which uses MODIS, SeaWiFS, and VIIRS sensors binned at an 8-day interval at a 4 km of spatial resolution. Originally developed for SeaWiFS, the product presents an accuracy of R = 0.88 and RMSD = 5.7 Em⁻²d⁻¹ compared to in situ measurements (Frouin et al., 2003).

All data are reprojected to stereographic projection centered at 65°N and 7°E and at 4 km spatial resolution using the nearest neighbor interpolation method. Because of the coarse spatial resolution, we extract the location time series as the average of unmasked grid cells within the 11 × 11 grid of the fields around the location. Therefore, the SST, SSS, MLD, and PAR represent the average conditions at the surrounding of the farm and come with a 44 km effective resolution. An example time series of all extracted data is shown in Figure 2i–l for Arendal 2019.

2.4. The calibration of presence probabilistic models

The algae abundance is converted to binary values needed to calibrate the probabilistic models. For the presence models, a threshold of 1 CellsL⁻¹ is used for creating the two classes. This means that samples with abundances above or equal to 1 CellsL⁻¹ are defined as class = 1, and samples with 0 CellsL⁻¹ are defined as class = 0. The models are fed with SST, SSS, MLD, and PAR as input. The models are evaluated at an operational level where past data are used for training (2006–2013) and employed in future data for testing (2014–2019).

The first step is preprocessing the data, starting with scaling (or normalization) of each input predictor:

(1) $$ {\mathbf{x}}_{\mathrm{scaled},\mathrm{pred}}=\frac{{\mathbf{x}}_{\mathrm{pred}}-{\overline{\mathbf{x}}}_{\mathrm{pred}}}{\sigma_{\mathrm{pred}}} $$

where $ {\mathbf{x}}_{\mathrm{scaled},\mathrm{pred}} $ is the scaled sample $ \mathbf{x} $ for an input predictor (pred), $ {\mathbf{x}}_{\mathrm{pred}} $ is the input predictor value, $ {\overline{\mathbf{x}}}_{\mathrm{pred}} $ is the average, and $ {\sigma}_{\mathrm{pred}} $ is the standard deviation computed on the training data (2006–2013). Since the scaled environmental data are significantly correlated, a second preprocessing step using the principal component analysis is used to convert the scaled predictors to four decorrelated principal components (c1, c2, c3, and c4), which are fed into the SVM model.

The SVM probabilistic algorithm consists on first calibrating a margin hyperplane in an n-dimensional space to separate two classes by using the nearest support vectors (samples) of both classes. The computation of the margin hyperplane depends on the input (principal components) and on the SVM hyperparameters that need to be fixed before optimization, such as the kernel function, the penalty factor ( $ C $ ), and the $ \gamma $ (for nonlinear kernels). The kernel function converts the input values to a feature space where the hyperplane is computed. Common kernel functions are linear, polynomial, and radial basis function (RBF). The margin hyperplane is optimized by the Hinge Loss function with a tolerance stop criterion of 0.001. The hyperparameter $ C $ controls the trade-off between maximizing the margin and minimizing the training errors, and $ \gamma $ controls the weight a single sample has on adjusting the hyperplane in an RBF kernel. A deeper explanation of SVM can be found in Cortes and Vapnik (1995), Platt (1999), Tan et al. (2008). The hyperparameters are tuned in a grid search using cross-validation procedure (Hastie et al., 2009) in the training dataset with two folds randomly split 100 times. The grid search uses the correlation between observed presence frequency and estimated probability—obtained in a reliability diagram—as a decision criterion (see Section 2.7). The tuned hyperparameters for all algae models were kernel = RBF, $ C=1 $ , and $ \gamma $ defined by the following equation:

(2) $$ \gamma =\frac{1}{n_{\mathrm{comp}}\times {\sigma}_{\mathrm{all}}^2} $$

where $ {n}_{\mathrm{comp}} $ is the number of principal components ( $ {n}_{\mathrm{comp}}=4 $ ), and $ {\sigma}_{\mathrm{all}}^2 $ is the scalar variance of all the components stacked as one vector. Since the number of samples with the harmful algae taxa detected can represent only a portion of the total dataset, we adjust the weight of the samples—used during the hyperplane optimization—inversely proportional to the class frequencies:

(3) $$ w\left(\mathbf{x}\right)=\frac{n_{\mathrm{sam}}}{2\times {n}_{\mathrm{sam},\mathrm{class}}} $$

where $ w $ is the weight given to sample $ \mathbf{x} $ , $ {n}_{\mathrm{sam}} $ is the total number of samples, and $ {n}_{\mathrm{sam},\mathrm{class}} $ is the number of samples for the class of $ \mathbf{x} $ .

Finally, the second step of the SVM probabilistic algorithm is to use the margin hyperplane to fit a probability function using the Platt (1999) method:

(4) $$ P\left(c\left(\mathbf{x}\right)=1|f\right)=\frac{1}{1+\exp \left( Af+B\right)} $$

where $ P $ is the probability of sample $ \mathbf{x} $ being class 1 ( $ c\left(\mathbf{x}\right)=1 $ ), the input $ f $ is the SVM output of each predicted sample corresponding to its orthogonal distance from the hyperplane, scaled proportionally from $ - $ 1 to 1 defined between the support vectors distance, and $ A $ and $ B $ are the parameters fitted using the maximum likelihood in the training dataset. The SVM is implemented in the Python programming language on Scikit-learn package (Pedregosa et al., 2011).

2.5. Model and data input uncertainties

Model uncertainty is estimated by training new models using two-thirds of the pseudo-randomly subsampled data from the training dataset. Subsampling is repeated 100 times, and all models are employed in the testing dataset. The reliability (see Section 2.7) is then estimated for all interactions and subtracted by their median.

Uncertainties in the input data from model reanalysis and remote sensing data are susceptible to yield errors that can deteriorate the outcome probability estimates. We assess the input uncertainty by adding a random perturbation to the input testing dataset and evaluating the impact on reliability. We generate a perturbed ensemble of input, as follows:

(5) $$ {\mathbf{x}}_i=\mathbf{x}+{\varepsilon}_i $$

where $ \mathbf{x} $ is the input vector, and $ \epsilon $ is the Gaussian white noise $ \in {\mathrm{\mathbb{R}}}^{n\times N}\approx \mathrm{\mathbb{N}}\left(o,{\sigma}^2\right) $ , with n being the size of the input vector and $ \sigma $ is the standard deviation of the input data. We set the $ \sigma $ as 0.4 ° C for SST, 10 m for MLD, 0.3 PSU for SSS, and 5.7 Em⁻²d⁻¹ for PAR. These $ \sigma $ values relate to the reported input errors (see Section 2.3). This experiment is repeated 100 times and the reliability changes are estimated for each interaction.

2.6. Calibration of probabilistic models for higher sanitary thresholds and HAs

The presence of a toxic species is not necessarily harmful to the environment or to human shellfish consumption. The NFSA establishes specific sanitary thresholds of taxa abundance to consider harmful (Table 1). Probabilistic models of HA are more desirable than presence models for shellfish farms because some taxa might be present at a low level without requiring action. We recalibrate the models for an increasing sets of thresholds for each taxon from the presence to HA (with a percentile range of 20%), shown in Table 1. Each threshold is used to define class 0 referring to below the threshold and class 1 above it (see Section 2.4). For example, we recalibrate the A. tamarense probabilistic models for abundances above 1 (presence model), 40, 80, 120, 160, and 200 CellsL⁻¹ (HA model). We estimate the reliability change along the percentiles and evaluate the feasibility of moving from the presence probabilistic model to the HA probabilistic model for each taxon. Note that harmful thresholds for Alexandrium spp. and A. spinosum are not defined by NFSA, so we resort to using 200 CellsL⁻¹ for Alexandrium spp. as it is used for A. tamarense, and 3600 CellsL⁻¹ for A. spinosum as it corresponds to the 99% percentile value of the database.

Table 1. Sanitary thresholds used for calibrating the probabilistic models for each taxon; from presence (CellsL⁻¹ > =1) to HA of each taxon

HA levels are determined by the NFSA. Percentage thresholds correspond to the percentile from the presence to HA.

2.7. Reliability assessment

Our models are evaluated by comparing the probability estimated by the models with the observed frequencies for different bins of the probability density function—an approach called reliability diagram (Bröcker and Smith, 2007). First, the estimated probabilities in the testing period (2014–2019) are ranked in 10 separate bins using the percentiles to define the bin widths. Probability values from percentiles 1%–10% are clustered in bin 1, and from 11% and 20% in bin 2, continuing up to bin 10. Then, the bins’ average probability and observed frequency are computed. For example, the presence probability of D. acuta is estimated for all weekly observations in the farm locations available from 2014 to 2019, following the spatial variability and time intervals of the NFSA monitoring program. We match the estimated presence probability with the observed presence (0 or 1), which is used to estimate the average presence probability and observed frequency in the bins as aforementioned. As an example, we verify that when our model predicts that there is a 10% chance that a species is detected, it happens 10% of the time in an independent validation period. Hence, our models are evaluated by their skill in modeling the frequency of that toxic algae presence—or HA—being detected in the monitoring. To quantify how well our model is doing, the Pearson correlation (R), Root Mean Squared Error (RMSE), and average bias (AB) between the averaged probability and observed frequency are calculated as follows:

(6) $$ R=\frac{\operatorname{cov}\left({F}_{\mathrm{bin}},{P}_{\mathrm{bin}}\right)}{\sigma_{{\mathrm{F}}_{\mathrm{bin}}}{\sigma}_{{\mathrm{P}}_{\mathrm{bin}}}} $$

(7) $$ \mathrm{RMSE}=\sqrt{\frac{1}{n_{\mathrm{bin}\mathrm{s}}}\sum \limits_{i=1}^{n_{\mathrm{bin}\mathrm{s}}}{\left({P}_{\mathrm{bin},i}-{F}_{\mathrm{bin},i}\right)}^2} $$

(8) $$ \mathrm{AB}=\frac{1}{n_{\mathrm{bin}\mathrm{s}}}\sum \limits_{i=1}^{n_{\mathrm{bin}\mathrm{s}}}\left({P}_{\mathrm{bin},i}-{F}_{\mathrm{bin},i}\right) $$

where $ {F}_{\mathrm{bin}} $ and $ {P}_{\mathrm{bin}} $ are the observed frequencies and averaged estimated probabilities of the bins, cov is the covariance function, $ \sigma $ is the standard deviation, $ {n}_{\mathrm{bins}} $ is the number of bins (10), and i refers to one bin. For evaluating the RMSE changes along different sanitary thresholds, the RMSE is normalized by the range of minimum and maximum probabilities the model can predict:

(9) $$ {\mathrm{RMSE}}_{\mathrm{norm}}=\frac{\mathrm{RMSE}}{P_{\mathrm{max}}-{P}_{\mathrm{min}}} $$

where $ {\mathrm{RMSE}}_{\mathrm{norm}} $ is the normalized RMSE, $ {P}_{\mathrm{max}} $ and $ {P}_{\mathrm{min}} $ are the maximum and minimum probabilities the model could estimate in the testing dataset.

2.8. Presence model sensitivity to predictors

Machine learning is often mistakenly perceived as a “black box” approach that provides outputs without revealing the underlying inference processes. While it may not provide explicit equations describing the inference processes, the models can be used to understand how the target responds to input variations. One of such method is estimating the sensitivity of the model, which should not be confused with the commonly referred term “sensitivity” associated with true positive rates or recall. In the sensitivity analysis, we fix the values of all input predictors except one, for which we span the variability and observe the response changes in the model’s output. In our study, we choose to fix three input predictors at their respective medians and vary the remaining one from the 2.5th percentile to the 97.5th percentile based on the total dataset. For example, when estimating the response of each taxon to SST, we utilize the medians of MLD, SSS, and PAR. We then vary the SST values from 3.7 °C to 18.8 °C and compute the corresponding probability response. Specifically, the medians for SST, MLD, SSS, and PAR are 11.96 °C, 9.27 m, 32.93 PSU, and 29.58 Em⁻²d⁻¹, respectively. It is important to note that the sensitivity is analyzed by the fixed values of the other predictors (their median), which can influence the amplitude of the response. Nevertheless, this approach allows us to evaluate the marginal response to variations in the input across all ranges of variability.

2.9. Presence probability maps

We generate presence probability maps for all the analyzed taxa in this study. To ensure consistency between the calibrated model and the probability maps, we employ an 11 × 11 average filter to the input, ensuring that the spatial variability present in the maps represents the one which the models are calibrated. In addition, we exclude grid cells located beyond a 30 km distance from the coastline, as the models are calibrated solely using coastal data, and their application is limited to coastal areas. Subsequently, we apply the probabilistic models on a weekly basis for all months and years, yielding what we refer to as the weekly probability maps. These weekly maps are then averaged for each year, resulting in an annual average of weekly probabilities, and then are averaged from 2006 to 2019.

2.10. Seasonal probability estimation

We select four distinct sampling locations, namely Arendal, Bømlo, Nærøy, and Vesterålen, to extract regional averaged seasonal time series from 2014 to 2019. While some regions look similar (e.g., Nærøy and Vesterålen), there is a slight delay as a result of the latitude difference. These regions are specifically chosen for their geographical spread, covering southern to northern Norway, and for having the longest available time series among the nearby shellfish farms. To derive the seasonal probability estimates, we calculate the average probability for each week of the year across all years available in the testing dataset. This approach allows us to capture the typical patterns of presence probability for all taxa under consideration. Additionally, for four taxa where the HA probability models demonstrate reliable performance, we also compute the seasonal HA probability. HA models with inadequate performance (see Section 3.1) are excluded from the computation of seasonal probability, ensuring that only models with satisfactory results are utilized in determining the overall seasonal probability estimates.

3. Results

3.1. Reliability of the probabilistic models

The presence probability estimates of all taxa are strongly correlated with the observation frequency (Figure 3, Table 2). The R varies from 0.69 to 0.98, the RMSE from 0.8% to 7.8%, and the AB from −5.2% to 5.7%. It is noteworthy that the probabilistic models estimate distinct percentage ranges for each individual taxon. The models for A. tamarense, D. acuta, P. reticulatum, and A. spinosum can estimate low chance of presence (P < 1%), while models for the Alexandrium spp., D. acuminata, D. norvegica, and Pseudo-nitzschia spp. cannot. None of the models can estimate high probabilities of presence—with a maximum probability ranging from 14% for A. spinosum to 69% for D. norvegica.

Figure 3. Reliability diagram. Comparison between the estimated presence probability and the observed presence frequency estimated in 10 bins for all taxa and their linear regression.

Table 2. Statistical results for presence models for the eight taxa studied

Values are estimated by comparing the probability average and the observed presence frequency of percentiles bin shown in Figure 3. The minimum and maximum probabilities returned by our model are reported. * Denotes significant correlation tested against an $ \alpha =0.05 $ .

The model uncertainty shows robust R for most taxa models, where subsampling the training data has led to deviations between $ - $ 0.1 and 0.1 (Figure 4a). The exception is P. reticulatum as R deviations are larger than 0.2. The RMSE is consistent to a certain degree among the taxa models (Figure 4b) with deviations between −1 and 1%. The AB follows a similar pattern as its deviations are in the −1 and 1% interval (Figure 4c).

The input data uncertainty is relatively low for most taxa models. The R remains relatively high for most taxa with values above 0.8 independent of input variables with artificial white noise added (Figure 5a). The exception is the P. reticulatum model, showing a general decrease in R down to 0.59. The RMSE is relatively low with values below 8.2% (Figure 5b). A significant increase in RMSE is observed for D. acuminata, being 0.4% higher compared to the reference model. AB shows low changes when white noise is added (Figure 5c). In general, no significant degradation of results is observed when adding white noise that is equivalent to the input error products, except for the P. reticulatum model.

As the sanitary thresholds increase (in % of HA), the quality of the models evolves differently (Figure 6). The models of Alexandrium spp., Alexandrium tamarense, D. acuta, and A. spinosum show little R decrease (remaining significant) and a low increase of $ {\mathrm{RMSE}}_{\mathrm{norm}} $ . The D. norvegica models show a decrease in R toward higher thresholds, where R is still significant up to the 80% percentile (3200 CellsL⁻¹) but decreases to R < 0.5 at HA. The D. norvegica models also show an increase in $ {\mathrm{RMSE}}_{\mathrm{norm}} $ . The models for D. acuminata, Pseudo-nitzschia spp. and P. reticulatum exhibit significant R only in the presence thresholds, and the $ {\mathrm{RMSE}}_{\mathrm{norm}} $ increases rapidly for higher sanitary thresholds. The AB for all taxa models tends to 0% while the thresholds increase (Figure 6c), but due to the range reduction (minimum and maximum predicted, not shown) rather than improvement of the models. It is noteworthy that the total number of observations above the sanitary levels (referred to as class = 1) decreases in higher sanitary thresholds for all taxa (Figure 6d), showing a maximum of 1974 presence observations for the Pseudo-nitzschia spp. and a minimum of 11 HA observations for the P. reticulatum.

Figure 4. Model uncertainty for the presence models. The model uncertainties are shown in R (a), RMSE (b), and AB (c) deviations of the median over 100 interactions of randomly subsampling two-thirds of the training dataset for training new models and applying them to the testing dataset. The x-axis is the model for each taxa.

Figure 5. Data input uncertainty for the presence models. The data input uncertainties are shown in R (a), RMSE (b), and AB (c). The blue bars are the reference models (shown in Table 2), and the orange bars are the average of 100 interactions of randomly adding white noise to the testing input dataset. Black lines are the 95% confidence interval. The x-axis is the model for each taxa.

Figure 6. Statistical changes along different sanitary thresholds. The changes of R (a), RMSE (b), AB (c), and total number of samples above the threshold (d) are shown for different sanitary levels of each taxa. The x-axis shows the relative percentile threshold from presence (CellsL⁻¹ > =1) to the HA of each taxa. The black dashed horizontal line in (a) corresponds to the significant level threshold for p-value < 0.05.

3.2. Presence models response to environmental input

The presence model response of each taxon shows the probability peak at different intervals of SST (Figure 7a). Among the Dinophysis spp., D. norvegica shows a probability increase in colder waters, D. acuminata shows the highest probability between 10 °C and 11 °C, and D. acuta probability peaks around 18 °C. The Alexandrium spp. and Alexandrium tamarense models show similar responses where the highest probabilities of both are in colder waters near 7.5 °C. Pseudo-nitzschia spp. probability increases toward warmer waters. P. reticulatum exhibits the highest probability near 11 °C, and the A. spinosum probability peaks around 14.5 °C.

Figure 7. Presence models sensitivity to (a) SST, (b) MLD, (c) SSS, and (d) PAR. For each sensitivity simulation, the other predictors are a fixed value (the median of the dataset). SST, MLD, SSS, and PAR medians are 11.96 ° C, 9.27 m, 32.93 PSU, and 29.58 Em⁻²d⁻¹.

The models’ responses to MLD display two regimes (Figure 7b). Alexandrium spp., A. tamarense, D. acuminata, P. reticulatum, and A. spinosum exhibit increasing probabilities toward shallower MLD (exponential distribution), with amplitudes varying from 4% to 43%. D. acuta shows a slight increase in probability toward MLD at around 40 m. Similarly, Pseudo-nitzschia spp. probabilities also increase toward shallow MLD, peaking around 27 m. In contrast, D. norvegica shows apparently low sensitivity to MLD, with probabilities varying from 17% to 22% across all MLD values.

The influence of SSS on the models of toxic algae reveals three configurations. Probabilities for Alexandrium spp., A. tamarense, D. acuta, and P. reticulatum are lower than 5% at 25 PSU, and increase to more than 10% at 34.5 PSU. D. acuminata and Pseudo-nitzschia spp. exhibit a sharper increase, with probabilities ranging from 15–20% at 25 PSU to 37–46% at 34.5 PSU. A. spinosum probabilities increase from 0% at 25 PSU to 6% at 34.5 PSU. In contrast, D. acuta shows minimal changes across the simulated SSS range. Finally, D. norvegica displays a sharp decrease in probability, from 47% at 25 PSU to 14% at 34.5 PSU.

Increased PAR demonstrates a positive influence on the probabilities of most taxa. The probabilities of D. acuminata, D. norvegica, and Pseudo-nitzschia spp. increase from values below 24% at 7 Em⁻²d⁻¹ to values above 36% at 50 Em⁻²d⁻¹. Alexandrium spp., A. tamarense, and P. reticulatum also show an increase, although to a lesser extent, from values below 5% at $ 7\hskip0.1em {\mathrm{Em}}^{-2}{\mathrm{d}}^{-1} $ to values above 8% at 50 Em⁻²d⁻¹. D. acuta is the only species showing a substantial decrease in probability with higher values of PAR, from 14% at 7 Em⁻²d⁻¹ to 2% at 50 Em⁻²d⁻¹. Finally, A. spinosum probability shows little response to PAR.

3.3. Presence probability maps

The Alexandrium spp. and Alexandrium tamarense annual average of weekly probabilities have a similar pattern, but the latter reaches lower probability values (Figure 8). The Alexandrium spp. annual average of weekly probabilities along the Norwegian coast varies from 6% to 18%, while the A. tamarense varies from 3% to 12%. For both taxa, probabilities are larger in the northern parts along the Norwegian coast. The probability decreases in the southern regions, and is lowest values in the Skagerrak Strait. The Dinophysis spp. annual average of weekly probabilities is also similar among the three species but show different amplitudes. D. acuta shows the lowest probability values, varying from 2 to 9%, D. acuminata varies from 17% to 29%, and D. norvegica varies from 12% to 35%. For all Dinophysis species, the highest probability is found along the southern Norwegian coast in the Skagerrak, and it decreases northward. Pseudo-nitzschia spp. probabilities are similar from the western to the northern Norwegian coast, but decreases in the Skagerrak Strait, ranging from 17% to 35%. P. reticulatum probabilities are nearly homogeneous along the entire Norwegian coast, with the lowest probability in the Oslofjord and far northern coastal waters, varying from 3% to 7%. Finally, A. spinosum shows the lowest annual average of weekly probabilities, varying from 0% in the Skagerrak Strait to 3% on the western and northern coast.

Figure 8. Spatial distribution of the annual average of the weekly probability (in %) predicted by the presence models for Alexandrium spp. (a), Alexandrium tamarense (b), D. acuta (c), D. acuminata (d), D. norvegica (e), Pseudo-nitzschia spp. (f), P. reticulatum (g), and A. spinosum (h). The prediction period corresponds from 2006 to 2019.

3.4. Seasonal presence and HA probabilities

The presence probability of Alexandrium spp. displays significant seasonal variability in Bømlo, Vesterålen, and Nærøy, peaking from week 10 to 30 (Figure 9). However, in the Arendal location, the seasonal probabilities remain low. A. tamarense also exhibits seasonal variability in all four regions, with a brief period of increasing probability from week 10 to 20 in Arendal, and a longer period from week 1 to 40 in the other regions.

Figure 9. Seasonal variability of all the presence probability of (a) Alexandrium spp., (b) Alexandrium tamarense, (c) D. acuta, (d) D. acuminata, (e) D. norvegica, (f) Pseudo-nitzschia spp., (g) P. reticulatum, and (h) A. spinosum. Seasonal probabilities are shown for the regions of Arendal (blue), Bømlo (orange), Vesterålen (green), and Nærøy (red).

Dinophysis spp. shows distinct seasonal periods of increased presence probability. D. acuta starts increasing from week 20 to 25 but decreases to values close to 0% from week 45 onward. On the other hand, D. acuminata and D. norvegica exhibit probabilities higher than 10% throughout the year, indicating their presence in all seasons, including winter. The timing of increased probability for D. acuminata is delayed as the sampled region is located northward, starting in week 11 in Arendal and week 20 in Vesterålen. The seasonal probability amplitude of D. norvegica varies significantly among the regions. Arendal shows probabilities higher than 50% from week 13 to 30, while the other regions only reach values up to 30% in shorter periods.

In most regions, except Nærøy, the presence probabilities of Pseudo-nitzschia spp. display high-frequency changes and a weak seasonal pattern. Similar to D. acuminata and D. norvegica, the probabilities remain higher than 10% throughout the year, indicating that Pseudo-nitzschia spp. may be detected year-round. P. reticulatum exhibits increased probabilities in short periods, shifting to later in the year as the sampled region moves northward. In Arendal, the window of increased probability is from week 15 to 30, while in Vesterålen, it extends from week 20 to 37. A. spinosum seasonal probability detection is generally low, varying from 0% to 9%. In Arendal, the increased probability period is from week 30 to 53, whereas in the other regions, it is from week 25 to 40–53.

The seasonal HA probabilities of Alexandrium spp., A. tamarense, D. acuta, and A. spinosum are characterized by shorter periods compared to their overall presence probabilities (Figure 10), indicating that while these taxa can be detected throughout the year, the occurrence of HA is limited to shorter time windows. For instance, in Bømlo, the seasonal period with elevated HA probability for A. tamarense (Figure 10f) is between weeks 10 and 30, while the presence probability remains above 0% throughout the entire year. This pattern of HA restriction is consistent across different taxa and in other regions as well. Furthermore, all the HA observed from 2014 to 2019 only occurred in the shortened period of increased HA probabilities.

Figure 10. Seasonal presence (in blue) and HA probabilities (in orange) of Alexandrium spp. (a, b, c, and d); Alexandrium tamarense (e, f, g, and h); D. acuta (i, j, k, and l); and A. spinosum (m, n, o, and p) from 2014 to 2019. Probabilities are shown for the regions of Arendal (column 1: a, e, i, and m); Bømlo (column 2: b, f, j, and n); Næroy (column 3: c, g, k, and o); and Vesterålen (column 4: d, h, l, and p). HA observations by the local monitoring from 2014 to 2019 are shown as red columns.

4. Discussion

4.1. SVM skill in modeling harmful algae probability

Few machine learning techniques for modeling the probability of toxic algae have been explored, including GLM (Anderson et al., 2011), decision trees (Bouquet et al., 2022), and GBM (Klemm et al., 2022). We demonstrate that SVM is a highly reliable approach for estimating the presence and HA probability of toxic algae in Norwegian coastal waters. By employing an RBF kernel, the SVM model is able to learn distinct responses to different environmental inputs. For instance, the response of D. acuminata to SST exhibited a bell-like shape, while its response to MLD followed an exponential pattern. The parameter $ \gamma $ plays a crucial role in controlling the smoothness of these responses, and careful fine-tuning is necessary to avoid overfitting the model to the training dataset. To ensure realistic and smooth responses without overfitting, the $ \gamma $ has to be set as a function of the input variables (Equation 2). Furthermore, most of the taxa models—except for P. reticulatum—demonstrate a low model uncertainty, indicating that the SVM converges to similar solutions despite the randomness effect of subsampling the data for training. Modeling the P. reticulatum presence proves more challenging because of the higher model uncertainty and inferior performance.

It is important to acknowledge that the SVM models calibrated with SST, MLD, SSS, and PAR cannot estimate a 100% probability, because we only consider a few of the very many potential input predictors. Other inputs such as prey availability, nutrients, and grazing (Kim et al., 2008; Smayda, 2008; Wells et al., 2020) may be critical but are not considered here since they are not available for a long enough period. Nevertheless, the SVM algorithm provides the minimum and maximum probability range associated with the available set of inputs. With longer time series of observations in the future, additional input predictors can be incorporated, reliability improved, and the range of predictable probability can be expanded.

The reliability of the probabilistic models is likely impacted by the input quality (see Section 2.3). We estimate the input uncertainty as relatively small—except for P. reticulatum—considering their errors assessed on their original spatial resolution and on a global scale ( $ \sigma $ in equation 5). However, when we average each farm time series in a 44 km resolution, we might increase the input errors that are not included in the estimation of data input uncertainty. For example, episodic freshwater input from small rivers may introduce high spatial variability of SSS near the coast (Frigstad et al., 2020). The SSS response for most taxa is highly variable in the 30–35 PSU range (Figure 7), and high spatial variability may mismatch the observation frequency in the farm (induced by the local SSS) with the probability modeled by SSS averaged at 44 km resolution. Although we have not estimated the effect of coarser resolution on the model skill, it should not be severe since we can still produce reliable models using a 44 km resolution.

When the sanitary thresholds are increased to HA levels for calibrating the probabilistic models, the R and RMSE remain relatively unchanged for Alexandrium spp., A. tamarense, D. acuta, and A. spinosum, while they substantially worsens for D. norvegica, D. acuminata, Pseudo-nitzschia spp., and P. reticulatum. The decrease in the number of samples with class = 1 available for training the SVM model likely contributed to this deterioration. For instance, there are 1586 observations of presence available for training the presence model of D. norvegica (threshold > = 1 CellsL⁻¹), while there are only 13 observations of HA available for the HA model (threshold >4000 CellsL⁻¹) of the same species. The decrease of observations representing positive instances as the threshold increases constrains the patterns that the SVM can learn and, importantly, the remaining data available for evaluating the model’s performance.

The small changes in R and RMSE for Alexandrium spp., A. tamarense, D. acuta models may be attributed to the small difference between the presence threshold (1 CellsL⁻¹) and the HA threshold (200 CellsL⁻¹). For Pseudo-nitzschia spp., D. acuminata, D. norvegica, and P. reticulatum HA models, the number of HA observations is too few to expect good skill. In the coming years, with longer periods of training data, we may be able to improve our estimation for these taxa. Furthermore, the conditions leading to HA are often more complex than those for presence as they depend on compound environmental factors.

4.2. Environmental influence: geographical and seasonal variability

The analysis of the simulated environmental responses provides insights into how each taxon responds to different environmental conditions. However, it is important to note that we provide an analysis of the response to one input predictor using fixed values (median) for the other input predictors. For example, the probability response to SST is observed at a PAR of 29.58 Em⁻²d⁻¹. If we set PAR to 0 Em⁻²d⁻¹, the SST probability response decreases close to zero, but the general variability of the response stays the same (not shown). Our primary interest lies in the relative probability, and how the responses may explain the geographical and seasonal variability of each taxon group.

The simulated responses indicate that the probability of presence for most taxonomic groups increases in proportion to PAR, except for D. acuta and A. spinosum. This observation aligns with previous studies that have established a positive correlation between PAR and the occurrence or growth rates of Alexandrium spp. (Anderson et al., 2012), D. acuminata (Kim et al., 2008), Pseudo-nitzschia spp. (Bates et al., 2018), and P. reticulatum (Paz et al., 2006). Comparative analysis between D. acuminata and D. acuta has revealed that the latter is more susceptible to photodamage under high light intensity (García-Portela et al., 2018). This susceptibility may explain the higher probability of D. acuta in relatively low PAR and its typical bloom period in autumn (Dahl and Johannessen, 2001; Naustvoll et al., 2012). A. spinosum growth rates measured in laboratory are low sensitive to light intensities (Jauffrais et al., 2013), which may contribute to its low probability response to PAR changes. It is worth noting that D. acuta and A. spinosum cannot grow in complete darkness like most other taxonomic groups, resulting in close-to-zero presence probabilities during winter months (Figures 9 and 10). Therefore, PAR may limit the presence of toxic algae during winter and contribute to the restriction of occurrence of D. acuta during the autumn season. However, PAR may not entirely limit the latitudinal expansion of the annual average probability of a few taxa. For instance, Alexandrium spp. and A. tamarense exhibit higher presence probabilities in the northern region, even though polar nights are more prolonged.

Blooms of Alexandrium spp. and Dinophysis spp. are commonly associated with stratified waters (Klemm et al., 2022; Reguera et al., 2012). The simulated response supports these associations, as the presence probability of D. acuminata, Alexandrium spp., and A. tamarense increases in shallower MLD—commonly correlated with more stratified waters. Note that MLD is also commonly associated with other important features related to HABs, such as nutricline depth and upwelling (Hällfors et al., 2011; Paulino et al., 2018; Peralta-Ferriz and Woodgate, 2015; Rial et al., 2023). The presence probability in shallower MLD can be up to four times higher than in deeper MLD, as observed for the response of A. tamarense. Besides, the probability of Pseudo-nitzschia spp., P. reticulatum, A. spinosum also shows an increase in shallower MLD. In contrast, the presence probability of D. acuta increases in shallower waters but does not peak at the lowest MLD values, and D. norvegica shows relatively low sensitivity to water stratification. Similar to PAR, the probability response to MLD may contribute to the seasonal detection patterns. MLD is typically deeper in winter and begins to shallow in spring due to surface heating and freshwater input (Peralta-Ferriz and Woodgate, 2015). Therefore, MLD is expected to contribute to the lower presence probabilities during winter for most species.

The evaluated taxa exhibit an increased presence probability in more saline waters within the salinity range from 25 to 34.5, except for D. norvegica. Therefore, the potential impact of salinity on the geographic distribution of toxic algae should be considered. D. norvegica shows the highest annual probabilities in the Skagerrak Strait, where relatively fresh waters are prevalent due to inflows from the Baltic Sea (Eldevik et al., 2009; Furevik et al., 2002). Note that higher D. norvegica presence toward relatively fresher waters, in the 25–34.5 PSU range, was also found in the Canadian coast (Boivin-Rioux et al., 2022). However, it should be acknowledged that salinity is a good tracer of water masses of different origins in the Norwegian Sea. The NCC water is fresher as it comes from river inflow and the Baltic Sea than the NwAC of tropical and Atlantic origin. This might have an influence on the occurrence of D. norvegica. Additionally, the highest probabilities of Alexandrium spp. and A. tamarense in northern Norway can be attributed to the lower freshwater input compared to southern Norway (Frigstad et al., 2020; Furevik et al., 2002).

The response to SST exhibits the greatest variability among the four inputs employed. It is important to note that the response to SST is of great concern, as it is commonly assumed that an increase in temperature may lead to increased prevalence and intensity of HABs (Wells et al., 2020). However, our modeled response from 3.7 to 18.8 °C suggests that temperature may primarily influence the selection of the most common taxon rather than causing a general increase of all toxic algae presence probability. Therefore, without considering changes in SSS, MLD, and PAR, an increase in temperature could favor taxa associated with warmer waters at the expense of those related to colder waters. Further investigation is needed to explore this aspect, which we intend to address in future studies. In terms of geographical distribution, SST may also explain why the highest probabilities of Alexandrium spp. and A. tamarense occur in colder regions of northern Norway, while D. acuta and D. acuminata are more prevalent in warmer southern Norway. Notably, D. norvegica shows an increase in probability with lower temperatures but is still more abundant in southern Norway, likely due to the strong influence of low salinity as previously mentioned. Finally, Pseudo-nitzschia spp. exhibit a wide range of SST tolerance and that might explain its increased presence probability along the entire Norwegian coast.

4.3. Improving local monitoring and mitigation actions

The current monitoring protocol assesses the occurrence and levels of harmful algae taxa on a weekly basis and measures their respective toxins in shellfish monthly. However, our analysis presented in Figures 8, 9, and 10 indicates that the presence and HA probabilities vary significantly across geographical regions and seasons for the different taxa. For instance, the annual presence probability of A. tamarense is five times higher in northern Norway compared to southern Norway, and the presence probability of D. acuta is up to six times higher in autumn than in spring. The current monitoring effort remains consistent throughout the main algae growth season, regardless of the likelihood of harmful algae risks. During periods of elevated presence or HA risk, it could be beneficial to increase the monitoring frequency, employ faster analysis methods for algae abundance and toxin concentration, and provide general guidance to the public regarding the consumption of wild shellfish. The aquaculture industry could better optimize their harvesting plans and thus contribute to increasing safety. Furthermore, climate decadal variability is large and predictable in the region (Smith et al., 2020). Refined information on the toxic algae evolution in the coming decade could be relevant for the farmers and the monitoring agency.

In our study, we calibrate probabilistic models for eight toxic taxa monitored along the Norwegian coast, ranging from the presence to HA. Models focusing on HA probability are more desirable as they are more closely related to toxic accumulation in shellfish, while also being able to pinpoint the most dangerous periods (Figure 10). HA probabilistic models of Alexandrium spp., A. tamarense, D. acuta, and A. spinosum can be applied as they are well correlated with their observation frequencies. However, the HA models for other taxa are not correlated with observation frequency and should not be employed. For D. norvegica, the model for abundances exceeding 3600 CellsL⁻¹ (80th percentile for HA) demonstrates satisfactory performance and can be utilized for tailored assessments. In the case of P. reticulatum, only the presence model exhibits good results. Nonetheless, it still provides valuable information for monitoring purposes. For instance, seasons with nearly 0% probability of presence (Figure 9g) suggest periods that are safer from P. reticulatum toxin contamination. The presence model for D. acuminata lacks the ability to estimate periods with a close-to-zero probability. However, it can still identify seasons with varying likelihoods of detection. As for Pseudo-nitzschia spp., the HA model needs improvement to provide useful information. Although the presence probability of Pseudo-nitzschia spp. demonstrates good correlation with observed frequency, this genus is commonly found in the Norwegian waters in all productive seasons (Hasle et al., 1996) and makes the presence probability model less relevant.

An even more preferable product would be the probability of toxin accumulation in the blue mussels, an approach not addressed in this study for two reasons. First, the modeling becomes more complex as the toxin accumulation depends on the algae abundance and toxicity, and on the mussel feeding and starving periods (Aasen et al., 2005; Smith et al., 2018; Röder et al., 2012; Lindahl et al., 2007; Nielsen et al., 2016; Svensson, 2003; Duinker et al., 2007). Second, toxin data are collected monthly, and thus, there is less data to train the machine learning models, roughly 25% that of algae abundance data. This amount of data is insufficient to produce skillful probabilistic models, as it must include training and testing datasets large enough to be evaluated on reliability diagrams. Nevertheless, we foresee the probability modeling of toxin accumulation in the blue mussels as a viable option in the future when more data becomes available.

5. Future perspectives

This study develops SVM probabilistic models for estimating the probability of eight toxic algae along the Norwegian coast. The models, ranging from the presence to HA, estimate well the probability of presence and particularly excel for HA levels of Alexandrium spp., A. tamarense, D. acuta, and A. spinosum. Feeding the probabilistic model with observations or predictions of SST, MLD, SSS, and PAR can provide crucial insights into periods and regions of increased risk, enabling local authorities and actors to devise enhanced monitoring strategies to prevent shellfish poisoning outbreaks and mitigate production losses in shellfish farms. Predictions can be provided by subseasonal-to-seasonal forecasts and seasonal-to-decadal predictions (Meehl et al., 2021) that are relatively skillful in the region (Doblas-Reyes et al., 2013; Passos et al., 2023; Langehaug et al., 2017; Wang et al., 2019; Bethke et al., 2021). These forecasts are probabilistic (provided as ensembles) and can be easily used by our probabilistic model. The model can also be used to infer long-term projections of toxic algae frequency occurrence by using the future projection of SST, MLD, and SSS (Carvalho et al., 2021; Davy and Outten, 2020). The integration of probabilistic models could empower decision-makers with evidence-based tools to proactively safeguard public health and ensure the resilience of coastal ecosystems, mitigating the impacts of toxic algae contamination and promoting sustainable and safe shellfish production and consumption. How decisions would be made based on the probabilistic output and how they should be communicated is yet to be planned based on continuous dialogues, which includes farmers and monitoring authorities.

Notably, our method can be extended to other coastal locations and help improve the safeguard of public health and reduce economic impacts. The machine learning modeling relies only on environmental drivers commonly available through remote sensing and modeled reanalysis, which are not limited to a specific region or national territories and can be obtained for various coastal waters. HAB observations used for training the model come from the Norwegian national program of toxic algae monitoring, which has similar infrastructure and data collection practices as in other regions. The calibrated probabilistic models for harmful algae presence perform well with eight taxa, suggesting a large possibility of expanding to other species monitored in other locations. Extending our proposed method to other coastal waters has the potential to enhance other monitoring programs and proactive mitigation actions to protect public health and the aquaculture industry.