Early Warning of Harmful Algal Bloom Risk Using Satellite Ocean Color and Lagrangian Particle Trajectories

Combining Lagrangian trajectories and satellite observations provides a novel basis for monitoring changes in water properties with high temporal and spatial resolution. In this study, a prediction scheme was developed for synthesizing satellite observations and Lagrangian model data for better interpretation of harmful algal bloom (HAB) risk. The algorithm can not only predict variations in chlorophyll-a concentration but also changes in spectral properties of the water, which are important for discrimination of different algal species from satellite ocean color. The prediction scheme was applied to regions along the coast of England to verify its applicability. It was shown that the Lagrangian methodology can significantly improve the coverage of satellite products, and the unique animations are effective for interpretation of the development of HABs. A comparison between chlorophyll-a predictions and satellite observations further demonstrated the effectiveness of this approach: r2 = 0.81 and a low mean absolute percentage error of 36.9%. Although uncertainties from modeling and the methodology affect the accuracy of predictions, this approach offers a powerful tool for monitoring the marine ecosystem and for supporting the aquaculture industry with improved early warning of potential HABs.


INTRODUCTION
Harmful algal blooms (HABs) occur in many coastal regions around the world and appear to be increasing in severity and extent (Hallegraeff, 1993(Hallegraeff, , 2003Grattan et al., 2016;Gobler, 2020). HABs have caused severe economic losses to aquaculture, fisheries, and tourism while creating major environmental and human health impacts (Anderson et al., 2000;Landsberg, 2002;Heisler et al., 2008). Toxic bloom-forming algae can cause wildlife mortality or human seafood poisoning, and even HAB species that do not produce toxins are able to cause harm through development of high biomass, leading to foams or scums, depletion of oxygen as blooms decay, or destruction of habitat for fish or shellfish by shading of submerged vegetation (Sellner et al., 2003). Such impacts from HABs pose a serious threat to aquatic ecosystems and can disrupt their associated food web (Fogg, 1969;Paerl, 1988). Therefore, considerable attention has been focused on methods to reduce the risks of HAB impacts (Sengco and Anderson, 2004;Anderson, 2009;Anderson et al., 2012).
Satellite ocean color sensors offer a means of detecting and monitoring HABs in the ocean and coastal zone. The potential value of remote sensing for HABs was first described by Mueller (1981), after an experimental ocean color sensor attached to an aircraft detected a bloom of Karenia brevis. As the instrument was developed to simulate the Coastal Zone Color Scanner (CZCS), launched in late 1978, this indicated the capability for satellite detection of blooms. Various approaches were further developed for detection and monitoring of HABs from satellite remote sensing (reviewed by Klemas, 2012;Blondeau-Patissier et al., 2014).
There are many advantages of satellite ocean color products compared to in situ monitoring: wide temporal and spatial coverage, and inexpensive; however, there are significant limitations. First, the coarse spectral resolution of current ocean color satellites is only sufficient to distinguish certain clear anomalies in bloom coloring, as the visible spectrum is mostly determined by optical properties of varying concentration of chlorophyll pigments. Second, data from satellite ocean color only reflect the state of the waters at the moment the measurements are taken, whereas the blooms are subject to physical forcing from tidal and wind-driven currents which will advect them away from this state. It is thus very challenging from satellite data alone to provide information regarding the future development and movement of HABs. Another limitation of ocean color remote sensing is that satellite observations are often hampered by weather conditions, such as clouds and sea fog, which can substantially reduce the number of valid ocean color pixels. As a result, satellite imagery may only partially capture features of HABs.
In recent years, the use of Lagrangian particle tracing (or other numerical models) for HAB monitoring has gained more interest (Olascoaga et al., 2008;Wynne et al., 2011;Son et al., 2015;Kwon et al., 2019;Li et al., 2020;Fernandes-Salvador et al., 2021). Lagrangian particle tracing models are useful for determining sources, trajectories, and destinations of drifting water parcels, with high temporal and spatial resolution. Thus, this approach could compensate for some of the limitations of satellite ocean color. However, this approach has so far been used to track particle locations along limited trajectories, preventing a synoptic view of the variability of water properties associated with HABs, e.g., algal concentration. Recently, a new methodology was proposed for synthesizing ocean color data (or in situ observations) with ocean circulation velocity fields from an operational model (Jönsson et al., 2009;Jönsson and Salisbury, 2016). This method could significantly improve our capability for monitoring of HABs.
Therefore, this paper aims to expand the Jönsson et al. (2009) scheme for synthesizing satellite observations and Lagrangian data to include extrapolation for better interpretation of the development of HABs. This improved scheme can not only fill gaps in HAB patches in the satellite images captured on cloudy days but also provides an early warning of harmful algal risk. An application example is presented showing the development of a Karenia mikimotoi bloom along the southern coast of England.

Satellite Data
In this study, the prediction scheme is based on the algorithm of Jönsson and Salisbury (2016), which combines simulated velocity fields with ocean color observations to create prediction of biological production in spatial and temporal scales. We applied this prediction scheme to reprocessed Sentinel-3A OLCI Level-3 products, which were downloaded from Plymouth Marine Laboratory (PML) ocean color archive. These products include chlorophyll-a (Chl) and remote sensing reflectance (R rs , sr −1 ) at 400,443,490,560,620,665,681,709,885, and 1020 nm with a spatial resolution of 300 m (Tilstone et al., 2020). The region of interest covers coastal and offshore areas in the southeastern England including the Celtic Sea and English Channel. The satellite passes the study region at around 12:00 a.m. (local time). Two case studies will be presented, covering September 10-15, 2019 and June 29 to July 5, 2019.

Lagrangian Particle Model
The particle tracking model PyLAG (Uncles et al., 2019) is used to produce the Lagrangian trajectories. PyLAG uses a fourth order Runge-Kutta scheme to advect particles, numerically integrated over a 100 s timestep. PyLAG is forced using hourly output from a hydrodynamic model with the horizontal turbulence statistics from the same model used to parameterise the diffusion term as random displacements.
The hydrodynamic forcing is from an operational setup of the Finite Volume Community Ocean Model (FVCOM) (Chen et al., 2003). This solves the prognostic equations on an unstructured grid, allowing higher resolution around complex coastlines and bathymetry, and lower resolution in the open ocean. Horizontal mixing is parameterised through the Smagorinsky scheme (Smagorinsky, 1963) and vertical turbulent mixing is modeled with the General Ocean Turbulence Model (GOTM) using a κ−ω formulation (Umlauf and Burchard, 2005). Lateral boundary conditions for the model are taken from the CMEMS North West Shelf data product (AMM7) and surface forcing from a Weather Forecast and Research (WRF) model which downscales output from the NOAA GFS global forecast model to provide 6-hourly forcing at 1 km resolution over the hydrodynamic model domain. Output of temperature and precipitation from the WRF model is also used to drive a neural network model to provide forecast river flows.

Merging Satellite and Particle Tracking
The algorithm that merges satellite products and velocity fields using Lagrangian particle tracking is described in previous studies (Jönsson et al., 2009(Jönsson et al., , 2011Jönsson and Salisbury, 2016). Here we summarise the main steps as follows (flowchart shown in Figure 1). Satellite products are firstly remapped to a uniform grid with a spatial resolution of 300 m. Virtual particles are then seeded randomly in each grid and advected for 7-10 days at hourly intervals using PyLAG. Next, satellite data are attached to the particle trajectories for the first 5-7 days, where the time difference between ocean color measurements and particle trajectories is limited to 30 min. The attached values are set to missing if satellite data are unavailable. Any particles leaving the model domain are removed to avoid errors. Further, an extrapolation procedure is conducted to predict the values of the satellite data (e.g., Chl) associated with each particle for the remainder of the particle drift, 2-3 days. The extrapolation is similar to the interpolation procedure in Jönsson and Salisbury (2016), but expanded to employ a linear extrapolation to calculate FIGURE 1 | Flowchart of ocean color prediction algorithm for water properties and early warning maps of HAB using Lagrangian particle trajectories and ocean color observations. a prediction based on any two valid satellite matchups. Note that the extrapolation will be stopped if only one valid matchup is obtained. In most cases there can be more than two matchups, thus an average value of all these extrapolations (v) is calculated using a weighting function of 1 divided by the days offset according to the following equation: where v i is the extrapolation value of the ith pair, n is the number of extrapolations, and d is the days offset.

Prediction of Remote Sensing Reflectance
The optical spectra of water bodies provide useful information on its constituents. Hence, the above approach is revised for the capability to predict R rs . It is worth noting that R rs is an apparent optical property (AOP), which highly depends on the light distribution. It would be questionable if R rs is estimated by linear weighting from extrapolations of all trajectories, as the light distribution for the particles along the trajectories could vary significantly. In this study, we propose a revised approach for prediction of R rs based on inherent optical properties (IOPs), which solely depend on optically properties of constituents in the water but are not affected by the changes of light field. The details are described by following steps. In the first step, the quasi-analytical algorithm (QAA) is employed to retrieve IOPs: absorption (a) and backscattering (b b ) coefficients of water (Lee et al., 2002). Then a new set of IOPs is predicted from the retrieved a and b b . Next, the predicted a and b b are used to reconstruct R rs with the following equation (Gordon et al., 1988).
where g 0 = 0.0949 sr −1 , g 1 = 0.0794 sr −1 , and R rs can be converted from r rs with: FIGURE 2 | Particle trajectories during September 10-17, 2019, which were used for an application of the prediction scheme.
Following these steps, we are thus provided with time-series images of Chl and R rs at 30-min intervals for the domain. These predictions are important as they will be used by algorithms to determine a quantitative estimate of HAB risk: this will be the focus of a future study. However, here we use R rs to predict the ocean color to enable visual forecasts of bloom development.

Animation of Map Sequences
To demonstrate water property changes over time frames, we combine individual images of Chl by interpolating intermediate timesteps for a smoother appearance using the Matplotlib Animation Python package. 1 To further diagnose more details on various water types, images were composited with redgreen-blue (RGB) bands from R rs (560), R rs (490), and R rs (443), respectively, and animated in an analogous way. These wavelengths cover the blue-green section of the optical spectrum, within which most ocean color variability is exhibited; hence this provides an enhanced view of the ocean rather than the true color.

RESULTS
The applicability of this prediction scheme is demonstrated by analyzing the predictions of water properties in two applications in the coastal and offshore regions of southeastern England.

The Celtic Sea and English Channel
Six days of ocean color data were used for this case study (September 9-15, 2019). A total of ∼1 million particles were seeded randomly and advected at hourly intervals for 8 days from September 9 to 17, 2019, including 48 h of predictions following the satellite period, at the same spatial resolution as the ocean color data (300 m × 300 m). Figure 2 shows a map of all modeled particle trajectories over the 8 days: the path of each particle is represented by a random colored line, which overlap each other Frontiers in Marine Science | www.frontiersin.org due to the high density of particles. Figure 3 shows maps of Sentinel-3 OLCI Level-3 Chl, from which we can identify the variable availability of ocean color during the 6 days due to cloud cover and orbit trajectories. There are some missing pixels in the study region but for most regions there was coverage for at least 50% of the time, e.g., for the region of English Channel, there were 3 days of ocean color observations available (September 9, 10, and 15) out of the 6 days. It is worth mentioning that from these maps we can observe that some algal blooms occurred along the coastal regions of England with high Chl (>∼4 mg/m 3 ). Figure 4 presents the time series predictions of Chl that shows the advection of surface waters during the first of the 2-day forecast. The most dynamic feature in this sequence is the movement of an algal bloom along the south coast of England, highlighted within the fixed ellipse in Figure 4. The value of this approach can be fully appreciated by viewing the animation of the sequence of predicted Chl maps that accompanies the figure. Advection effects are stronger in regions of the eastern English Channel with significant Chl variabilities over time, consistent with stronger tides in this region. The coastal waters generally demonstrate sharper Chl gradients than open ocean water, which is expected when considering tide, river runoff, and other forces dominating the shelf. The time series and animation of the predicted enhanced ocean-color results is shown in Figure 5: the colors in the maps indicate different water types, e.g., yellow colors in river estuaries may indicate high concentration of  sediments, and dark brown colors in the eastern English Channel could be due to a suspected K. mikimotoi bloom.
To further demonstrate the efficiency of the prediction approach, the conditions of water quality were also studied by plotting changes over time. Three transects were selected along latitudes 50 • N, 51.3 • N, and 52.5 • N (Figure 6), covering possible bloom regions (e.g., estuaries, south-east England coasts). Figure 7 shows the time series and animation of predicted Chl along the three transects. It was found that Chl varied over a wide range each day, usually with a periodic pattern mainly due to different stages of the tidal cycle that influence the advection of water parcels. In Figure 7B, the upward slope indicates a net westward advection of water along the transect.
An error estimation for the prediction algorithm was performed by comparing the predicted Chl on the second forecast day (September 17, 2019 12:00 p.m.) against the satellite-observed Chl scene that day (12:29 p.m.). The scatter plot of cloud-free pixels in the two datasets are shown in Figure 8. The predicted Chl agrees well with satellite Chl with an r 2 of 0.81 and a

High Spatial Resolution in English Channel
Monitoring of water bodies at high spatial resolution is of paramount importance, especially in dynamic coastal waters, where the water constituents can vary dramatically over a small distance. To evaluate the ability of the prediction algorithm for high spatial resolution, a further period was studied in the English Channel, June 29 to July 5, 2019 when a large coccolithophore bloom was observed. In this case, many more particles were seeded (∼10 million particles). To demonstrate the techniques using computer time feasible for an operational monitoring system, the region of interest was limited to the English Channel and only 16 h of predictions were generated. The spatial resolution of the resulting prediction images is 100 m with a temporal resolution of 1 h. Figure 9 shows and animates the changes of Chl distributions over time. Small scale structures are interpreted here (e.g., some fine structures of algae patches were revealed). To gain a precise understanding of water types, the enhanced ocean-color images with high spatial  resolution are shown and animated in Figure 10. The (harmless) coccolithophorid bloom occurred in the northern English Channel, covering thousands of square kilometers with milky blue. The prediction method here successfully discriminated the locations of the bloom and its advection over time, further demonstrating the value of using this prediction method for monitoring algal blooms. FIGURE 10 | Frames of predicted enhanced ocean-color RGB in the English Channel with high spatial resolution (100 m) showing more details of water properties including the movement of a bright coccolithophore bloom near 53 • N (animation: https://rsg.pml.ac.uk/shared_files/junl/paper_animation/example2/rgb.gif).

DISCUSSION
It is quite clear, as seen both in our findings and in earlier studies (e.g., Jönsson et al., 2009;Jönsson and Salisbury, 2016), that biological processes in coastal are strongly affected by physical advection of water parcels, and that combining satellite derived products and simulated velocity fields using particle tracking provides a powerful approach for such analyses. While earlier works by Jönsson et al. (2011) and Jönsson and Salisbury (2016) are primarily focused on estimating the rate of change in satellite properties, we present results showing that the method can be expanded to predict locations and concentrations of the properties. This novel application can be used to better predict HAB events and to identify the extent and location of blooms more precisely, especially in regions with strong tidal cycles where surface currents are predictable. We find our results to be generally physically coherent and that any patches of high phytoplankton biomass are predicted and advected in a reasonable way. The nature of the method is such that we expect significant methodological errors to be seen as noise or spurious changes in our predicted fields (see Jönsson et al., 2009 for further discussion). We see for example, as expected, generally sharper Chl gradients in coastal waters compared to the open ocean in the Celtic Sea and English Channel. The method is also able to successfully predict the location of a coccolithophore bloom. The skill of the method discussed here is in line with what earlier studies have reported from the California Current (Jönsson and Salisbury, 2016) and Gulf of Maine (Jönsson et al., 2009) when estimating rates of change, which suggests that the extrapolation approach has a similar ability to provide useful information.
While the prediction tool presented in this study has the potential to improve regional HAB monitoring by more precisely assessing how the current state will develop in the near future, we readily agree that there are limitations. It is inevitable that there will be inaccuracies due to imperfections in the modeling and methodology. Errors in the prediction method could arise from many aspects, e.g., imperfect velocity of fields, extrapolation of model velocities in time and space, numerical integration of the trajectory path, and the omission of vertical velocity.
In this study, Lagrangian particle trajectories were computed according to the effect of both advection and diffusion. However, because each trajectory only contains a limited number of particles, the final locations of simulated particles could not be representative of all possible particle locations due to diffusion. The problem could be especially acute for some areas with large horizontal current shear. Thus, it would be important to keep in mind that uncertainty from diffusion always exists in these areas. The problems could be solved by seeding many more particles at each initial location, though the additional burden of computing their trajectories is currently unfeasible for this real-time application. A practical solution could be to introduce a probability density function indicating possible locations for each traced particle. However, it is beyond the scope of this study and substantial studies of this question are anticipated in the future. Regarding errors from advection, this was assessed in the previous study of Jönsson et al. (2009) by comparing the changes in Chl and sea surface temperature (SST) between the start and end positions of particle trajectories advected between two satellite images. The results showed that the advective errors are small relative to total changes due to other factors.
The extrapolation procedure will also introduce errors to the predictions. The predicted values were estimated via linear extrapolation with a weighting function. To reduce error from the extrapolation, when the extrapolated value for a trajectory exceeded 50% of the mean value of all extrapolations, this trajectory was disregarded and was not included for the calculation of final prediction. Considering uncertainty from the IOPs retrieval algorithm, errors in predicted IOPs are inevitable, but errors in R rs could possibly be compensated by reconstruction from these IOPs. Furthermore, it is also worth noting that although individual trajectories have errors, the statistics obtained from thousands of particles are still very informative.
Still, while not perfect, the resulting predictions provide an enhanced set of information for managers and stakeholders to include when assessing the risk of HABs that has not been available until now. The prediction tool is also agnostic to any of its component modules. We are able to easily leverage new improved ocean circulation models and satellite products into the framework to provide as accurate predictions as possible.

CONCLUSION
This study set out to develop a prediction scheme for monitoring HABs by merging satellite observations and Lagrangian particle tracking. Two case studies in regions along the coast of England have shown that the Lagrangian methodology is effective for the interpretation of satellite data for early warning of HAB risk. The accuracy of the predictions relies on many factors. Particularly, the uncertainties from modeling and methodology, e.g., velocity of fields, extrapolation of model velocities in time and space, numerical integration of the trajectory path, and the omission of vertical velocity, may have big impacts on determining the accuracy of final predictions. Notwithstanding these imperfections, this work offers a powerful tool for monitoring and observing the state of the marine ecosystem. The synoptic, time-resolved quantification is invaluable to our understanding of HABs developments. Animated sequences generated by this method promote greater understanding and usage of satellite ocean color data for communicating with aquaculture farmers. Future research will extend the approach to predict quantitative risk maps for key high-biomass HAB species. Hence our future priority is to seek further applications of this new technique to support the aquaculture industry with improved early warning of potential HABs.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JL developed the prediction scheme, implemented the code, and wrote the first draft of the manuscript. PM proposed initial ideas and helped shape the research. BJ provided suggestions on code and algorithm development. MB helped carry out numerical simulations. All authors discussed the results and commented on the manuscript.

FUNDING
This study was supported by the Interreg Atlantic Area Programme project "Predicting the Impact of Regional Scale Events on the Aquaculture Sector" (PRIMROSE, Grant No. EAPA_182/2016) and UK BBSRC project "Safe and Sustainable Shellfish: Introducing local testing and management solutions" (reference BB/S004211/1). BJ was in part funded by NASA 80NSSC21K0565 and a Simons Foundation grant (549947, SS).