The three demonstrated HAB modelling systems are from three different partners in PRIMROSE : Plymouth Marine Laboratory (PML), Institut Français de Recherche pour l’Exploitation de la Mer (IFREMER), and Instituto Español de Oceanografia (IEO), and will be referred to by the respective institute initials throughout. Whilst the domains of these models intersect, they are primarily focused on different areas: the Western Channel and Celtic Sea (PML), the Western Channel and the eastern Bay of Biscay (IFREMER) and NW and N Iberia (IEO) ( Figure 1 ). All the model systems are formed of an operational hydrodynamic model with a separate particle tracking model which is run offline utilising the output from the hydrodynamic model. None of the particle tracking models include explicit biological behaviours (e.g. diurnal migration, life cycle, nutrient limitation) and are purely passively advecting particles.

We have identified various bloom events of different HAB species to understand factors involved in the design of a HAB forecast system (see Table 1 ). When different model configurations with different spatial and temporal resolutions were available, i.e. blooms in the Western Channel in 2003 and 2020 and in Dinophysis spp. bloom in Iberia in autumn 2013, different particle tracking simulations were run forced by those models and compared. Two events in the Western Channel in 2010 and 2015 had remote sensing data available over a period of weeks, and this allowed the investigation of the period over which an advection-only forecast might be applicable, and how the location of the triggering area affects the results. Examples from 2013 and 2020 of Iberian Dinophysis spp. blooms demonstrate the integration with in-situ measurements of HAB species, biotoxins and closures at aquaculture monitoring sites. Finally, an event in 2020 in the Channel shows the factors to consider in a nearshore advection event and the importance of wind forcing, which is also discussed in the other events in the Western Channel.

To understand the spatial development of high biomass blooms in the Channel we have used remote sensing measurements of chlorophyll-a, in particular, daily composites retrieved from the Ocean Colour component of the European Space Agency (ESA) Climate Change Initiative (CCI) project (Sathyendranath et al. 2019). The remote sensing product has a spatial resolution of 4 km and combines observations from multiple sensors (MERIS, MODIS Aqua, SeaWiFS LAC & GAC, VIIRS, OLCI). These have been compared with the particles advected in the Lagrangian models to understand the change in spatial extent and location of the blooms. In-situ plankton observations from various sources have been included to allow species identification for the chosen blooms.

Several metrics have been used to quantify the model performance with respect to the satellite observations (in the case of high biomass blooms), and to compare the results of the different model configurations.

For the high biomass blooms that can be identified from satellite observations, the polygons representing the blooms have been automatically identified as the areas above a threshold of 5 mg/l, corresponding to the threshold used in Mishra and Mishra (2012) to define a ‘severe bloom’ on their index: all points on the grid of remote sensed data were determined whether they were above or below the threshold, and those above were connected to form polygons outlining the areas which were entirely above the threshold. The same has been performed on gridded data of particle counts obtained from the Lagrangian model for purposes of comparison with the satellite data, using a threshold number of particles scaled to the total number of particles run. In all cases, days with high cloud cover (>30%) were omitted, and for the days the calculations were performed, the model data was masked to omit data corresponding to areas covered by clouds in the observations. The centre of mass is defined as the mean longitude and latitude, weighted by chlorophyll concentration (particle count) of the observations (model). These centres of mass were compared to the centroids of the observed and modelled polygons and found to give very similar results, so centroids are not shown. Another metric is the percentage of model particles which lie within the observation polygons for a particular day. The dispersion is calculated for both the observed and modelled polygons by taking the change in total polygon area between two snapshots and dividing by the time (the time spacing between snapshots varies depending on which days have been omitted due to cloud cover). The final metric is the percentage overlap between the two polygons, defined as the percentage of the observed polygon which intersects with the modelled polygon.

The D. acuta and D. acuminata blooms that occurred in Portugal and Galicia in 2013 were not high biomass blooms that could be detected from remote sensing. In this case, we compared the results of the different model configurations by using as metric some of the Dispersal Kernel indices described in (Woillez et al. 2009). In particular, the Centre of Gravity (CG), which represents the mean location of the particles at a certain time; the Positive Area (PA), which is the sum of the area units that contain, at least, one particle; the Equivalent Area (EA), which is the area that the particles would occupy if they were distributed at homogeneous densities; the Coefficient of Variation (CV) defined as CV=PA/EA and the mean distance between the mean start and mean end particle positions. These indices, among others, have been used in previous studies to estimate dispersion patterns associated with Lagrangian models e.g. (Huret et al. 2010; García-García et al. 2016).

The PML setup was compared to two historical blooms which occurred in the Western Channel region in 2010 and 2015. The blooms were selected on the following criteria: blooms within the available model period, availability of recorded HAB species from in-situ sampling and with enough cloud free remote sensing data to be able to compare the model predictions to the development apparent in the observations. The model setup was identical to the operational setup, with underlying currents and diffusion coefficients from the FVCOM hydrodynamic model, the only change being the lateral boundary forcing came from the reanalysis version of prior lower resolution CMEMS model (AMM7 - O’Dea et al. 2012), since historical forecasts for AMM15 were not available.

