First Results of Phytoplankton Spatial Dynamics in Two NW-Mediterranean Bays from Chlorophyll-a Estimates Using Sentinel 2: Potential Implications for Aquaculture

: Shellﬁsh aquaculture has a major socioeconomic impact on coastal areas, thus it is necessary to develop support tools for its management. In this sense, phytoplankton monitoring is crucial, as it is the main source of food for shellﬁsh farming. The aim of this study was to assess the applicability of Sentinel 2 multispectral imagery (MSI) to monitor the phytoplankton biomass at Ebro Delta bays and to assess its potential as a tool for shellﬁsh management. In situ chlorophyll-a data from Ebro Delta bays (NE Spain) were coupled with several band combination and band ratio spectral indices derived from Sentinel 2A levels 1C and 2A for time-series mapping. The best results (AIC = 72.17, APD < 10%, and MAE < 0.7 mg / m 3 ) were obtained with a simple blue-to-green ratio applied over Rayleigh corrected images. Sentinel 2–derived maps provided coverage of the farm sites at both bays allowing relating the spatiotemporal distribution of phytoplankton with the environmental forcing under di ﬀ erent states of the bays. The applied methodology will be further improved but the results show the potential of using Sentinel 2 MSI imagery as a tool for assessing phytoplankton spatiotemporal dynamics and to encourage better future practices in the management of the aquaculture in Ebro Delta bays.


Introduction
Shellfish are filter-feeding organisms that feed on different types of suspended particles in the water column, thus their production is mainly related to phytoplankton availability [1].Spain is the leading producer and consumer of bivalves in Europe, Catalonia being the most important producer area in the Spanish Mediterranean, with most of the production concentrated in the Ebro Delta (Figure 1).The most important species for aquaculture are the Mediterranean mussel (Mytilus galloprovincialis) and the Pacific oyster (Crassostrea gigas), but other bivalves such as clams (e.g., Ruditapes philippinarum) and cockles (e.g., Venus verrucosa) are also harvested.Bivalve culture is mainly developed inside its two bays, Alfacs and Fangar, representing 1.8% and 6.5% of their respective surfaces [2].Since 1990, an official monitoring program carried out by the Regional Government of Catalonia establishes a weekly analysis of the phytoplankton community and water physicochemical parameters at different locations of both bays (12 samples per week).However, the sampling procedure is temporally and spatially limited, so global extrapolations are subject to large uncertainties.
Temporal phytoplankton dynamics are highly influenced by the nutrient input from rice fields trough the irrigation network [3].Furthermore, freshwater inputs have a great physicochemical impact in both bays, increasing water column stratification and dominating over wind on a seasonal scale [4].Therefore, freshwater input imposes a double layer circulation system like typical estuarine circulation patterns.However, when channels are closed (from October to April), the water renewal time of the bays increases, forming retention areas that can become accumulation zones.Both scenarios may favor phytoplankton growth.On shorter time scales (days to weeks), the wind is the main controlling factor of water mixing [5] by breaking the vertical stratification.Therefore, water circulation patterns, and hence phytoplankton temporal and spatial variability, depend on freshwater inputs, meteorology, and coastal geomorphology [6].Remote sensing allows obtaining information of marine and continental processes at different spatiotemporal scales [7].Chlorophyll-a (chl-a) is the main photosynthetic pigment present in algae and an optically active seawater constituent; thus, it is commonly used as indicator of phytoplankton biomass and has significant implications on remote sensing [8][9][10][11].The estimation of chl-a concentration from remotely sensed data requires the development of algorithms with a maximal sensitivity to chl-a and minimal to the rest of constituents present in the water [12].Different authors have proposed several methodologies to estimate chl-a from satellite remote sensing imagery (see a review in [13][14][15]); for instance, a classical approach is developing relationships between band-ratios (namely color indices) or their combinations [14].Several ratio-based and 3-band combination algorithms have been proposed, including the common Blue to Green ratios, the Ocean Color-based algorithms [16,17], and those including the Red edge [18][19][20][21], which take advantage of pigment's absorption maxima (i.e., at 665 nm) [22,23].Other approximations are based on spectral band difference by using band triplets from the Red and Near Infrared (NIR), such as the Fluorescence Line Height (FLH) [24], the Maximum Chlorophyll Index (MCI) [24], and the Maximum Peak Height (MPH) [25].The properties of coastal waters, however, are controlled by complex interactions and fluxes of material between land, ocean, and atmosphere, which makes challenging to achieve reasonable estimates of water-leaving radiance (removing atmospheric contributions from a signal received at the TOA), and to obtain a robust relationship between water quality and satellite-based parameters [26] (integrating the remote sensing and in situ measurements).Although a large amount of satellite data is available for remote sensing of chl-a (e.g., SeaWiFS, MODIS, MERIS), the fast dynamics of phytoplankton in coastal areas, both temporally and spatially, cannot be fully resolved because of either their coarse spectral, spatial and/or temporal resolution.Currently, the increased frequency (up to five days under ideal conditions) and higher spatial resolution (10 to 60 m 2 ) of Sentinel 2 together with its spectral band configuration has opened a new potential to remote sensing of chl-a in coastal zones of small geographical extension, and hence as an alternative for phytoplankton monitoring in coastal areas.
The overall purpose of this study was to analyze the potential of Sentinel 2 multispectral imagery (MSI) data as a support tool for the future management of shellfish cultures through the monitoring of phytoplankton biomass in the Ebro Delta bays.Thus, this paper is a first attempt to assess chl-a in a shallow coastal environment with Sentinel 2 imagery data, a free public resource.The objectives of this study were to 1.
Generate 20 m 2 resolution chl-a maps from Sentinel-2 MSI imagery covering the whole system; 2.
Understand the spatiotemporal phytoplankton biomass dynamics by using the derived chl-a maps and to relate them to environmental variables and the rice farming year; 3.
Assess the applicability of the results to shellfish aquaculture management in the area.