A bloom of K. mikimotoi in the Western Channel, identified from FerryBox data (Smythe-Wright et al. 2014; Barnes et al. 2015) occurred in summer 2010. The bloom was evident on satellite chlorophyll data ( Figure 2 , column A) for a considerable period (July 5^th - August 4^th). Regions of high chlorophyll (threshold 5mg/l as above) in the remote sensing were identified to define the start locations for the model runs. Since the bloom was long lived, and to disentangle advection from evolution, two separate start times 2010-07-05 ( Figure 2 , column B) and 2010-07-15 ( Figure 2 , column C) were used, which correspond to points of clear remote sensing images. The second release was chosen for the 15^th of July because at this date the particles released during the initial seeding (the 5^th of July) had significantly diverged from the observed bloom.

The model shows limited ability in representing the evolution of the bloom at the early stages; particles move southwards in the model and this only matches one part of the observed bloom, since it extends considerably further to the east and north. This is reflected in the rapid drop in particles in polygon and polygon overlap metrics ( Figure 3 ), though the spatial split means the centres of mass diverge in a more linear manner. Initiating a second release 10 days later results in a better performance for the first 5 days (improved across all statistics, Figure 3 ) after which the model shows faster nearshore advection than the observations close to Brittany ( Figure 2 ). In this latter part of the second release the model is no better than the original release, with observed and model polygons having low overlaps and increasing centre of mass separations. Dispersion for both runs is similar to the observations, but with considerable variability ( Supplementary Material: Figure 1 ). Sensitivity to initial release time on shorter time scales was also studied by running experiments at hourly intervals across a six hour period but found not to significantly affect the results ( Supplementary Material: Figure 2 ).

The second event considered is a bloom in September 2015 of the small HAB dinoflagellate Prorocentrum cordatum (minimum) that can form extremely dense blooms. P. cordatum is detected quite commonly in some monitoring areas along the southern coast of England, but with variable cell densities from year to year (Turner et al. 2017), data available in (UK Food Standards Agency 2020) and it is also commonly found in the summer in the western English Channel at the L4 time series station (Widdicombe et al. 2010). Observations from the Western Channel Observatory (WCO) in September 2015 report concentrations of >5,000,000 cells per litre (Hiscock et al. 2016). The bloom is clearly visible on the remote sensing data as a high chlorophyll-a patch between the 10^th of September and the 3^rd of October ( Figure 2 , column D) and our particle tracking simulations of the evolution of a potential HAB were run. Again, two separate starting points based on a high chlorophyll threshold in the remote sensing were used to initiate the particle tracking model: 2015-09-10 ( Figure 2 , column E) and 2015-09-26 ( Figure 2 , column F). The performance of the model in the 2015 event was broadly similar to the 2010 event. The initial seeding failed to recreate the north and eastward evolution of the bloom, leading to swift drops in the overlap metrics and a separation in the centre of mass too ( Figure 3 ). The model was also reseeded later in the bloom (26^th September) and showed better agreement with the advection over the first few days (reduced centre of mass separation and less steep decline in particles in polygon percentage). However, by the fifth day the statistics were similar for both runs. The simulations show that particles reach the English southern coast, and this is in accordance with in situ observations: P. cordatum cells were measured in HAB monitoring samples taken along the southern coast of England in late summer 2015 (Turner et al. 2017), and high numbers (>1,000,000 cells per litre) were recorded at the end of September in St. Austell Bay, Cornwall (Hiscock et al. 2016; UK Food Standards Agency 2020).

The IFREMER model was used to simulate two events: one historical (2013) described in this section and one during the PRIMROSE project in 2020 (see section 3.6). A well documented bloom of K. mikimotoi occurred in the Channel in 2003, which resulted in mortality of many wild fish species. (Vanhoutte-Brunier et al. 2008) used an Eulerian passive tracer model to simulate this event with an older MARS3D configuration and concluded that apparent transport towards the coast must have been driven by progressive occurrence of favourable conditions for growth, and not by advection (since the tracer took a northward path). In this work, the new operational IFREMER model is used, which is an upgrade of the MARS3D system with several key improvements: increase in spatial resolution, thirty vertical layers (instead of twelve), a wider geographical domain and more rivers. Tracer transport is simulated with a Lagrangian instead of an Eulerian approach. The IFREMER model was initiated from and compared to the tracers released in (Vanhoutte-Brunier et al. 2008) ( Figure 4A ) In the IFREMER model ( Figure 4B ), from the initial release on 23 June, particles are transported westwards (see particle positions on the 7^th of July) and then eastwards (see position of the particles on the 21^st of July). The remote sensing data for the period shows the westward extension of the bloom matched by the particles ( Figure 4C ), but also bloom activity to the east. Unfortunately, cloud cover after the 13^th of July hampered satellite imagery, which therefore could not provide information on the subsequent transport of the bloom. Comparing the overlap with the bloom area and the distance between the centroids of the shapes, the new model demonstrates that the separation between the centroids was much less than the estimate for the (Vanhoutte-Brunier et al. 2008) model ( Figure 4 and Supplementary Material: Figure 3 ), since it captures the western movement, but the overlap of the polygons is not much greater since it does not capture the expansion to the east. The new run is also not able to reproduce the bloom reaching the French coast of Brittany by the beginning of August. The French Phytoplankton and Phycotoxin Monitoring Network (REPHY) counted 405 000 cells l^-1 in the Saint-Brieuc Bay on the 15^th of August 2003 (REPHY 2021). As suggested by (Vanhoutte-Brunier et al. 2008), our simulation also indicates that the Karenia cells found at Saint-Brieuc result from a coastal event related to local growth and not to the advection of the bloom observed offshore.

In the second event explored with the IFREMER model, a bloom of Mesodinium rubrum occurred in the Bay of Biscay during PRIMROSE. Mesodinium spp. is associated with blooms nearshore, in estuaries and embayments (Johnson et al. 2013; Trowbridge et al. 2017) and can persist for considerable periods (Crawford et al. 1997). A bloom was observed in spring 2020 at two local harbours (Oleron and Sables, Figures 5A, B respectively) and created considerable public concern due to water discolouration during the period of confinement due to COVID. The presence of M. rubrum was confirmed by the French HAB monitoring system REPHY, and this bloom was followed by the detection of Dinophysis spp. The ad-hoc operational model was run during this event to investigate the potential connection between the two sites. Particles were released at two locations on two different start dates: at Sables harbour where discolored waters were observed first to assess if there was an advective connection to the later observation at Oleron (05/04/2020 12:00). The second drift was initiated from Oleron, to estimate displacement of the bloom. Locations of particle release at Oleron were estimated from a satellite index of bloom (Normalized Difference Chlorophyll Index, NDCI as defined by (Mishra and Mishra 2012), values over 0.5 are considered as a blooming area) computed from a satellite product (Copernicus Sentinel 1A, 18/04/2020 10:56) ( Figure 5D ). Hydrodynamics were provided by a MARS3D configuration (horizontal resolution of 2.5 km) forced by ARPEGE-HR solutions (hourly, 0.1°, Meteo France).

Particles released at Sables were not advected to Oleron area ( Figure 5C ), which suggests the bloom was spread over a larger area. In the same figure panel, we can see that particles released from Oleron strongly drifted offshore after 8 days, possibly driven by a short period of strong offshore winds after the 18^th of March ( Figures 5E, F ).

The autumn 2013 Dinophysis acuta event on the northern Portuguese coast and Galician rías was the most intense since 2005 (Díaz et al. 2016), and D. acuta appeared together with D. acuminata. The HAB caused the closures of most of the harvesting polygons in the area for several months including the Christmas period, one of the most profitable of the year. The along-shore transport of D. acuta from the Portuguese shelf to the Galician rías has been described as the cause of sudden autumn blooms (Escalera et al., 2010; Pitcher et al., 2010; Díaz et al., 2016). D. acuta blooms usually occur in the summer in the shelf off Aveiro and get transported northwards. Several events of D. acuta autumn transport to the Galician Rías have been described, particularly those in 2005 and 2013 (Díaz et al., 2016). Table 2 summarizes the timeline of the sequence of closures from Portugal to the Galician bivalve harvesting polygons in 2013.

Lagrangian modelling studies of the 2013 event have been previously carried out (Ruiz-Villarreal et al. 2016; Silva et al. 2016; Moita et al. 2016) to try to elucidate if advection could explain the transport of D. acuta and the timing of closures from Portugal to the Galician rías. All model exercises predicted northwards transport, although no detailed analysis of the impact of the model configuration in the predicted particle transport was performed. In this paper we revisit this HAB event to evaluate the impact of using hydrodynamic models of different resolutions and configurations to force the Lagrangian model. We used four different model configurations: the IEO RAIA in the parent 4 km grid (IEO RAIA) and the child 1 km grid (IEO RAIA nested); the IEO BIO 3.5 km grid and the MeteoGalicia 2 km operational forecast configurations. Model set ups differ in the extent of the model domain (see Figure 1 ) and also in forcing at the open boundary (RAIA is forced with a climatology and IEO BIO and MeteoGalicia with Mercator), and river input.

Figure 6B shows the release of particles at the Portuguese polygon close to Aveiro the 17 ^th of September 2013 (blue particles), the date for which the maximum concentration of D. acuta was detected at Aveiro ( Table 2 ). One week later (the 24 ^th of September 2013), the particles were mostly located to the south and offshore the initial release position (red dots), due to the dominant upwelling winds ( Figure 6A ). This dispersion to the south is higher in the coarser resolution models ( Figure 6B , IEO BIO and IEO RAIA). Note that in the case of the finest model resolution ( Figure 6B , IEO RAIA nested) the particles did not displace further south because of the size of the model domain (see Figure 1 ).

After the 24^th of September 2013, strong downwelling conditions prevailed in the area ( Figure 6A ), which explains the northward transport of particles that reached the south of Ría de Vigo on the 30^th of September 2013 (green particles in Figure 6B ) in the case of the IEO RAIA configuration, and further north in the case of the IEO BIO, MeteoGalicia and the high-resolution configuration (IEO RAIA nested). The average transport northward is higher in the case of the mid-resolution model MeteoGalicia (as shown by the trajectory of the CG (black line) in 6B) because the transport to the south under upwelling conditions the previous week was less intense than in the coarse resolution models IEO RAIA and IEO BIO. The same happened in the case of the high-resolution model IEO RAIA nested, but in this case it was affected by the limits of the model domain.

Finally, on the 4^th of October, when the polygons started to close in the north of Portugal and the southern Galician Rias, we can see that all the models show particles (yellow in Figure 6 ) located at the mouth of every ria and even further northwards in the case of the IEO BIO, MeteoGalicia and IEO RAIA nested models. Hence, the results shown here are compatible with the sequence of events compiled in Table 2 .