Study Sites
The Ebro Delta is one of the largest (320 km 2 ) deltas in the northwestern Mediterranean Basin.The climate is Mediterranean temperate with warm dry summers and cool wet winters, annual mean temperature ranges between 5 and 33 • C, and annual precipitation from 500 to 600 mm, being maximum in autumn and minimum in summer.The study area was located in the two bays of the Ebro Delta (Figure 1).Fangar Bay, with an area of 12 km 2 , is connected to the sea by a narrow mouth (ca. 1 km wide) that is currently closing.Maximum depth is ca. 4 m, which makes it very sensitive to environmental variations.Water temperature ranges between 6.5 and 32 • C, salinity varies from 9 to 37 PSU (Practical Salinity Unit) and renewal time is about four days when channels are open [6].Alfacs Bay, covering an area of 56 km 2 and connected to the sea by a channel of 2.5 km wide, has an average depth of 3.13 m (maximum depth is 7 m).Water temperature ranges between 8 and 32 • C, the salinity varies from 26 to 37 PSU, and the renewal time is about 15 days with open channels [6].The hydrology of both bays is highly influenced by freshwater inputs from the irrigation network (Figure 1).Freshwater and nutrient inputs from the river allowed the development of prosperous fishery and farming activities.The production of bivalves in Ebro Delta bays constitutes a major economic activity in the area (Figure 1), together with agriculture, since 210 km 2 of the delta plain are devoted to rice production (Figure 1).Ebro Delta (Figure 1).Fangar Bay, with an area of 12 km 2 , is connected to the sea by a narrow mouth (ca. 1 km wide) that is currently closing.Maximum depth is ca. 4 m, which makes it very sensitive to environmental variations.Water temperature ranges between 6.5 and 32 °C, salinity varies from 9 to 37 PSU (Practical Salinity Unit) and renewal time is about four days when channels are open [6].Alfacs Bay, covering an area of 56 km 2 and connected to the sea by a channel of 2.5 km wide, has an average depth of 3.13 m (maximum depth is 7 m).Water temperature ranges between 8 and 32 °C, the salinity varies from 26 to 37 PSU, and the renewal time is about 15 days with open channels [6].
The hydrology of both bays is highly influenced by freshwater inputs from the irrigation network (Figure 1).Freshwater and nutrient inputs from the river allowed the development of prosperous fishery and farming activities.The production of bivalves in Ebro Delta bays constitutes a major economic activity in the area (Figure 1), together with agriculture, since 210 km 2 of the delta plain are devoted to rice production (Figure 1).

In situ Data: chl-a
Eight water samplings campaigns were carried out from April 2016 to August 2017 coinciding with the Sentinel 2A satellite pass (Table 1).Different sampling grids were used (see Table 1) for different days, and not both of the bays were sampled every day.Integrated water samples were collected using the Lindahl methodology [27] (N = 106).In addition, on 25 July 2017 and 4 August 2017, surface water samples were collected (N = 40) with polypropylene bottles.Seawater samples were kept in a portable cool box until arrival to the laboratory.In the laboratory, three different methods were used to measure chl-a concentration, in vivo fluorimetry [28] (hereafter in vivo), and after acetone extraction both in a fluorometer (corrected chl-a; hereafter FL) [29], and in a spectrophotometer (chl-a; hereafter SP) [30].For all samples (N = 106), chl-a was estimated in vivo, and in 58 of them, chl-a was measured after acetone extraction.Briefly, water samples (550-1000 mL) were filtered using fiberglass filters (GF/F), and filters were submerged in 10 ml of acetone inside 15 mL labelled conical centrifuge tubes.After 24 hours in the fridge (4 °C), they were sonicated for 5

In situ Data: chl-a
Eight water samplings campaigns were carried out from April 2016 to August 2017 coinciding with the Sentinel 2A satellite pass (Table 1).Different sampling grids were used (see Table 1) for different days, and not both of the bays were sampled every day.Integrated water samples were collected using the Lindahl methodology [27] (N = 106).In addition, on 25 July 2017 and 4 August 2017, surface water samples were collected (N = 40) with polypropylene bottles.Seawater samples were kept in a portable cool box until arrival to the laboratory.In the laboratory, three different methods were used to measure chl-a concentration, in vivo fluorimetry [28] (hereafter in vivo), and after acetone extraction both in a fluorometer (corrected chl-a; hereafter FL) [29], and in a spectrophotometer (chl-a; hereafter SP) [30].For all samples (N = 106), chl-a was estimated in vivo, and in 58 of them, chl-a was measured after acetone extraction.Briefly, water samples (550-1000 mL) were filtered using fiberglass filters (GF/F), and filters were submerged in 10 mL of acetone inside 15 mL labelled conical centrifuge tubes.After 24 h in the fridge (4 • C), they were sonicated for 5 min (ultrasonic processor) and centrifuged for 10 min at 4000 rpm at 4 • C. The chl-a concentration was then measured in a SHIMADZU UV-1800 spectrophotometer (Shimadzu Corporation, Kyoto, Japan) and/or in a TURNER Trilogy ® fluorometer (Turner Designs, San Jose, CA, USA) (Table 1).The datasets generated during the current study are available from margarita.fernandez@irta.caton reasonable request.2a,b).Grid 2: specific sampling grid designed for ground truth of Sentinel images (green dots in Figure 2a,b).Grid 3: sampling grid of the project "Model of water circulation in Fangar Bay from the European Maritime and Fisheries Fund (EMFF)" (white dots in Figure 2b).