Figure 6C illustrates the position of the particles in the first and second weeks after their release close to Aveiro the 24^th of September (blue particles). In this case, downwelling conditions were dominant ( Figure 6A ), and just one week later (30^th of September, red particles), the particles are close to the mouth of the ria de Vigo for the IEO RAIA model ( Figure 6C , second column), up to ria de Arousa for the IEO BIO and the MeteoGalicia models ( Figure 6C , first and third columns) and up to ria de Muros in the case of the IEO RAIA nested domain ( Figure 6C , fourth column). The 4^th of October (two weeks later, green particles), the particles are occupying the mouths of all the rías, except Muros, when using the IEO RAIA model as forcing. With the IEO BIO configuration, particles move further north (up to Malpica) and, in the case of the MeteoGalicia and IEO RAIA nested models, the particles reach up to the Artabro Gulf. These results are also compatible with the sequence of events described in Table 2 and do not differ much from those for the 30^th of September and 4^th of October in Figure 6B , although in the latter case, the density of particles reaching the Galician Rías would be higher, since the transport was always in the same direction and confined onshore. This is clear in the trajectory of the CG (black line), which shows a northward transport for all the configurations, with the MeteoGalicia model showing the longest transport and the IEO RAIA the shortest. Figure 6 illustrates that the mechanism of northward transport from Portugal persisted over several weeks.

The effect of the different hydrodynamic model configurations on the particle transport is summarized by means of some Dispersal Kernel indices in Table 3 for the 17^th and the 24^th of September 2013 releases. These indices give an idea of the particle dispersion. For instance, if we consider the simulations starting the 17^th of September, we see that the IEO RAIA configuration is the one producing the highest dispersion of particles (highest PA=13,750 km²), followed by the IEO BIO (PA=13,575 km²) and the IEO RAIA nested (PA= 6,850 km²) and MeteoGalicia (PA= 6,800 km²) models. This is clearly reflected in Figure 6B , where we can see that the IEO BIO and IEO RAIA configurations produce a southward transport of particles first, followed by a northward transport, which results in a larger area occupied by the particles. The MeteoGalicia configuration, which is run in a model domain as large as IEO RAIA ( Figure 1 ), does not produce the initial movement of particles southwards, which results in a smaller occupied area and an overall transport further away from the mean initial location in the northward direction (compare the trajectories of the CG, black lines in Figure 6 ). The particles are also more homogeneously distributed when using the IEO RAIA configuration (highest EA). The CV indicates the existence of more homogeneous densities (lower values) or aggregates (higher values), meaning that the IEO BIO configuration produced more aggregates (CV=6.1), followed by the IEO RAIA (CV=4.51) and the MeteoGalicia (CV=3.94) configurations. In terms of mean distance between the mean particle origin and the mean particle end, it is the MeteoGalicia model the one that shows higher values (159.8 km), followed by the IEO BIO (65 km) and the IEO RAIA (42.5 km). The indices for the IEO RAIA nested model are influenced by the fact that the particles get accumulated in the southern limit of its domain (see Figure 1 ) during the first days of the simulation, which results in a PA similar to the MeteoGalicia model, but with much higher aggregation (CV=10.34).

For the simulations starting the 24^th of September 2013, the IEO BIO shows the highest dispersion (PA=12,650 km²) due to a number of particles being transported southwards. In the remaining configurations, this southward transport is minimal (IEO RAIA, IEO RAIA nested) or non-existent (MeteoGalicia) and the dispersion area increases with resolution. The degree of aggregation is similar in all the configurations (CV around 3), although slightly higher for IEO BIO. The longest mean transport (229.4 km) corresponded to the MeteoGalicia configuration (which did not show southward transport at all), followed by the highest resolution model (IEO RAIA nested, 160.33 km), for which the mean distance was similar to the IEO BIO model (157.63 km).

During PRIMROSE, a transboundary event involving a proliferation of Dinophysis acuta and Dinophysis acuminata that lead to harvesting polygon closures in North Portugal and the Galician Rías in September 2020 was investigated using the IEO early warning system, forced by the ROMS MeteoGalicia forecast model. Table 2 summarizes the time sequence of toxic phytoplankton and biotoxin measurements during this event, together with the dates when closures were enforced by the regulatory authorities. Figures 7A, B show the evolution of the concentration of D. acuminata and D. acuta, respectively, during the month of September 2020, and Figure 7C depicts the concentration of DSP biotoxin. The DSP toxin was already above regulatory threshold permissible level (160 μg/kg) in the first week of September in Portugal (caused by D. acuminata), as well as in some polygons in Ría de Pontevedra, which were already affected for several weeks due to the presence of both species of Dinophysis ( Table 2 ). The second week of September (second column Figure 7 ), some polygons at the mouth of Ría de Vigo exceeded the DSP threshold, but it is during the third week (third column) when high concentrations were measured inside the ría causing most of the polygon to be closed due to the accumulation of both species. In the following week (fourth column), D. acuminata and D. acuta blooms closed polygons in the Rías de Arousa and Muros.

The IEO early warning system was used to investigate whether: a) there was a potential transport from the Portuguese shelf to the Galician Rías and this caused the closures of harvesting areas and b) the sequence of closures from south to north in the Galician Rías was caused by exchange between adjacent rías.

The upwelling index evolution from the 4^th of September, when the first polygons were closed in North Portugal, shows that winds were upwelling favourable and relaxed around the 8^th of September ( Figure 8B ) during neap tides ( Figure 8A ). The Lagrangian model forecasts for the 4^th of September do not suggest northward transport on the Portuguese shelf ( Figure 7D ), and neither do the forecasts for the subsequent days before the closures in Ria de Vigo after the 16^th of September (not shown), probably due to the rather calm wind conditions. Polygon closures extended from Vigo to Arousa from the 16^th of September onward, just a few days after the occurrence of neap tides and still under rather calm wind conditions ( Figures 8A, B ). The results of the Lagrangian forecasts for the 16^th of September at Ría de Vigo ( Figure 7E ) show that particles mostly remained confined inside the ría, and that the exchange with Ría de Pontevedra and the rías further northwards was unlikely, based on the model results.

Just a few weeks after the polygon closures, conditions changed and most of the closed polygons in Rías de Vigo, Arousa and Muros reopened between the 9^th and the 17^th of October. Figures 8A, B show that this was a period of neap tides under strong upwelling conditions. Upwelling conditions are favourable for the offshore transport, and hence the outflow and exportation of toxic cells, which would produce a “cleaning effect”. Figures 8C, D show that D. acuminata and D. acuta concentrations reduced during that period, except for Ría de Pontevedra, where D. acuminata persisted. The reduction in the DSP biotoxin throughout October is clear in Figure 8E for all the rías, except Pontevedra. The IEO early warning system predicted conditions of transport outside the rías (see Figure 8F for Ría de Vigo), and concluded that polygon re-openings were likely. Forecasts of trends to open/close of harvesting areas are another service provided by the forecasting system, already demonstrated in ASIMUTH project (Maguire et al., 2016; Ruiz-Villarreal et al., 2016; Silva et al., 2016).

A potential large bloom was spotted in the Channel region in June 2020 during PRIMROSE and since it incurred on both PML and IFREMER model domains, model runs using both the PML and IFREMER (2.5 km) model frameworks were undertaken. The full operational setups as described in section 2.1 were used. The identification of the bloom came from the algorithmic Karenia spp. risk identification from satellite images based on algorithm by (Kurekin et al. 2014) introduced in the section 2.1.1. The HAB risk map in the PML HAB Risk portal¹ shows the composite image over several days. However, for sensible comparison to the higher time frequency involved with advection, the respective model runs were compared to the algorithm output from individual passes of the satellite products ( Figure 9 ). Unfortunately, whilst the initial satellite images were clear and the algorithm was able to identify a HAB, the following days were affected by cloud coverage. In any case, the eastern part of the bloom remained in the visible part of the images and thus, the satellite could be used to assess the skill of the two models.

Both models were started from points the operational HAB risk algorithm deemed as high risk. The PML model shows a much greater spread of particles from the western starting positions, where they reach the higher velocity currents around the Channel Islands, and also a closer incursion towards the coast from the eastern starting points. The significant difference driving this is that the offline forcing of the IFREMER Lagrangian model in this case is saved at a daily resolution, so tidal frequency is not resolved (though will include underlying tidal features in the hydrodynamical model).

Several case studies were focused on high biomass blooms such as those of K. mikimotoi in the Western Channel with the objective of evaluating how a short term advective only forecast performs. K. mikimotoi is frequently present in the summer phytoplankton community in the Channel (Widdicombe et al., 2010) and it is thought to be responsible for occasional mass finfish and benthic mortalities (Barnes et al., 2015). In sections 3.1, 3.2 and 3.6, Lagrangian particle trajectories were compared to satellite imagery following the metrics described in section 2.3.1.

Remote sensing products have several challenges for direct use in model comparison. Firstly, remote sensing is not a perfect measure of a bloom, especially for toxic species. K. mikimotoi for example is known to produce a thin layer which may not show up on satellite imagery (McManus et al. 2008; Brand et al. 2012), and nearshore chlorophyll measurements will be affected by sediments (Le et al. 2013). Also, Karenia species cell colouring imply that they have reduced reflectance, which makes it hard to identify species and will likely under count concentrations, especially once cell concentrations are above 10⁴ cells 1^-1 (Hu et al. 2005; Cannizzaro et al. 2008). Whilst studies have used chlorophyll products to study K. brevis blooms e.g. (Redalje et al. 2008), the Normalized Fluorescence Line Height (NFLH) product was shown to give better results (Hu et al. 2005; Cannizzaro et al. 2008; Tomlinson et al. 2009). In a review of various different approaches, (Soto et al. 2015) found further improvements when combined with other products, determining the Rrs - NFLH method using NFLH and the 555nm reflectance to perform best for identifying Karenia blooms. However, for these events coverage with chlorophyll-a products was better than the NFLH and 555nm product, and when comparing to available scenes from NFLH and Rrs - NFLH ( Supplementary Material: Figures 4, 5 ), the chlorophyll product shows a very similar bloom evolution. In this study, since we have identified the species from available in-situ observations and we are concerned with the bloom shape and movement rather than absolute concentrations, the use of a chlorophyll product is less problematic.

Another issue related to remote sensing is that data availability is variable due to cloud cover. This can be mitigated with the use of composites of several days together or interpolation between days of data availability. Here we are concerned with short time scale advection and hence, products composited over a longer period of time will average out this movement. Interpolated products should not have this problem, but must use some extra technique to fill in the gaps, for example using a hydrodynamic model [(Lin et al. 2021), this issue] or machine learning (Vandal and Nemani 2019) to produce intermediate images, which may itself introduce spurious movement. Here we decided to be conservative and only use products which are either daily composites or individual scenes, omitting days with high cloud cover (>30%), to ensure we are comparing to actual bloom development.