Sentinel 2 Data
A set of 47 Sentinel 2A L1C images (i.e., not cloud covered) were downloaded from Copernicus Open Acces Hub (https://scihub.copernicus.eu/).Thirteen of them (six from Alfacs and seven from Fangar) within the period April 2016-August 2017 were selected for calibration and validation purposes (Figure 3).The remaining images between January 2017 and January 2018 (Table 1), 18 from Alfacs and 16 from Fangar, were used for time-series estimation (Figure 3).Although the calibration/validation (CalVal) image sets covered mainly spring and summer, the time-series was estimated for a full year in order to include the full rice growing season.Additional meteorological data, including daily air temperature ( • C), wind direction ( • ), wind speed (m/s), and precipitation (mm), were obtained from the Illa de Buda meteorological station (Station Id. 11043, located at 1 m above sea level) of the Catalan Meteorological Service, http://www.meteo.cat(Figure 1).

Atmospheric Correction: ACOLITE
Sentinel 2A L1C imagery were atmospherically corrected with ACOLITE processor.It bundles the atmospheric correction algorithms and processing software developed by the Royal Belgian Institute of Natural Sciences (RBINS) for aquatic applications of Landsat (5/7/8) and Sentinel 2 (A/B) satellite data.The Dark Spectrum Fitting (DSF) [31], used here, computes the atmospheric path

Atmospheric Correction: ACOLITE
Sentinel 2A L1C imagery were atmospherically corrected with ACOLITE processor.It bundles the atmospheric correction algorithms and processing software developed by the Royal Belgian Institute of Natural Sciences (RBINS) for aquatic applications of Landsat (5/7/8) and Sentinel 2 (A/B) satellite data.The Dark Spectrum Fitting (DSF) [31], used here, computes the atmospheric path

Atmospheric Correction: ACOLITE
Sentinel 2A L1C imagery were atmospherically corrected with ACOLITE processor.It bundles the atmospheric correction algorithms and processing software developed by the Royal Belgian Institute of Natural Sciences (RBINS) for aquatic applications of Landsat (5/7/8) and Sentinel 2 (A/B) satellite data.The Dark Spectrum Fitting (DSF) [31], used here, computes the atmospheric path radiance based on multiple targets in the scene or sub-scene, with no a priori dark band, allowing an aerosol correction.ACOLITE includes a sun glint correction, which uses the short-wave infra-red (SWIR) bands to estimate a glint signal [32] and to establish the threshold to determine which pixels need to be corrected (0.05 by default).Sentinel 2A B11 and B12 bands (SWIR at 1604 nm and 2202 nm) were used for sun glint correction.The thresholds were set manually image-by-image after a SWIR analysis that was carried out considering the response of Sentinel 2A B11 over water pixels compared to non-water pixels.For each day and bay land/water mask and sun-glint correction thresholds were defined, ranging between 0.0215 and 0.1.Therefore, the atmospheric correction procedure output included, for each image, uncorrected (a), partially corrected (b), and fully corrected atmosphere (c and d) reflectance data.
(a) Rhot: top of atmosphere reflectance (TOA) derived from the original input file.(b) Rhorc: Rayleigh corrected reflectance.This is the Rhot with removed and corrected reflectance for two-way Rayleigh transmittance.An additional pre-processing step was made to avoid high reflectance pixels by fixing a maximum threshold (Rhorc reflectance at 492 nm or 560 nm ≥ 0.11) above which pixels were assigned as invalids.(c) Rrs: remote sensing reflectance (sr −1 ) for water pixels (Rrs = Rhow/π).(d) Rhow: surface reflectance for water pixels.

Chlorophyll-a Estimation Algorithms
Seven different spectral algorithms band-ratio and band-combination based, were applied to each product resulting of ACOLITE processing (Rhot, Rhorc, Rrs and Rhow).Briefly, I.
BG: The Ratio between Blue and Green spectra uses the reflectance at 490 nm (blue) and 560 nm (green).At 490 nm carotenoids absorb light strongly, while at 560 nm the absorption of all photosynthetic pigments is minimal (i.e., green reflection).This algorithm was initially proposed by [33].R stands for Rhot, Rhorc, Rrs, or Rhow reflectance. [ II. BR: The Blue-Red ratio is based on the two chl-a maximal absorption peaks. [ III. RG: The Green-Red ratio is based on the minimal and maximal absorption peaks of chl-a, thus avoiding the use of the blue bands [23,34].
IV. NR: The ratio between Red and NIR assumes that the absorption by non-algal particles, yellow substances and the backscattering are insignificant when compared to chl-a absorption at red wavelengths (665 nm) [35].Between 700 and 720 nm, the absorption due to water constituents is minimal. [ V. NDCI: The Normal Difference Chlorophyll Index developed by Mishra et al. [36] for turbid productive waters uses the information of the reflectance peak centered at 700 nm, which is maximally sensitive to variations in chl-a concentration in water.Furthermore, a wide spectral absorption peak between 665 nm and 675 nm is generally associated to the absorption by chl-a pigments.The normalization through the NDCI eliminates uncertainties in the estimation of the remote sensing reflectance, seasonal solar azimuth differences, and atmospheric contributions.
VI. DO5: Dall'Olmo and Gitelson [37] presented a three-band model using Red and NIR bands.It assumes that (i) the absorption by coloured dissolved organic matter (CDOM) and total suspended matter (TSM) beyond 700 nm is approximately equal to that at 665-675 nm and the difference between them can be neglected; (ii) total chlorophyll, CDOM, and TSM absorption beyond 730 nm is almost 0; and (iii) backscattering coefficient of chl-a is spectrally invariant [36,37].
VII. MCI: The Maximum Chlorophyll Index allows the detection of red tides and aquatic vegetation [24].For Sentinel 2, it uses the band 5 (705 nm), perfectly located to detect high biomass water bodies against the baselines of the bands 4 and 6 (665 and 740 nm).
In Equation ( 7), k is the thin cloud correction factor fixed at 1.005 for thin clouds. [

Model Calibration and Validation
Sentinel 2A (Level 1C and 2) images and all in situ chl-a of coinciding days were used for model calibration and validation.In order to reduce the effect of noise from the sensor and the time-difference between the image (20 m 2 resolution) and water samples acquisition, reflectance was averaged over a 3 × 3 pixel-box centered at the in situ measurements.However, not all of the nine pixels per in situ sampling location could be used as there might be outliers coming from different sources such as bottom contamination, different affection of sun glint and adjacency or infrastructures as rafts or harbor jetties interfering in some pixels.For this reason, a pre-processing step was carried out on each spectral band used and for all atmospheric correction levels.For each day and bay, considering together all in situ sampling locations, outliers were detected and removed by Tukey's fences method (Boxplot).The criteria flagged as invalid a pixel if in one of the five spectral bands (see Equations ( 1)-( 7)) the reflectance value was classified as outlier.A second step to clean the remaining outliers was carried out applying the same methodology to each 3 × 3 pixel-box centered at in situ sampling sites, individually.To ensure the possibility of using the averaged reflectance of 2-9 pixels, without corrupting the methodology, standard deviation (SD) of the average at each sampling location was computed against the number of pixels used for the average.
After outlier deletion, the seven algorithms were computed using the averaged reflectance at each chl-a sampling location.Model calibration was done with 70% of the data (with ordinary least of squares fitting, OLS) and the remaining 30% was used for model validation.Models were calibrated and validated in two different ways: (i) using only those samples where in situ chl-a was measured by the three methodologies (i.e., in vivo, FL, SP) and (ii) for each methodology including all the available data.In both cases, model development was carried out considering all possible combinations of ACOLITE-derived imagery together with two different scenarios (individually or both bays together).
Model performance was assessed graphically by plotting observed and predicted values, and efficiency was measured with the Akaike Information Criterion (AIC), the Averaged Percentage Difference (APD), and the Mean Absolut Error (MAE).AIC combines fit and parsimony (number of parameters) of models, with the best fitting model having the lowest AIC.MAE and APD were applied following the criteria of [38], who suggested that these metrics account better for accuracy of the models over non-Gaussian distributions by not amplifying outliers and precisely reflecting the error magnitude.Models with lowest AIC, MAE, and APD, in this order, were considered better.Although the coefficient of determination (i.e., R 2 ) and Normalized Root Mean Squared Error (RMSE) are widely used goodness-of-fit measures, they are not recommended for non-Gaussian distributions [38].Thus, both measures were only included to allow the comparison with previous works.

Time-Series Estimation
The best model was selected to construct chl-a time-series maps with the available Sentinel 2A images in 2017.Pixel-stability was assessed by using an unsupervised classification cluster analysis (2 classes) based on the inter-pixel slope of the averaged time-series chl-a and the coefficient of variation (CV; Equation ( 8)) of chl-a of the same set of images.
where σ stands for standard deviation (SD) and X for the average.

Workflow
The proposed workflow (Figure 4) started with the selection and download of Sentinel 2A L1C images.The images were processed with ACOLITE after the SWIR analysis, including a resampling of all bands to 20 m 2 , image cropping to the region of interest (Ebro Delta bays), and the atmospheric correction to derive Rhorc, Rrs, and Rhow reflectance (see Section 2.4).After ACOLITE processing, for the spectral bands of interest, outliers were detected and removed.Then, for each image of the calibration set, the spectral algorithms were computed, and the resulting values were extracted at each chl-a sampling location.Models were calibrated and validated with Rhot, Rhorc, Rrs, and Rhow imagery together with ground truth data.The best algorithm and methodology were selected, applied to all the available images in 2017, and the pixel-stability analysis was carried out.Finally, the resulting time series was then used to analyze spatiotemporal patterns of chl-a (as indicator of phytoplankton biomass dynamics), thus covering different seasons and the full rice farming cycle.

In situ Data: chl-a
Overall, chlorophyll-a concentration varied among seasons and sites, with different spatial distribution patterns in both bays.In Alfacs Bay, chl-a showed a spatial gradient trend defined All statistical analyses were performed with R version 3.5.2; the packages foreign 0.8.71, xlsx 0.6.1,xlsxjars 0.6.1,ncdf 1.16.1, and raster 2.8.19 were used to load external data with different formats.Packages rgdal 1.4.3,spatstat 1.59.0, and maptools 0.9.5 were used to work with geospatial data (create masks, band math calculator, and pixel extraction).Packages FSA 0.8.24,NCStats 0.4.7,nlstools 1.0.2, and minipack.lm1.2.1 were used to evaluate the model performance (ROC curves and associated statistical parameters).

In situ Data: chl-a
Overall, chlorophyll-a concentration varied among seasons and sites, with different spatial distribution patterns in both bays.In Alfacs Bay, chl-a showed a spatial gradient trend defined generally by higher concentrations from the central zone with higher concentration values, to the inner area, with minimum chl-a concentrations in the shellfish rafts (Figure 1).In Fangar Bay, maximum chl-a concentrations were found in the mouth and minimum concentrations in the inner part of the bay, which showed similar values to those in the shellfish rafts.Table 2 summarize chl-a results per bay.In general, Alfacs Bay showed higher chl-a concentrations.
Among the different laboratory methodologies used to measure chl-a concentration, in vivo results showed moderate correlation values with both FL (Pearson's r = 0.60, N = 55, P < 0.001) and SP (r = 0.62, N = 43, P < 0.001), while these two methods (FL and SP) were highly correlated (r = 0.93, N = 43, P < 0.001).The average percentage difference (APD) between methodologies was 98% between in vivo and FL, 56% between in vivo and SP and 20% between FL and SP.Surface and integrated water column (sampled in Fangar Bay on both 25 July and 17 August) in vivo chl-a concentrations were strongly correlated (r = 0.80, N = 40, P < 0.001), with an APD of 7.6%.

Atmospheric Correction and Outlier Removal
The averaged reflectance at the sampling locations for the different atmospheric correction products (i.e., Rhot, Rhorc, Rhow, and Rrs) for each CalVal date and bay (Figure 5), and at each Sentinel 2A band, showed less reflectance from uncorrected to full corrected levels, this being more pronounced for shorter wavelengths.Fangar Bay showed higher averaged reflectance than Alfacs, when comparing the same day, and for all different Level products.
Averaged reflectance of a 3 × 3-pixel box centered at the in situ sampling points was used as the reflectance at each location; however, outlier pixels were removed.After outlier detection, 18 sampling points were completely removed and were not used in the CalVal process.Sixteen of the 18 removed points corresponded to Fangar Bay and were mostly located within the shellfish rafts, the mouth of the bay, and the inner area.Two points were removed from Alfacs Bay, both located in the harbor on 20 June 2016.Final available chl-a data are summarized in Table 3.In order to evaluate the impact of outlier pixels on the reflectance estimation, it was assessed the reflectance SD relative to the number of valid pixels (2 to 9), at each sampling site, for all type of Sentinel 2A products and all the bands used in algorithm calculation.Pearson's correlation coefficient, in absolute value, ranged between 6.11 × 10 −3 and 0.23, thus, reflectance values were similar, independent of the size of the pixel-box around the sampling point (from 2 to 9), and outliers can be removed without introducing significant errors.

Model Calibration and Validation
All variable combinations resulted in 252 models (see Table S1); chl-a methodology (in vivo, FL and SP) × bay (Alfacs, Fangar, and both bays together) × Sentinel 2A images (Rhot, Rhorc, Rrs and Rhow) × spectral algorithm (BG, BR, RG, NR, NDCI, DO5, MCI) (see Table S1).Overall, considering all the possible models, the algorithms performed better when applied to Rhorc images, although Red-to-Green (RG) and, especially MCI, showed less sensibility to the atmospheric correction and similar results were achieved with Rhot, Rhorc, Rrs, or Rhow reflectance.The best results were obtained combining Rhorc images with spectrophotometer chl-a measures (SP).Within the "Rhorc_SP" models, the best performing algorithms were BG (Blue-to-Green ratio) for Alfacs Bay and for both bays together, and the NDCI (Normal Difference Chlorophyll Index) algorithm returned the best results for Fangar Bay (Table 4).Close results to BG were achieved in Alfacs Bay and both bays together with RG (Red-to-Green) band ratio, while worse results in both cases were obtained with Maximum Chlorophyll Index (MCI).In Fangar Bay, despite differences among the performance of the different algorithms were smaller than in Alfacs Bay (Table S1), NIR-to-Red (NR) band ratio and MCI performed similar to NDCI, while BG performed worse.
Chlorophyll-a was not measured by the three methodologies (in vivo, FL, and SP) in all the

Model Calibration and Validation
All variable combinations resulted in 252 models (see Table S1); chl-a methodology (in vivo, FL and SP) × bay (Alfacs, Fangar, and both bays together) × Sentinel 2A images (Rhot, Rhorc, Rrs and Rhow) × spectral algorithm (BG, BR, RG, NR, NDCI, DO5, MCI) (see Table S1).Overall, considering all the possible models, the algorithms performed better when applied to Rhorc images, although Red-to-Green (RG) and, especially MCI, showed less sensibility to the atmospheric correction and similar results were achieved with Rhot, Rhorc, Rrs, or Rhow reflectance.The best results were obtained combining Rhorc images with spectrophotometer chl-a measures (SP).Within the "Rhorc_SP" models, the best performing algorithms were BG (Blue-to-Green ratio) for Alfacs Bay and for both bays together, and the NDCI (Normal Difference Chlorophyll Index) algorithm returned the best results for Fangar Bay (Table 4).Close results to BG were achieved in Alfacs Bay and both bays together with RG (Red-to-Green) band ratio, while worse results in both cases were obtained with Maximum Chlorophyll Index (MCI).In Fangar Bay, despite differences among the performance of the different algorithms were smaller than in Alfacs Bay (Table S1), NIR-to-Red (NR) band ratio and MCI performed similar to NDCI, while BG performed worse.
Chlorophyll-a was not measured by the three methodologies (in vivo, FL, and SP) in all the samples; thus, models had different sample size.In order to avoid the influence of sample size on results, the models were also fitted using only those chl-a samples measured by the three methodologies (Table S2).There were not significant changes associated to N, but changes on model performance were more related to the range of chl-a covered by the samples (e.g., the lower variability of chl-a in Fangar Bay).
Different algorithms performed better in Alfacs and Fangar Bay.The low number of available SP data and the good results obtained calibrating and validating the model including both bays together suggest the use of "Rhorc_SP" configuration (Figure 6) until more data are available.Despite BG performance in Fangar Bay was worse than the achieved with other algorithms (i.e., NDCI, NR, and MCI), probably it was due to the lack of variability towards higher concentrations and the weight of few extreme values over a small dataset.In fact, the linear distribution of chl-a SP in Fangar fit with the trend of Alfacs (Figure 7).Also, the trend line using data of both bays or using data only from Alfacs Bay was almost equal, denoting that Fangar samples were in agreement with the global trend described (Figure 7).These results support the idea of using both bays together and reinforce the assumptions for applying the same model to both bays.According to previous results, chl-a time-series was generated from the Blue-to-Green ratio (BG) model using partially corrected images (Rhorc) for both bays together, and chl-a measured by spectrophotometer (SP).Two pre-processing steps were applied to reduce the sources of error on the bottom of Rayleigh corrected reflectance.First, although images were selected according to cloud absence, in two of them, small areas at the extremes of Alfacs Bay were contaminated by clouds.There, the threshold applied to the Rhorc images removed pixels associated with thick clouds (Figure 8), but the thinnest clouds were not successfully detected, and the ground information was not restored in either case.The second pre-processing step consisted in the generation of a mask to avoid areas where BG did not responded only to chl-a, but to other sources such as bottom reflectance or macrophytes (Figure 9).The clustering used to make the mask highlighted the boundaries where maximum changes occurred, namely, shallow waters with bottom or seagrass contribution, hard structures such as rafts, and semi-static objects like the ships in the Alfacs Bay harbor.Finally, before chl-a time-series estimation, a 20 m buffer (i.e., 1 pixel) was applied, around each raft, created in order to delete mixed border pixels.According to previous results, chl-a time-series was generated from the Blue-to-Green ratio (BG) model using partially corrected images (Rhorc) for both bays together, and chl-a measured by spectrophotometer (SP).Two pre-processing steps were applied to reduce the sources of error on the bottom of Rayleigh corrected reflectance.First, although images were selected according to cloud absence, in two of them, small areas at the extremes of Alfacs Bay were contaminated by clouds.There, the threshold applied to the Rhorc images removed pixels associated with thick clouds (Figure 8), but the thinnest clouds were not successfully detected, and the ground information was not restored in either case.The second pre-processing step consisted in the generation of a mask to avoid areas where BG did not responded only to chl-a, but to other sources such as bottom reflectance or macrophytes (Figure 9).The clustering used to make the mask highlighted the boundaries where maximum changes occurred, namely, shallow waters with bottom or seagrass contribution, hard structures such as rafts, and semi-static objects like the ships in the Alfacs Bay harbor.Finally, before chl-a time-series estimation, a 20 m buffer (i.e., 1 pixel) was applied, around each raft, created in order to delete mixed border pixels.Both for Fangar and Alfacs bays, one-year chl-a time-series were processed (Figure S3).Overall, during winter and early spring, higher concentrations of chl-a were observed in Alfacs Bay.From April to October, chl-a concentrations were comparable between bays; after, in Fangar Bay, chl-a concentrations decreased more sharply.Despite the differences in chl-a concentration, the general trend was similar in both bays almost all the year, but chl-a peaks differed.In Alfacs Bay, chl-a peaked during February-April and October-November, achieving maximum concentrations in March.In Fangar Bay, chl-a peaked on May and September-October, being the most productive along the year in this last period.Minimum chl-a concentrations were found in winter in both bays, January and November for Alfacs and Fangar Bay, respectively.The coefficient of variation (CV) of chl-a (Figure Both for Fangar and Alfacs bays, one-year chl-a time-series were processed (Figure S1).Overall, during winter and early spring, higher concentrations of chl-a were observed in Alfacs Bay.From April to October, chl-a concentrations were comparable between bays; after, in Fangar Bay, chl-a concentrations decreased more sharply.Despite the differences in chl-a concentration, the general trend was similar in both bays almost all the year, but chl-a peaks differed.In Alfacs Bay, chl-a peaked during February-April and October-November, achieving maximum concentrations in March.In Fangar Bay, chl-a peaked on May and September-October, being the most productive along the year in this last period.Minimum chl-a concentrations were found in winter in both bays, January and November for Alfacs and Fangar Bay, respectively.The coefficient of variation (CV) of chl-a (Figure 10a,b) along the year showed, in general, higher CV in Fangar than Alfacs Bay.In Fangar, higher variability was observed in the mouth of the bay, associated to higher chl-a concentrations, while lower CV values were found at the inner area and within the shellfish rafts, where lowest values of chl-a were found.Conversely, in Alfacs Bay, higher CV was observed in the eastern half of the bay, especially in the inner area and the eastern half of the rafts polygon, with lower averaged chl-a concentration.The harbor area and its surroundings, including the western half of rafts and the mouth of the bay, showed lower values of CV, related to higher concentrations of chl-a.The time-series (Figure S3) was revised according to the four different rice-paddies irrigation network scenarios (i.e., Closed channels in winter, semi-closed channels in spring, opened channels in summer and semi-opened channels in autumn)and aquaculture production.The closure of the discharging channels (closed, semi-closed, and semi-opened) propitiated a more eutrophic environment, reaching higher chl-a concentrations than during the opened channels stage, this phenomenon being more evident in Alfacs Bay.During the closed and semi-closed stage (from January to April), chl-a tended to increase in both bays, but the increment was much more pronounced and long-lasting in Alfacs.During these months rice paddies are dry and so, the supply of water from the channels is minimum.Regarding the chl-a within the shellfish rafts, while Fangar Bay showed similar chl-a concentrations inside and outside the rafts (more homogeneous bay), in Alfacs Bay, lower concentrations of chl-a were observed inside the rafts' polygon.During the opened channels stage (from April to September), chl-a concentration decreased in Alfacs and remained the same in Fangar Bay, compared with the prior period.However, from late July, both bays showed an increase of chl-a concentration that lasted until the end of September, when chl-a dropped sharply, achieving values close to 0 mg/m3 in both bays.The opened channels stage is characterized by high freshwater inputs with the maximum occurring in September-October.Despite shellfish filter more actively during the warm months, no significant differences in chl-a were observed between the rafts  The time-series (Figure S3) was revised according to the four different rice-paddies irrigation network scenarios (i.e., Closed channels in winter, semi-closed channels in spring, opened channels in summer and semi-opened channels in autumn)and aquaculture production.The closure of the discharging channels (closed, semi-closed, and semi-opened) propitiated a more eutrophic environment, reaching higher chl-a concentrations than during the opened channels stage, this phenomenon being more evident in Alfacs Bay.During the closed and semi-closed stage (from January to April), chl-a tended to increase in both bays, but the increment was much more pronounced and long-lasting in Alfacs.During these months rice paddies are dry and so, the supply of water from the channels is minimum.Regarding the chl-a within the shellfish rafts, while Fangar Bay showed similar chl-a concentrations inside and outside the rafts (more homogeneous bay), in Alfacs Bay, lower concentrations of chl-a were observed inside the rafts' polygon.During the opened channels stage (from April to September), chl-a concentration decreased in Alfacs and remained the same in Fangar Bay, compared with the prior period.However, from late July, both bays showed an increase of chl-a concentration that lasted until the end of September, when chl-a dropped sharply, achieving values close to 0 mg/m3 in both bays.The opened channels stage is characterized by high freshwater inputs with the maximum occurring in September-October.Despite shellfish filter more actively during the warm months, no significant differences in chl-a were observed between the rafts The time-series (Figure S1) was revised according to the four different rice-paddies irrigation network scenarios (i.e., Closed channels in winter, semi-closed channels in spring, opened channels in summer and semi-opened channels in autumn)and aquaculture production.The closure of the discharging channels (closed, semi-closed, and semi-opened) propitiated a more eutrophic environment, reaching higher chl-a concentrations than during the opened channels stage, this phenomenon being more evident in Alfacs Bay.During the closed and semi-closed stage (from January to April), chl-a tended to increase in both bays, but the increment was much more pronounced and long-lasting in Alfacs.During these months rice paddies are dry and so, the supply of water from the channels is minimum.Regarding the chl-a within the shellfish rafts, while Fangar Bay showed similar chl-a concentrations inside and outside the rafts (more homogeneous bay), in Alfacs Bay, lower concentrations of chl-a were observed inside the rafts' polygon.During the opened channels stage (from April to September), chl-a concentration decreased in Alfacs and remained the same in Fangar Bay, compared with the prior period.However, from late July, both bays showed an increase of chl-a concentration that lasted until the end of September, when chl-a dropped sharply, achieving values close to 0 mg/m 3 in both bays.The opened channels stage is characterized by high freshwater inputs with the maximum occurring in September-October.Despite shellfish filter more actively during the warm months, no significant differences in chl-a were observed between the rafts and their surroundings in neither of the two bays.Finally, the semi-opened channels stage (from October to December) started with a strong increase of chl-a at both bays in October, which lasted until late November, when chl-a concentrations dropped close to 0 at both bays.During December, Alfacs Bay recovered chl-a concentrations similar to those of the opened channels stage, but Fangar Bay kept low chl-a values.The semi-opened channels stage implies that water is still being discharged in the bay, but contributions decrease with time.Most of the shellfish are harvested during summer, so the bivalve grazing pressure is reduced the last months of the year.Although similar chl-a concentrations were found between the rafts and the rest of the bays, lower values of chl-a tended to aggregate in the middle area of the Alfacs rafts' polygon, and the Northern rafts of Fangar Bay.

In situ chl-a Data
Three different laboratory methods for chl-a quantification from water samples were compared.
Chlorophyll-a concentration measured by spectrophotometer (SP) after acetone extraction was better correlated with satellite data.The in vivo method is only used as a fast qualitative proxy of chl-a due to its sensibility to errors with unknown uncertainty (i.e., overestimation due to non-phytoplanktonic contribution), while extracting the pigment with a solvent (i.e., alcohol-based or acetone) and measuring it with the fluorometer or spectrophotometer is the common procedure in remote sensing of chl-a [21,39,40].Regarding the use of surface or integrated water samples for ground truth chl-a quantification, the vertical distribution of the phytoplankton biomass might have a significant impact on the remote sensing signal.In Fangar Bay, significant differences were not found between surface and integrated water column chl-a concentrations.This finding is in agreement with Ramón et al. [41], who found homogeneous chl-a concentration by depth in a 10 month study (1 sample per month) in Fangar.These results suggest the use of an integrated water column chl-a for remote sensing model calibration and validation in coastal shallow waters, but further research should include data of both bays under different scenarios to prove the validity of this assumption during the year.

Atmospheric Correction and chl-a Estimation Algorithms
ACOLITE was used for atmospheric correction of Sentinel 2A L1C images using the Dark Spectrum Fitting (DSF) based on the SWIR bands.The results showed that TOA contributed over 50% for all MSI over surface reflectance of water pixels (Rhow).This might be related with non-negligible water reflectance in the SWIR band.According to [42], the invalidity of the SWIR black pixel assumption could lead to an overcorrection of the reflectance (SWIR reflectance for water pixels was up to 10 times larger when solar zenith < 42 • ; i.e., spring and summer).In this study, it was not possible to validate the atmospheric correction with field radiometric measurements, but the drop of reflectance of Rhow images compared to Rhot in the blue bands was noticeable.However, a strong reflectance peak was observed in the green part, independent of the level of atmospheric correction applied.Similar results were obtained by [43] using ACOLITE without sun glint correction in an estuarine area, and they also found higher water reflectance in all bands in areas with higher concentration of total suspended matter.In the Ebro Delta, Fangar Bay always showed larger reflectance at all spectrum compared to Alfacs Bay.Fangar is shallower and thus is more susceptible to wind-driven mixing and sediment resuspension.However, the increased reflectance of Fangar Bay might be also related to bottom reflectance or contamination due to adjacency effects.As suggested by other authors [44,45], adjacency effects have significant impact in coastal waters due to typical lower reflectance in relation with their neighbourhood surfaces (i.e., sand beaches and rice paddies), increasing the apparent brightness.This effect might be more pronounced in Fangar Bay due to its geomorphological characteristics (smaller, shallower, and more closed).
The simplified atmospheric correction procedure, which normalizes the TOA signal for Rayleigh effects, was preferred in favor of a full aerosol atmospheric correction given the large uncertainties associated with water leaving reflectance over turbid waters.Ref. [31] found a similar performance between the median spectra derived for full atmospheric correction and only Rayleigh correction.Our models showed better performance using Rhorc images instead of Rhow.For both bays together, the best performing algorithm was the BG ratio.Common ocean color algorithms based on the ratio of blue and green bands do not perform well in optically complex coastal waters (less sensivity to chl-a concentration changes) [20,23,36,46].However, Ref. [46] suggested that, in case 2 oligotrophic waters ([chl-a] < 4 mg/m 3 ), the use of blue and green wavelengths is more appropriate.In the present study, averaged chl-a concentrations (measured in situ) were 1.50 and 2.99 mg/m 3 , with maximums of 2.60 and 5.60 mg/m 3 , for Fangar and Alfacs Bays, respectively, and the results achieved were in agreement with [47,48].The chl-a estimates were reasonably well derived (MAE = 0.63 and APD < 10%) using the BG ratio.

Model Calibration and Validation
In order to reduce noise and minimize temporal-gap effects, ground truth chl-a data were averaged over a 3 × 3-pixel box, centered on the sampling point.Despite the averaged reflectance is commonly used, it might be a poor measure of central tendency if the set of pixels used for its calculation contains outliers.Here, outliers were removed before averaging the 3 × 3-pixel box.Although the number of valid pixels differed among locations and dates, the number of valid pixels (between 2 and 9) was not correlated with the mean standard deviation, thus demonstrating the suitability of the results.After model selection, the BG applied over Rhorc images was preferred and used to make a "pixel-stability mask" to identify and reject those areas where the values obtained with the integration of the remote sensing and the model were not responding to the changes of chl-a.Based on k-means clustering, the applied methodology allowed us to distinguish the boundaries where maximum changes occurred, thus defining the edges for the delimitation of the mask.In Alfacs Bay, better results were achieved that were able to differentiate each raft individually (rafts more separated than in Fangar) and masking the shallow waters (confined only to the margins of the bay).
The applied model was based on algorithms specifically tuned for Alfacs and Fangar bays.Despite the good results achieved for the CalVal dates, the suitability of the model depends on the ratio between the range of remote sensing and the range of the available ground data and its representativeness along different seasons or scenarios.In our study, not many samples were available, but their spatial distribution covered a wide range of in-day scenarios at each bay.However, most of the samples for CalVal purposes included only the seasons of summer and spring so the application of the models over winter and autumn was subjected to higher uncertainty.Indeed, the range of chl-a measured in situ included low number of samples with concentrations under 1 mg/m 3 , which are highly representative in the winter season.This fact coupled with the linearity of the developed models, increase the error related to low chl-a concentrations, tending more rapidly to negative numbers (e.g., Figure S1b).

Spatiotemporal chl-a Dynamics
The time-series of chl-a covered all year 2017, including the different channel stages at both Fangar and Alfacs bays.Overall, the temperature increased from winter (T mean ~14.4 • C) to summer (T mean 27.7 • C), and in autumn, the temperature (T mean ~22.5 • C) was similar, even higher, than in spring (T mean ~20.1 • C).The most frequent winds during the year came from the NW sector, predominantly in the morning, with strong influence of southern winds (spring and summer), switching to SSE (winter and summer) and to SSW (spring and autumn) from noon onward.Within the dates included in this study, highest intensities were registered in March and December, both related to direction of 300-360 • (NW and NNW).The rainiest month along year was January (72.4 mm in seven days), followed by March (36.5 mm in 14 days), and within the days included in the time-series, it rained on 4 August between 4:30 and 5:00 am (accumulated precipitation of 0.1 mm) and on 23 September from 4:30-6:30 am (accumulated precipitation of 0.3 mm).
Regarding the variability of chl-a inside each bay (in terms of CV), Fangar Bay showed higher heterogeneity along time, but more homogeneity along space than Alfacs Bay.Fangar Bay is smaller and shallower, which makes it more susceptible to environmental variations, making the changes faster and affecting more of the bay's area.In this study, it has been observed that in Fangar maximum variations were associated to more energetic areas with more chl-a (mouth), while in Alfacs, higher variabilities were associated to less energetic areas where chl-a dynamics depend largely on the wind-driven mixing (inner NE area).These findings are related with the renewal time of the bays (higher in Alfacs) and linked to the capacity for developing larger phytoplankton populations (higher chl-a concentrations).In this sense, Alfacs Bay characteristics (larger and deeper bay more perpendicular to NW and N winds with higher water residence times) allow nutrients to sink and get stored in the sediment of the bay and, at the same time, allow them to be released and suspended during more time (increased nutrient availability for phytoplankton).Conversely, quicker changes in Fangar Bay make chl-a to be diluted faster by the Mediterranean water inputs (less productive waters).
In relation to chl-a concentrations, besides the seasonal temperature-driven dynamics, wind was the environmental parameter more related to the maximum variations of chl-a inside the bays.In terms of temporal dynamics, overall, chl-a increased more with prolonged NW and N strong winds episodes occurred in 17 March both bays (Figure S1e,f), 25 July both bays (Figure S1u,v), 13 September Fangar Bay (Figure S1ab), 23 October both bays (Figure S1ag,ah), and 12 November Alfacs Bay (Figure S1ai).The highest accumulations of chl-a at both bays occurred on March (Figure S1e,f), October (Figure S1ae-ah) and early November (Figure S1ai), when channels were closed or semi-closed.Conversely, weaker winds from southern components were related to decreases in chl-a concentrations as happened in 15 July both bays (Figure S1s,t), 23 September both bays (Figure S1ac,ad) and 22 November both bays (Figure S1aj,ak).Reduction of chl-a concentration in both bays was enhanced after rainy events as in 23 September (Figure S1ac,ad) and 22 November (Figure S1aj,ak).These results suggest that wind plays a major role in the nutrient load of the water column.On one hand, mixing the water re-suspending the sediment, thus making the nutrients available (wind-induced turbulence).On the other, enhancing water renewal which increases flushing cells ratio and does not allow time enough for development of large phytoplankton populations [49] (wind-enhanced circulation).Therefore, higher chl-a concentrations are expected to occur when the estuarine circulation is weakened and the turbulence increases.This effect might be enhanced at the end or after drainage of the irrigation channels stage (August-November) which increase the nutrients stored in the bays.In the time-series presented in this study, this occurred in August at both bays (Figure S1w-aa) and 23 October both bays (Figure S1ag, ah), especially when winds blow from land (NW and N).In general, the observed trend was in agreement with previous studies [49], which found high chl-a concentrations of chl-a in October in Alfacs Bay and high concentrations from July to November in Fangar Bay.
In terms of spatial dynamics (in-day scenarios), high chl-a concentrations were related more frequently to the mouth area of both bays.There, the exchange with the Mediterranean Sea leads to a more instable scenario in which, despite water renewal might be higher, increased turbulence favors phytoplankton growth prevailing over the wind's regime.High concentrations in the mouth of Fangar Bay were related with more chl-a within the central channel of the bay (northern face of rafts), while in Alfacs Bay there was not so clear relation.In this bay, highest chl-a concentrations were also found in the inner area (NE), which is more retentive and concentrate more nutrients [5].
Because of all the aforementioned, shellfish aquaculture in the bays benefits from increased chl-a concentrations compared to the open sea.However, the retentiveness that characterizes the bays become a double-edged sword due to the high temperatures that water reaches during summer (>30 • C), which negatively affects the feeding rate of shellfish, becoming lethal when temperatures above 28 • C are maintained for days [41].In order to develop a feasible method for aquaculture management by means of remote sensing monitoring, temperature must be included as one of the main factors, together with chl-a, controlling shellfish growth and conditioning the sustainability of the cultures.
In this article, the first results have been presented, and measures to enhance aquaculture can be proposed.However, the feasibility of implementing them is subjected to the availability of bio-geophysical models considering longer time-series, which would allow to make a more integrated and robust approach.Including more data (parameters considered, increased number of data, wider dynamic range) and integrating them into the models would lead to carry out studies in line with [27,30], which coupled remote sensed chl-a with other environmental parameters to establish shellfish farming suitability index, to determinate the load capacity of the bays, and to rezoning the rafts' locations.

Conclusions
Moderate spatial resolution (10-60 m 2 ) Sentinel 2 imagery offers a new opportunity for remote sensing of water quality at small coastal geographic areas.In the Ebro Delta bays, the main Spanish Mediterranean shellfish production site, Sentinel 2 imagery has demonstrated the potential to become a suitable tool for resolving the fast dynamics of phytoplankton in the area (in terms of chl-a concentration), within short space and time-frames.The monitoring using satellite remote sensing improves the standard in situ sampling-based methodology, allowing moving from punctual to full coverage, thus enabling holistic analyses (time-series) to enhace coastal management (e.g., aquaculture).
After testing different levels of atmospheric correction, it is not feasable to use uncorrected atmosphere images (TOA), but the full correction of the atmosphere is still highly uncertain.The results obtained suggest the use of Rayleight corrected Sentinel 2 imagery together with a simple Blue-to-Green ratio for chl-a estimation, until full correction is completely solved/validated.With this configuration, APD < 10% and MAE < 0.7 mg/m 3 were achieved, being able to derive credible chl-a maps of the bays, including the preservation of some information within the rafts polygons.
Despite the aforementioned success, remote sensing of small complex coastal geographic areas still faces several challenges.The main limitations found in this study were (i) full atmospheric correction accounting for aerosol, Rayleigh, sun glint, and adjacency effects and (ii) uncertainties associated to shallower areas contaminated by bottom reflectance, contributions of seagrasses to the total chl-a concentration, and validity of the results out the range of derivation of the model (location of ground truth data, wider range of chl-a concentrations, and seasonality).Further research should be directed to solve these shortcommings by improving the atmospheric correction and gathering more field data covering higher number of scenarios.With these, a sensivity test should be conducted for algorithm bounding, and, ideally, specific tunned models should be developed for each scenario (bay/season/water optical properties).

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/15/1756/s1,Table S1: Model performance (in situ chl-a range and N, spectral algorithms' range, intercept, slope, RMSE, APD, Pearson's r, R 2 , AIC and BIC) with all variable combinations (chl-a method, bay, atmospheric correction level and spectral algorithm) using all available in situ chl-a data.Table S2: Model performance (in situ chl-a range and N, spectral algorithms' range, intercept, slope, RMSE, APD, Pearson's r, R 2 , AIC and BIC) with all variable combinations (chl-a method, bay, atmospheric correction level and spectral algorithm) using only the samples for which chl-a was measured by the three methodologies (in vivo, FL and SP).Funding: This research was funded by Agència de Gestió d'Ajuts Universitaris i de Recerca (grant number 2018 FI_B00705), the European Maritime and Fisheries Fund (ARP005/17/00076) and Fisheries Directorate from Generalitat de Catalunya.

Figure 1 .
Figure 1.Location of the study bays, meteorological station, mussel rafts, coastal lagoons, irrigation fields, and the discharging channels in Ebro Delta.

Figure 1 .
Figure 1.Location of the study bays, meteorological station, mussel rafts, coastal lagoons, irrigation fields, and the discharging channels in Ebro Delta.

Figure 3 .
Figure 3. Temporal distribution of Sentinel 2 images used in this study for calibration and validation and time-series development (TSD) at Alfacs and Fangar bays.

Figure 3 .
Figure 3. Temporal distribution of Sentinel 2 images used in this study for calibration and validation and time-series development (TSD) at Alfacs and Fangar bays.

Figure 4 .
Figure 4. Workflow to derive chl-a time-series from Sentinel 2A multispectral imagery (MSI) data at Ebro Delta bays for aquaculture management purposes.

Figure 4 .
Figure 4. Workflow to derive chl-a time-series from Sentinel 2A multispectral imagery (MSI) data at Ebro Delta bays for aquaculture management purposes.

Figure 5 .
Figure 5. Daily averaged reflectance spectra per bay for each band of Sentinel 1C and 2A products on Calibration and Validation dates.Alfacs Bay: solid line.Fangar Bay: dashed line.

Figure 5 .
Figure 5. Daily averaged reflectance spectra per bay for each band of Sentinel 1C and 2A products on Calibration and Validation dates.Alfacs Bay: solid line.Fangar Bay: dashed line.

Figure 6 .N
Figure 6.Calibration (a) and Validation (b) results of the Blue-to-Green ratio (BG) algorithm, chl-a SP, and both bays together over the set of Rhorc images.The 95% prediction (dashed line) and the 95% confidence interval (dotted line) are also shown.

Figure 6 .
Figure 6.Calibration (a) and Validation (b) results of the Blue-to-Green ratio (BG) algorithm, chl-a SP, and both bays together over the set of Rhorc images.The 95% prediction (dashed line) and the 95% confidence interval (dotted line) are also shown.

Figure 6 .
Figure 6.Calibration (a) and Validation (b) results of the Blue-to-Green ratio (BG) algorithm, chl-a SP, and both bays together over the set of Rhorc images.The 95% prediction (dashed line) and the 95% confidence interval (dotted line) are also shown.

Figure 7 .NFigure 7 .
Figure 7.Linear relationship between all available chl-a SP and BG per bay.The 95% prediction (dashed line) and the 95% confidence interval (dotted line) are also shown.3.4Chlorophyll-a Time-Series

Figure 8 .
Figure 8. Masking Rhorc high reflectance over blue and green bands (threshold = 0.11).Exempla of cloud presence in Alfacs Bay on 22 November.

Figure 8 .
Figure 8. Masking Rhorc high reflectance over blue and green bands (threshold = 0.11).Exempla of cloud presence in Alfacs Bay on 22 November.

Table 1 .
Summary of chl-a samples coinciding with Sentinel 2A pass.Grid 1: routine sampling for the official water-quality monitoring program (see orange dots in Figure

Table 2 .
Descriptive statistics of chl-a concentration (mg/m 3 ) per bay and measuring method, during the study period.FL = Fluorometer; SP = Spectrophotometer.N is the number of samples, SD is the Standard Deviation, and CV is the Coefficient of Variation.

Table 3 .
Summary of in situ chl-a data set used in the calibration and validation process of the Sentinel 2 derived data.

Table 4 .
Summary of the best performing models per bay for the calibration dataset."Algorithm" refers to the spectral algorithm applied to Rhorc images and calibrated with chl-a spectrophotometer (SP).

Table 4 .
Summary of the best performing models per bay for the calibration dataset."Algorithm" refers to the spectral algorithm applied to Rhorc images and calibrated with chl-a spectrophotometer (SP).