Despite these concerns, remote sensing products provide a spatial and temporal coverage of potential HAB events that cannot be achieved by alternative methods. For example, in-situ sampling requires structured field campaigns e.g. (Velo-Suárez et al. 2010) and has to be fortunate enough to coincide with a bloom [e.g. Jordan et al. 2021, this issue]. Opportunistic measurements, such as FerryBoxes provide only limited spatial resolution (Qurban 2009; Hartman et al. 2014), which is not sufficient to understand the advective component of mid-channel blooms. In future it is possible these limitations will be overcome by responsive sampling by unmanned platforms [e.g. as demonstrated by Ruiz-Villarreal et al. 2022, this issue].

If we assume that the remote sensing data provided an accurate representation of bloom development, then in both examples from the PML model (sections 3.1 and 3.6) the latter part of the bloom evolution was better forecast than the initial movement. In both cases the efficacy of the forecast dropped rapidly in the first five days after initialisation, as demonstrated by the various metrics applied for comparing particle model results and satellite data. The difference between the initial and latter movement could represent a period of rapid growth followed by advection, which would fit with rapid onset in the order of a few days followed by a maintenance phase that has been observed in some HAB species (Redalje et al. 2008; Aoki et al. 2012), and with the apparent change in surface area of the bloom during this period, though the latter could also be a consequence of diffusive processes. Onset here refers to onset visible on the remote sensing data; there may be prior activity at the subsurface or with a low chlorophyll content and hence not detected by the satellite, but this would be beyond the scope of a remote sensing triggered forecast system. If this scenario of growth during the early stage of a bloom followed by advection was the case, it would suggest that a forecast system based only on the advection of passive particles such as the one described here would not account for all the processes involved, and that adding some biological behaviour to the particles would improve the representation of the bloom evolution. In this respect, the individual-based modelling of the development and transport of the strong 2006 K. mikimotoi bloom in Scotland reported in (Gillibrand et al. 2016) showed the potential benefit from better parameterisation of temperature dependence of both growth and mortality, albeit over a longer time frame. Their simulations clearly also showed the importance of advection in HAB transport and they conclude that their results would be improved by improvement of spatial and temporal resolution of the underlying hydrodynamical model, which agrees with our conclusions. Additionally, models aimed at predicting the onset of blooms using sea surface temperature and wind indices (Cusack et al. 2015; Karki et al. 2018) could be included to indicate the reliability of a short term advective forecast at a particular time.

Another possibility is that the advection was poorly forecast in the initial stages compared to the latter, due to some change in the physical conditions. Wind conditions can play an important role both for onset - for mid channel blooms low wind speeds are a factor in K. mikimotoi bloom formation (Gentien et al. 2007) - but also for its impact on advection itself. Direct and indirect wind forcing (the latter via wind waves) are key to predict surface advection (Röhrs et al. 2021). The PML model does not account for the Stokes drift, the resultant Lagrangian transport from surface waves, and this can be a significant component of drift. The highest frequency components drive the largest drift effect at the surface of the water, and are often wind waves (Röhrs et al. 2014; Tamtare et al. 2021). The wind speeds during both years were modest, though in the latter portion of both years after the second release, wind speeds were lower than during the initial period ( Figure 10 ), which could suggest that either the advection is less well predicted during higher wind or the bloom dynamics are different.

We have presented several examples in which differences in model resolution and in other details of model configuration strongly influence particle dispersion and consequently HAB transport forecasts. The analysis of the autumn 2013 D. acuta event (section 3.4) clearly illustrates that different model configurations result in different trajectories of Lagrangian particles. We were able to compare four different model set-ups that differ in spatial resolution, model domain and model forcing (river input and open boundary forcing). Although there are some published results of particle tracking simulations of that 2013 event (Moita et al. 2016; Ruiz-Villarreal et al. 2016; Silva et al. 2016), the hydrodynamical conditions driving particle transport are only described in detail in (Ruiz-Villarreal et al. 2016), who report a northwards current on the shelf coinciding with the relaxation of upwelling winds around neap tides, which can be responsible of the along-shore transport of D. acuta from Portugal to Galicia. Surface currents and surface salinity at relevant dates during this event in our four different model configurations are plotted in Figure 11 , where it is noticeable that shelf circulation off the rias and representation of rivers differs. Shelf currents are directed northwards in all configurations although there are differences in speed and current location on the shelf, note especially the strong northwards velocities in the MeteoGalicia results. These differences in model configurations result in the differences in particle transport evident in Figure 6 and in Dispersal kernel indices, which we have seen in section 3.4. It is interesting to note that, in general, particles get advected more often to the north coast (longer mean transport) in the higher resolution configurations (MeteoGalicia and IEO RAIA nested), and this is especially clear when comparing IEO RAIA (4 km resolution) and IEO RAIA nested (1 km resolution). It is difficult to evaluate against observations the extent of the northwards transport in this event. The only available information is that D. acuta was not observed in the north coast HAB monitoring stations in the Ria de Ares, which could suggest that northwards advection was not so strong. The dependence on the resolution of the forcing hydrodynamic model has already been recognised in the literature, e.g. (Kvile et al. 2018; Nooteboom et al. 2020), who also report higher resolution appears to lead to greater transport, possibly from smaller scale effects being represented rather than diffused out. The impact of changing model resolution is also clear in the IFREMER 2003 model run (section 3.2), although the use of Eulerian and Lagrangian transport models in the two simulations does not allow us to isolate the impact of resolution in those results.

Differences in shelf circulation and dispersion of particles off the rias among the different IEO/MeteoGalicia model configurations could be partly attributable to the differences in river input, which are evident in surface salinity ( Figure 11 ). IEO RAIA and IEO BIO configurations consider river run off from all rivers in Galicia and in the Portuguese shelf. MeteoGalicia configuration is forced by model predictions from SWAT for Galician rivers including Miño, but on the northern Portuguese shelf it only includes a climatological monthly average run-off from Douro river. Moreover, in these hindcast simulations, IEO RAIA and IEO BIO configurations are forced by real daily river run-off data, while MeteoGalicia includes river run-offs from the SWAT forecasts, which tends to underestimate peak flows (Venâncio et al., 2019). In addition, IEO RAIA and IEO BIO configurations are free runs, resulting in the freshwater budget on the shelf being well represented, while in MeteoGalicia configuration, the re-initialization with MyOcean model (Costa et al., 2012) implies that the freshwater budget at scales longer than the forecast cycle comes only from the rivers considered in the MyOcean model (Minho and Douro). The lower influence of river plumes in the dynamics could explain the fact that the higher resolution MeteoGalicia (2km) configuration, run in a model domain similar to IEO RAIA and with open boundaries from an operational model (as IEO BIO), seems to favour the underestimation of southwards transport of particles and the overestimation of northwards currents as seen in Figure 11 , with the consequence that particles are mainly transported northwards in that model configuration.

The comparison of the IFREMER to PML model in the summer 2020 event (section 3.6) suggests that temporal resolution of the underlying hydrodynamic model used to force the Lagrangian model also has an effect. Whilst other studies have found temporal resolution to be potentially less important than spatial resolution (Qin et al. 2014; Kvile et al. 2018), here it meant that tidal currents were not included in the IFREMER Lagrangian simulations, and hence the physical forcing of the particle tracking was not represented equally. In spite of this, the limited evidence from the remote sensing suggests that the IFREMER model has done a better job in this case, at least on the eastern portion of the bloom, which was not advected onshore. This could indicate that tides were not the most relevant driver of particle transport in this particular case. There is also a difference in approach to diffusion, with the IFREMER model using an advective only setup in the horizontal and the PML model having spatially varying diffusion based on the underlying hydrodynamic model. However, for the short term forecasts, which are the aim of these operational systems, it is not expected that diffusion will have a great impact.

Our results showing divergence between models highlight the importance of the spatial and temporal variability of the model that forces the particle transport model. All these forcing models have lateral boundary conditions provided by coarser resolution models, the setup of which has an impact on the forecasts. In the September 2013 event in NW Iberia, we used different model set-ups and we could show that differences in model resolution, domain and river input impacted the transport model. Forcing at the open boundary was also different in models (climatology, one-way nesting or a larger operational model) and the comparison of IEO BIO and IEO RAIA results, showing that IEO BIO demonstrating higher variability offshore, indicates that forcing with an operational larger scale model impacts the resultant transport. A detailed sensitivity study of the impact of changes in lateral forcing in forecast models and in the associated particle transport is beyond the scope of this paper, since setting up operational boundary conditions for different upstream models is a complex process. Our results confirm however that it is necessary to evaluate the differences that model-set up (including resolution, freshwater input, and open boundary conditions) may induce in the dynamics predicted by the different hydrodynamic models, since the description of currents, frontal structures and well mixed/stratified areas may differ between circulation model setups and consequently have an impact on particle transport.

Early warning of the presence, location and subsequent evolution of HABs is the objective of a HAB alert system. A Lagrangian model based HAB alert system comprises three components: an initial detection component, for example from in-situ measurements or remote sensing; a Lagrangian model component to turn this into a forecast, the efficacy of which has been the focus of this study; and a dissemination component, to inform aquaculture producers and managers in charge of protection of human health of the risk that a HAB may affect aquaculture areas.

For HAB alert systems of high biomass species such as Karenia spp., which usually appear offshore in areas not routinely sampled by existing in-situ monitoring (concentrated in aquaculture production areas on the coast), detection of HAB risk to initiate the Lagrangian model comes from satellite imagery. The approach of Karenia spp. risk detection from remote sensing and then execution of particle tracking simulations is in use in the eastern Gulf of Mexico (Stumpf et al., 2009), and was demonstrated in the Scottish coast during ASIMUTH (Gillibrand et al., 2016). In this manuscript, we have illustrated in Section 3.6 that this approach has utility for assessing the advection of potential HAB species in the English Channel, especially K. mikimotoi, which appears frequently in the summer (Widdicombe et al., 2010; Barnes et al., 2015). This HAB forecast method is affected by the accuracy of the satellite products to detect a bloom and discriminate the particular species, which can be harmful or not. In the example in Section 3.6 although the alert was activated by the Karenia spp. risk and models were run to provide forecast of an eventual HAB, no bloom of Karenia spp. at the in situ monitoring sites were observed see (Atkinson et al. 2021). HAB detection algorithms are being constantly improved to distinguish harmful from non-harmful high chlorophyll blooms (Cannizzaro et al. 2008; Tomlinson et al. 2009; Sourisseau et al. 2016), but there are still considerable challenges in species discrimination (Kurekin et al. 2014; Feng et al. 2022) and in distinguishing high turbidity from high chlorophyll concentrations of HAB species (Martinez-Vicente et al. 2020), this issue, and this will remain a source of error independent of the errors from the Lagrangian model.

We have also shown in this study (Sections 3.1) that remote sensing combined with Lagrangian particle modelling is useful to follow blooms of other potential HAB species like Prorocentrum cordatum (minimum), which can form highly dense monospecific blooms. Blooms of P. cordatum (minimum) have been associated with anoxic/hypoxic events causing fish kills (Heil et al., 2005). Although there is debate about its toxicity (Heil et al., 2005; Turner et al., 2017), P. cordatum is currently monitored in bivalve aquaculture areas in Europe. Interestingly the HAB risk map aimed at Karenia spp. also indicated high likelihood of a bloom during this event², this could indicate a misidentification of species in this case and work discriminating the two, as has been done for K. mikimotoi and another Prorocentrum species (P. donghaiense) in (Feng et al. 2022), might be valuable for future HAB prediction. Finally, we have also demonstrated in section 3.3 that the Lagrangian particle simulations are a tool for tracking the evolution of high biomass blooms of the non-toxic Mesodinium rubrum, prey of Dinophysis spp.

For low biomass blooms like those of Dinophysis spp., initiation in the alert system relies only on the observations of the monitoring systems. If the monitoring systems detect the presence of the toxic species, it is an indication of the risk that a HAB could affect aquaculture sites. However, the monitoring only takes measurements of toxic species and toxins in aquaculture sites (located near the coast), but the bloom might be developing offshore, and then be transported to the sites. Therefore, the strategy chosen for HAB alerts of these species in the Galician early warning system is the constant release of particles in selected areas in Galicia as well as in Portugal as described in section 2.1.3. Our evaluation of the Galician system in September 2020 (section 3.5) confirmed that the Lagrangian-hydrodynamical coupled simulations provide predictions of favorable conditions of along-shore advection, exchanges between rías or flows in and out of the rías, and this is useful for characterising Dinophysis spp. HAB transport, confirming (Ruiz-Villarreal et al. 2016) conclusions. However, the results of the coupled hydrodynamical-Lagrangian simulations did not explain fully the sequence of events, which might be attributed to model uncertainties (misrepresentation of shelf dynamics, i.e river plumes, resolution not being high enough to resolve rías and fluxes between them, etc.), but also to the possibility of the HAB event being caused by local growth of D. acuminata and D. acuta and not by advective alongshore transport of D. acuta like in other autumns. The analysis of this HAB event clearly shows that although models are an important element of the early warning system, a HAB alert can only be issued in light of the measurements of the HAB monitoring system. When transnational alongshore transport is relevant, as in the presented simulations of Dinophysis spp. in Galicia and N Portugal (sections 3.4 and 3.5), transboundary exchange of information on HAB species and toxins between different HAB monitoring systems is crucial.

The final part of a HAB alert system is the dissemination to end-users. Hydrodynamic model simulations in the different areas are currently distributed via THREDDS (Thematic Realtime Environmental Distributed Data Services), a web data server that provides metadata and data access for scientific datasets using different remote data access protocols including Open Geospatial Consortium (OGC) standard protocols. The availability of model output in the THREDDS server allows us to obtain data via OPENDAP, HTTPserver or NetcdfSubset, but also via OGC Web Map Service (WMS) or Web Coverage Service (WCS) Interface Standards. However, a THREDDS server is not the best option for distributing non-gridded data, such as the Lagrangian model outputs. At IEO, the Lagrangian particle output files are transformed into shapefiles in order to be served through a GeoServer, which serves data using standard OGC protocols. In the GeoServer, the latest Lagrangian particle output from the daily run of the Lagrangian particle tracking simulations is served for visualization in viewers and data portals, including PRIMROSE data viewer. A similar approach is followed at IFREMER for distributing Lagrangian particle trajectories.

Presentation is a part of dissemination. Two differing presentations of Lagrangian trajectories for use in HAB forecast systems are shown in Figure 12 . The IEO model runs regularly using the same starting locations, which have been identified as possible sources (see section 2.1.3). These are presented in a portal with the individual tracks shown, which allows rapid visualization of trajectories ( Figure 12A ). The PML model is initiated from remote sensing HAB algorithms, which as above come with an associated uncertainty (Kurekin et al. 2014). Also, all the advection models include a diffusion term which adds an element of stochasticity to the results which represents sub-grid scale processes. The PML output combines both sources of uncertainty (from the detection algorithm and from the diffusion of the particles). In this way, each particle is associated with the probability from the detection algorithm, then kernel density estimates are fitted to each particle release to give the probability distribution from advection component. These two probabilities are combined using a simple Bayesian approach to produce a single map of HAB probability ( Figure 12B ). Whilst this attempts to include the available information on uncertainty from all inputs to the model, the result is less clearly a drift product, which makes it less easily interpretable, possibly resulting in less uptake by, and thus value, to end users. However, given the uncertainties of the drift, either presentation is improved by an accompanying interpretation. These forecasts are not usually presented in isolation; they more often form part of HAB bulletins or HAB report web sites, such as those produced in ASIMUTH (Cusack et al. 2016; Maguire et al. 2016; Ruiz-Villarreal et al. 2016; Silva et al. 2016) and PRIMROSE projects (Davidson et al. 2021; Fernandes-Salvador et al. 2021), this issue. HAB early warning system presentations add in local measurements, general oceanographic conditions and expert advice, along with advective forecasts, which help to mitigate the problems of interpretation and user uptake identified above.