Accurate Estimation of Chlorophyll ‐ a Concentration in the Coastal Areas of the Ebro Delta (NW Mediterranean) Using Sentinel ‐ 2 and Its Application in the Selection of Areas for Mussel Aquaculture

: Multispectral satellite remote sensing imagery, together with appropriate modeling, have been proven to provide chlorophyll ‐ a maps that are useful to evaluate the suitability of coastal areas for carrying out shellfish aquaculture. However, current approaches used for chlorophyll ‐ a estima ‐ tion in very shallow coastal areas often fail in their accuracy. To overcome this limitation, an algo ‐ rithm that provides an accurate estimation of chlorophyll ‐ a concentration in the coastal areas of the Ebro delta (North Western Mediterranean) using atmospherically corrected Sentinel 2 (S2) remote sensing reflectances (Rrs) has been calibrated and validated. The derived chlorophyll ‐ a maps cre ‐ ated have been used in a dynamic carrying capacity model that covers areas from very rich waters inside the embayment to the more oligotrophic waters in the open sea. The use of carrying capacity models is recommended to evaluate the potential of marine coastal areas for bivalve mollusk aqua ‐ culture. In this context, the depletion of chlorophyll ‐ a is an indicator of negative environmental im ‐ pact and thus a continuous monitoring of chlorophyll ‐ a is key. The proposed methodology allows estimation of chlorophyll ‐ a concentration from Sentinel ‐ 2 with an accuracy higher than 70% in most cases. The carrying capacity and the suitability of the external areas of the Ebro delta have been determined. The results show that these areas can hold a significant mussel production. The meth ‐ odology presented in this study aims to provide a tool to the shellfish aquaculture industry.


Introduction
Remote sensing enables the evaluation of the spatial and temporal variability of water quality worldwide, overcoming the lack of data from new, remote, or large marine areas. This has led to several applications of remote sensing in bivalve mollusk aquaculture, including spatial planning [1][2][3][4][5], early detection of harmful algal blooms [6][7][8][9][10][11], and detection of microbial contamination [12,13]. Among the different water quality parameters of interest that can be measured by sensors from satellites, chlorophyll-a concentration is commonly used in bivalve mollusk aquaculture because it is considered the best proxy of phytoplankton biomass [14], which is one of the ecosystem component Essential Ocean Variables (EOVs) [15,16]. The retrieval of water quality parameters (such as chlorophyll-a) using satellite remote sensing imagery and appropriate modelling is still challenging in coastal waters, especially in case-2 waters [17,18]. Although many approaches have been proposed, there is not yet any standardized approach. In case-2 waters, inorganic and/or organic sediments represent an important or dominant contribution to the optical properties [19], requiring high accuracy and precision in the atmospheric correction algorithms to successfully retrieve water constituents. Recently, several advances concerning enhanced modelling have been proposed to solve atmospheric correction [20], but their performances differ depending on the scenario (sun and observation geometry, atmospheric, optical) and site-specific conditions. The lack of a common approach and performance makes it necessary to continue validating different atmospheric correction approaches as well as water quality retrieval methods with in situ data accounting for a wide variety of water types and environmental conditions. Therefore, the estimation of chlorophyll-a concentration in coastal waters still has some difficulties in comparison to oceanic waters due to the more complex optical properties and to high spatial variability. Atmospheric correction and the correction of scale effects are necessary to accurately estimate chlorophyll-a concentration in coastal waters [21], where the retrieval of chlorophylla concentration through remote sensing has many applications, such as the evaluation of the environmental status of these water masses [17] or spatial planning for bivalve aquaculture [1].
Global aquaculture production in 2020 reached 17.7 million tons of mollusks, with bivalve aquaculture accounting for a high percentage of total aquaculture production in several countries such as New Zealand, France, Spain, republic of Korea, Italy, and Japan [22], with China being the main producer [23]. Bivalves do not require artificial feeds; they graze on phytoplankton communities present in marine waters [24]. Bivalve aquaculture, together with the cultivation of macroalgae and some herbivorous fish species, utilizes low trophic marine species that can provide alternatives for food production with lower environmental impacts [25]. Bivalve aquaculture is carried out in coastal waters, mainly in estuaries and embayments where phytoplankton biomass is higher than in offshore waters. The global production of bivalves is expected to continue growing [26]; however, seawater temperature in shallow coastal areas is currently reaching thresholds that are limiting the survival of some bivalve species.
Heat waves have been affecting, in recent years, shellfish species in different geographical areas such as the Pacific Northwest in 2021 [27,28], the eastern English channel in 2018 [29], and the coast of California in 2004 [30]. In the Mediterranean Sea, the warming trend in Sea Surface Temperature (SST) during the period 1992-2020 was 0.038+/−0.002 °C/year [31]. Marine heat wave (MHW) events in the Mediterranean Sea are increasing in intensity, duration, and frequency [32]. High seawater temperature in shallow coastal areas during MHW events produces high rates of mortality of Mytilus galloprovincialis [33], one of the main species in shellfish aquaculture in the Mediterranean Sea. These issues revealed the need for seeking new shellfish growing areas located in cooler waters, which requires the assessment of the production potential of new culture areas, with a sustainable approach. In this context, the coupling of remote sensing data and carrying capacity models has demonstrated a great potential for monitoring shellfish aquaculture and supporting management strategies [34,35].
Carrying capacity models are used to determine the total bivalve biomass supported by a given ecosystem as a function of the water residence time, primary production, and bivalve clearance rate [36,37]. Among them, the Dynamic Energy Budget (DEB) theory based on Kooijman [38] has been applied to model shellfish growth [39][40][41][42] in order to evaluate food limitation for different bivalve species such as Macoma balthica, Mya arenaria, Cerastoderma edule, Mytilus edulis, and Crassostrea gigas [43], and for Mytilus galloprovincialis [44], and to predict the impact of environmental changes on shellfish populations [45]. Used alone, or coupled with other biogeochemical [46], ecophysiological, and box models [47], DEB models have been proven to be a valid methodology to determine the quality of an area for shellfish aquaculture [34,35]. DEB models allow one to calculate the amount of energy used for shellfish growth and reproduction through the ingestion of the available food, thus they benefit from remote sensing estimates [48][49][50][51][52].
The Ebro delta, located in the western Mediterranean Sea, is a natural protected area where agriculture, fisheries, aquaculture, and tourism are the main economic activities. The aquaculture industry in the Ebro delta produces 3400 Tn of mussels and 300 Tn of oysters per year (average 2010-2021 [53]). Shellfish aquaculture is carried out in the two coastal embayments located in the northern (Fangar Bay) and the southern (Alfacs Bay) shores. The hydrodynamics in these embayments is highly sensitive to wind variability, which affects circulation patterns, water exchange times [54], and phytoplankton dynamics [55]. Drainage channels discharge freshwater from the rice fields into the embayments; the flow of freshwater is variable during the rice growing season from May to September, and it is almost negligible during the rest of the year. Shallowness and large water renewal times, coupled with the exacerbation of MHW events, is limiting the development of aquaculture in these embayments. Previous studies have focused on estimating chlorophylla concentration inside the embayments from remote sensing data (i.e., Sentinel-2 Multi-Spectral imagery), either for monitoring phytoplankton biomass [56], or to evaluate the impact of coastal storms on water quality [57]. None of them have attempted to model the carrying capacity of the embayments or to evaluate the suitability of the nearby Mediterranean waters for shellfish aquaculture. These are crucial topics for the future of shellfish aquaculture in the region of the Ebro delta and other Mediterranean areas.
This study focuses on the development of a methodology to create, calibrate, and validate an algorithm that provides an accurate estimation of chlorophyll-a concentration in the coastal areas of the Ebro delta (NW Mediterranean) using atmospherically corrected Sentinel 2 (S2) remote sensing reflectances (Rrs) and a tool to calculate the carrying capacity of the potential new aquaculture areas and to evaluate its suitability. For this purpose, the DEB model has been implemented using two-year estimates of chlorophyll-a concentration derived from Sentinel-2 imagery and additional site-specific auxiliary data. The proposed approach is based on open-source data and open-source software tools. The use of remote sensing data will help the low trophic aquaculture industry to expand their activities to new suitable areas.

Study Area
The Ebro delta extends 25 km from the mainland and forms two embayments, Alfacs and Fangar ( Figure 1). Fangar Bay, located in the north, has an extension of 12 km 2 and a mean depth of 2 m, it contains 16 × 10 6 m 3 of seawater, and the mouth of the Bay is 1 km wide. Alfacs Bay is in the south, its extension is 50 km 2 , the mean depth is 4 m, it contains 200 × 10 6 m 3 of seawater, and its mouth is 3 km wide. Seawater temperature inside the embayments presents wide fluctuations in comparison to open waters. Both bays receive freshwater inputs from the agriculture drainage channels, favoring stratification during calm conditions [58]. The average salinity is similar in both bays: 34.6 in Fangar Bay and 34.8 in Alfacs. Freshwater inputs are 228 × 10 6 m 3 /year and 365 × 10 6 m 3 /year in Fangar and Alfacs, respectively. Water renewal time is approximately 1-2 days in Fangar and 10 days in Alfacs Bay [58]. Chlorophyll-a concentration increases during wind events due to horizontal mixing and bottom resuspension [55]. Both embayments currently hold mussel and oyster production; the location of mussels and oyster rafts is shown in Figure 2.

In Situ Measurements
Sampling cruises were conducted, coinciding with a Sentinel-2 (S2) pass when the weather forecast was favorable. A total of 17 sampling cruises were carried out during the period 2 September 2020 to 27 October 2021, 9 of them covering the northern area and 8 for the southern area of the Ebro delta (western Mediterranean). Ten sampling stations were visited during each cruise, including waters inside and outside the bays. The location of the 20 sampling stations is shown in Figure 1. At each station, bottom depth, and the Secchi disk depth (ZSD) were measured; the profiles of temperature and salinity were obtained using a SeaBird19plus CTD and water samples were taken for chlorophyll-a analysis. Upon arrival at the laboratory, water samples were filtered through Whatman ® glass microfiber filters, Grade GF/F (0.7 μm) 47 mm, using low vacuum. The filters were maintained at -80°C until analysis. Filter contents were extracted in acetone 90%, the absorbance of the extract was measured in a Shimadzu UV-1800 UV/Visible Scanning Spectrophotometer, and the chlorophyll-a concentration was calculated using the formula Chlorophyll-a = 11.85 E664 − 1.54 E647 − 0.08 E630 [59].

Maps of Chlorophyll-a from Sentinel-2
The Sentinel-2 (S2) constellation consists of two satellites (S2A and S2B) operated by The European Space Agency (ESA). Each satellite has on-board the MultiSpectral Instrument (S2-MSI). The S2-MSI Level-1 (L1C, Top of atmosphere) imagery includes 13 spectral bands centered at different wavelengths from 443 nm to 2200 nm with different spatial resolutions (10, 20, and 60 m). All available S2-L1C images between 1 October 2019 and 30 September 2021 were downloaded from the Copernicus Services Data Hub (https://cophub.copernicus.eu/dhus/#/home, accessed on 14 July 2022). The only orbit completely covering the study area was orbit #51. Every day the data from this orbit were available, the four tiles needed to cover the area were downloaded. This amounted to a total of 560 images (140 images per tile), covering a surface of 100 × 100 km 2 for each tile.
Each image was visually checked for clouds and shadows over the region of study, leading to the rejection of 350 images (62.5%), implying the rejection of 46% of dates with chlorophyll-a in situ data. Due to project administrative deadlines, only images until September 2021 were used to accomplish the study.
The pre-processing of S2 imagery was performed through the Graphic Processing Tool (GTP) of SNAP v8.0 [60] and Rcore 3.6 [61]. For atmospheric correction of S2-MSI L1C imagery, the Case 2 Regional Coast Colour (C2RCC) and the C2X-COMPLEX (C2XC) processors (C2-Nets) [62] included in SNAP were applied on all valid S2 images. These processors are based on a multi-sensor per-pixel artificial neural network (NN) method, and differ in their training ranges of inherent optical properties [63]. The parametrization for the atmospheric correction of each image included: (a) pressure (hPa) from NCEP/DOE Reanalysis II data provided by the NOAA PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov (accessed on 16 July 2022) [64]; (b) atmospheric ozone in Dobson units (DU) from the Aura OMI NASA dataset [65], downloaded for each location and date from https://oceancolor.gsfc.nasa.gov/ (accessed on 17 July 2022) [66]. Surface seawater temperature and salinity were obtained from the sampling cruises. For land/water segmentation, the valid pixel expression was set as a threshold on the shortwave infrared (SWIR) band of S2-MSI L1C images centered at 1600 nm (B11). The threshold was defined independently for each image with the triangle threshold method. From each C2-Net, remote sensing reflectance (Rrs) of visible and near infrared bands and the pigment absorption product (apig) were generated. C2-Nets flags, which include codes for quality control of pixels, were also exported. For C2RCC and C2XC, independently, images from the same date were merged when more than one tile were available.
A procedure to select match-ups between chlorophyll-a concentration measured from water samples and S2 Rrs was performed for the two C2-Nets. This procedure involves several steps. Firstly a 3×3 pixel window, centered at the coordinates of in situ measurements, was extracted for each date and sampling location, and C2-nets were quality-checked in all extracted pixels by applying the recommended flags [67]. Then, flagged pixels as well as pixels with negative Rrs at bands B1, B2, B3, and B4 were removed from the analysis, as recommended by Cui et al. [68]. The number of remaining pixels within each pixel window was checked; windows with less than 5 remaining pixels were removed from the analysis. Outliers were defined through Boxplot analysis applied to each available pixel window and spectral band (B1-B7 and B8A). Finally, the remaining pixel windows with less than 5 valid pixels were removed from the analysis.
A set of band combinations in the form of spectral indices were computed for all valid Rrs match-ups (Table 1). These spectral indices include visible and red edge S2 spectral bands (B1 to B6; ~443 nm to ~740 nm), exploiting specific chlorophyll-a absorption peaks in the blue and red spectral regions. The red edge bands are used for mitigating the effect of absorption by non-algal particles, yellow substances, and backscattering. For each C2-Net, chlorophyll-a was modelled with all the computed spectral indices and the apig band. Seventy percent of the data were used for model calibration (cal) and 30% for model validation (val). The cal/val datasets were generated randomly. Models were developed for raw and log-transformed Rrs/apig and chlorophyll-a data. Linear, linear piecewise (1 breakpoint), polynomial (2nd, 3 rd , and 4th order), logarithmic, power, and exponential models were tested. The entire process was iterated 100 times with varying cal/val datasets. For each model (C2-Net x Band Combination x Type of fitting), performance was evaluated by means of the following statistics: mean average error (MAE), root mean squared error (RMSE), average percentage difference (APD), BIAS, and Pearson's r. GIT (3 band model) [71] (1/Rrs(B4) − 1/Rrs(B5)) × Rrs(B6) Normalized Difference Chlorophyll Index (NDCI) [72] (Rrs(B6) − Rrs(B5))/(Rrs(B6) + Rrs(B5)) apig (apig) [73] apig To generate the maps of chlorophyll-a, all flagged pixels (atmospheric correction flagging) were removed, as were pixels with negative Rrs in any of the spectral bands of each valid S2 image processed with the C2-Nets. The best performant model was selected according to lowest MAE, APD, RMSE, and BIAS, in this order, and it was applied to all valid images and pixels. Pixels with negative chlorophyll-a values were set to 0 mg/m 3 and pixels with unreliably high chlorophyll-a concentration (>30 mg/m 3 ) were removed. Pixels corresponding to the mussel farm structures or rafts were masked out using an available shapefile to mitigate/avoid pixel mixing problems.

Carrying Capacity Model
For the analysis and extraction of data from the generated images, QGIS 3.24 was used [74]. The area of interest corresponding to the masked images was divided into 2 different regions inside and outside the two coastal embayments (Figure 2). The studied area comprises 37.5 km 2 inside and 37.8 km 2 outside Alfacs Bay (southern area) and 2.7 km 2 inside and 19.7 km 2 outside Fangar Bay (northern area). The external area was divided into 12 polygons, 6 outside each embayment ( Figure 2). The area inside each polygon was 5.5, 7.5, 7, 6.8, 4.8, 6.1 km 2 for A1-A6 and 1.8, 2.5, 2.6, 2.9, 3.8, 6.1 km 2 for F1-F6. Monthly chlorophyll-a concentration was averaged for each area of study, combining the available images for the period. The external areas were divided into different parts (polygons) to evaluate the gradient in chlorophyll-a concentration at different distances from the mouth of the bay. Each polygon was divided into three parts to evaluate the impact of the distance to the shore on the chlorophyll-a concentration.
The chlorophyll-a concentration, retrieved from the S2-L1C images treated with the algorithm calibrated and validated in this study, was used to determine the suitability of the different areas and their carrying capacity. To determine the suitability of each polygon for mussel aquaculture, the DEB model was applied. The model allows us to calculate the amount of energy used for growth and reproduction through the ingestion of the available food using the version from Rosland et al. [75] and the notation from Kooijman [76]. Mussel ingestion rate is proportional to the surface of the mussel, expressed as structural volume. (1) where is the ingestion rate; the values used for each parameter were obtained from van der Veer et al. [77] and are shown in Table 2. Clearance rates ( ) were obtained from Galimany et al. [78]. The measured in situ seawater temperature (°C) at 1 m depth was converted to Kelvin degrees (K). In Equation (2), f is the Michaelis-Menten Equation to scale the ingestion rate to the food concentrations (X); this term scales the amount of food ingested as a function of the food (chlorophyll-a) available. In Equation (4), L is mussel length (cm), where δM is the dimensionless shape coefficient. Td in Equation (5) is the Arrhenius temperature function. The model (Equations (1)-(5)) was computed in Python.
In a previous study [79], the DEB model was found to better reflect field observations of mussel ingestion in Alfacs Bay, in comparison with the Scope for Growth (SFG). Carrying capacity was calculated by applying the recommendations of the Aquaculture Stewardship Council [80]. This standard provides indicators for a first approach to characterize the potential of an area. In the present study, the carrying capacity was defined as the number of rafts/ropes/mussels that can be grown in each area and was calculated as a function of the food availability measured through the chlorophyll-a concentration and the renewal time of the water masses. Chlorophyll-a concentration was obtained from the S2 images treated with the algorithm in this study, and the renewal time was calculated using the current speed at the different external areas.
In Formula (6), ℎ is the renewal time of the water in each area, ℎ is the total amount of chlorophyll-a in the whole volume of each area, and ℎ is the chlorophyll-a concentration of each area obtained from the Satellite images.
Different values of clearance rate, , , were applied by selecting the corresponding values for each time of the year. The percentage of chlorophyll-a available for the cultivated mussels was assumed to be 0.75, and 0.25 was the remnant percentage. A standardized length of 3 m per rope that holds 1250 mussels was also used. These are the usual characteristics of the mussel ropes used in the Ebro delta. For the external areas, the average renewal times were calculated from the Copernicus Marine Service [81] using data from the years 2019 to 2021. For all areas, the volume was calculated for 3 m of depth as standard length of mussel ropes. The renewal times used for the areas inside the embayments were obtained from the bibliography [82]. We have not included primary production time because renewal times are short in the external areas.
The following diagram (Figure 3) shows the relation between in situ data, remote sensing data, and the different models used in this study.

In Situ Measurements
The results from the in situ measurements in October 2021 were not used for the calibration and validation of the satellite images but are included here to provide a wider view of the in situ environmental conditions of the area. The average depth of the 16 sampling stations was 7.7 ± 0.2 m (Mean ± Std. Error) and the range was 3.7-14.8 m. The average Secchi depth was 3.5 ± 0.2 m (Mean ± Std. Error) and the range was 1.3-10.4 m ( Figure  S1). The average chlorophyll-a concentration in the water samples was 2.4 ± 0.1 mg/m 3 (Mean ± Std. Error) and the range was 0.1-8.9 mg/m 3 (Figure 4). There was a statistically significant difference in chlorophyll-a concentration and Secchi disk depth (p ≤ 0.001) between the interior and the exterior sampling stations from the whole studied area (Kruskal-Wallis One Way Analysis of Variance on Ranks). In the southern area, there was no statistically significant difference in the mean values of chlorophyll-a concentration among the different sampling stations (p = 0.118) and dates (p = 0.251) (Two Way Analysis of Variance). In the northern area, there was a statistically significant difference (p ≤ 0.001) between sampling stations and dates. The range in seawater temperature was 10.9-28.2 °C in the northern area outside Fangar Bay and 11.7-29.2 °C in the southern area outside Alfacs Bay. Minimum values were registered in January 2020 and maximum values in August 2020 ( Figure 5). Higher temperature is measured in the southern area in comparison to the northern area.

Maps of Chlorophyll-a from Sentinel-2.
Models performed better with C2XC than with C2RCC, particularly coupled with band combinations including spectral bands in the blue or green, and red or red edge regions. The best performance models for each C2-Net and band combinations are shown in Table 3. The best model for chlorophyll-a estimation was found with the C2XC processor, applying the REB1 band ratio (see Table 3) with a 2nd polynomial fitting ([chlorophyll-a] = −0.615 + 10.88 × REB1 + 8.704 × REB1 2 ). The calibration and validation plots, as well as the fit between observed and estimated chlorophyll-a concentration values, are presented in Figure 6. Maps of chlorophyll-a concentration from 56 different dates (Figure 7) along the period of study were generated using the methodology described in Section 2.3. A good agreement between estimated and observed chlorophyll-a concentration was observed within the sampling dates and the dynamic range of in situ data. However, unreliable high chlorophyll-a concentration was retrieved in some scenarios, such as the image of the 5 February 2020 (Figure 7), captured ~2 weeks after an extreme coastal storm. In those scenarios, dark brownish waters were visible in the RGB image composite, suggesting greater total suspended matter concentration and turbidity, leading to more optically complex waters.

Carrying Capacity Model
Monthly averages of chlorophyll-a concentration for the areas AI, A1-A6 together and FI, F1-F6 together are shown in Figure 8. Higher values of chlorophyll-a concentration were observed inside the two coastal embayments in comparison to the external area. The highest values occurred in Fangar Bay during autumn (13.6 mg/m 3 ) and summer (8.5 mg/m 3 ). In the north (F1-F6), the maximum averaged chlorophyll-a value per area and date (13.8 mg/m 3 ) was detected in November 2019 in F1; in the south (A1-A6) it was detected in A2 in December 2019 (12.1 mg/m 3 ). There are statistical differences between the different areas (Chi-square = 122, DF = 3, p < 0.001). The differences are between FI vs. AE, FI vs. FE, FI vs. AI, AI vs. AE, and AI vs. FE (Tukey test p < 0.05), but not between FE vs. AE (Tukey test p > 0.05). The averaged concentration of chlorophyll-a per date and polygon outside the embayments is shown in Figure S2. The minimum of the averaged values per image, external area, and date (0.12 mg/m 3 ) was measured in July 2021 in A6 (southern area). In the same image, clouds where present in some of the areas, while in other areas of the same image the chlorophyll-a concentration could be measured. There are statistical differences between the different polygons (Tukey test p < 0.05); the results of the statistical analysis are shown in Figure S3. The gradient in chlorophyll-a concentration in the external areas at different distances from the bay's mouth can be observed in Figure S4.
In the northern area outside Fangar Bay, apart from the gradient observed in chlorophyll-a at different distances from the mouth, a second gradient is observed inside each of the external areas at different distances from the shore, the closest areas to the shore being the richest in chlorophyll-a ( Figure S3). A gradient is also observed in the southern areas, but in the opposite way. In polygons A1-A4, the areas closest to the shore contain less chlorophyll-a than those farther from the shore, except in polygons A5-A6, where the same pattern observed in the northern area occurs.
Ingestion rates ( ) from October 2019 to September 2021 for the different polygons of the southern area and from November 2019 to September 2021 for the northern area are shown in Figure 9. In the external southern area, the mean value of was 8.7 ± 0.3 Jd −1 and the range was 1.6-14.  The carrying capacity (CC) calculated for each area is shown in Table 4 as the number of mussel rafts that can be held; a standardized raft in this study contains 1100 ropes of mussels that are 3 m in length. The renewal time for each area and the clearance rate time (h) calculated for the number of mussels and the chlorophyll-a stock of the area are also shown in Table 4. Renewal times are shorter than the clearance rate calculated because 0.75 of the chlorophyll-a stock was selected as available for CC. The renewal time used in Table 4 is the averaged renewal time for the year in each area. Because the renewal time is variable over the year, the CC was also modeled for different renewal times (Figure 10). We can observe that CC tends to stabilize at higher renewal times, depending on the internal areas below 50 rafts and the external areas below 200 rafts. The shaded red area in Figure 10 corresponds to the renewal time range described in previous studies and the red line is the CC calculated for this renewal time. In the external areas we observed the same exponential decreasing result.

Discussion
The values of in situ Secchi disk depth and chlorophyll-a concentration measured in this study (2019-2021) differed slightly from those found in a previous study [82] conducted during the period 1982-1983. The Secchi disk depth (ZSD), at that time, was measured at some of the same sampling stations used for our study. The values obtained in [82] were 0.7-1.9 m deeper than our measurements when comparing the average values for each station. The same authors measured chlorophyll-a concentration and found yearly averages for the inner areas of 3.2 mg/m 3 for AI and 3.44 mg/m 3 for FI; these values are slightly lower than those found in our study: 5.2 mg/m 3 for AI and 4 mg/m 3 for FI. Therefore, the differences in Secchi disk depth can be partially explained by differences in chlorophyll-a concentration. Secchi disk depth (ZSD) is influenced by the phytoplankton biomass and other optical variables such as the colored dissolved organic matter (CDOM) and suspended particulate matter (SPM). Resuspension of sediments and the brownification of coastal waters influence ZSD, and therefore its relationship with chlorophyll-a concentration is not a straight inverse relation. The same authors detected maximum chlorophyll-a values in June-July (3.5-11 mg/m 3 ) and September-October (25 mg/m 3 ), and values in the range of 1-3.7 mg/m 3 for the rest of the year. The maximum chlorophyll-a concentration in water samples measured in our study was 8.9 mg/m 3 in December 2020 (AI). Higher values were observed from the satellite images. It is important to highlight that the values of chlorophyll-a retrieved from the Sentinel-2 images for the different polygons located at different distances from the mouth of the bays reflect the gradients observed in the in situ measurements of chlorophyll-a concentration.
Different methods exist for atmospheric correction in coastal waters, and some of them (C2-Nets, iCOR) are available in open-source tools such as SNAP. This feature opens the door to non-specialized users to generate atmospherically corrected S2 images by using this tool and following the proposed methodology. The model proposed here for mapping chlorophyll-a concentration, and the pixels quality control procedure, can also be implemented using the band-math's functions provided by this tool.
The proposed methodology allows one to estimate chlorophyll-a concentration from Sentinel-2 with an accuracy higher than 70% in most cases (Table 3 and Figure 6). The best performant model was achieved with C2XC, which showed more consistent and accurate Rrs estimates than C2RCC in prior research [20]. However, empirical algorithms, such as the one used, can be expected to perform well only inside their range and for the area they were derived for. They are also limited in their ability to discriminate between non-unique signals from parameters that may be covariant, for example TSM and chlorophyll-a concentration [83]. The developed model is less reliable in conditions of high concentration of TSM and/or CDOM, which may occur after strong winds or storms (increased water turbulence, sediment resuspension, land runoff), as observed on 5 February 2020 ( Figure  7). In these cases, the selected model tends to overestimate chlorophyll-a concentration because the red edge reflectance increases at high TSM values overlapping the absorption of chlorophyll-a in the blue region of the spectrum, which is also affected by greater CDOM concentration (CDOM strongly absorbs light up to 500 nm) [73,84]. In these cases, other empirical or semi-analytical algorithms may be more accurate. Further research should focus on the development of a multi-algorithm blending approach, which has proven to be more suitable across different types of water optical properties [85], albeit it will involve a more complex procedure. The development of an algorithm switchingbased method requires the accurate definition of the optical water type pixel by pixel, accounting for the spectral shape, magnitude, and distinctive Rrs spectral features [86,87] and considering the uncertainty of the Rrs retrieval from C2-Nets across different optical water types [20]. In addition, in situ data covering all different scenarios are desirable, which may be difficult to reach due to complex logistics and meteorological limitations (i.e., cloud coverage), as occurred in our study.
The gradient observed in chlorophyll-a concentration from the inner areas towards the open sea is also observed in the values of ingestion rate ( ). The values of integrate the effect of chlorophyll-a concentration and seawater temperature on mussel physiology. Lower seawater temperature in summer in open waters in comparison to the inner areas could compensate the effect of having less phytoplankton available in these waters, obtaining similar or even better growth when high seawater temperature decreases mussel filtration rates. When we observe the mean values for the period of study in each area, in the north, we can see that the gradient decreases from 9.1 Jd −1 to 6.7 Jd −1 from area F1 to area F5, and it increases slightly to 6.9 Jd −1 in area F6. The gradient is less stepped in the southern area, from 9.3 Jd −1 in A1 to 8.3 Jd −1 in A6. We can employ the value of to rank the suitability of the different polygons for mussel aquaculture; the most suitable would be those with the highest , which are A1 and F1 followed by A2-5. The less suitable would be F3-6. We have based our model on the chlorophyll-a concentration measured in each area. The DEB model reflects well the decrease in food ingestion during summer months due to limiting seawater temperatures [88]. The use of chlorophyll-a as a proxy does not consider changes in food quality due, for example, to the presence of certain phytoplankton species with a different nutritional value. Some authors have detected gaps in DEB model parametrization and proposed to improve seston characterization [89,90], which should be tackled in further research.
A shellfish farm may exceed the ecological carrying capacity when the removal of phytoplankton biomass exceeds the renewal, resulting in a phytoplankton depleted water mass. In the proposed CC model, the renewal time of each area is shorter than the clearance rate time, thus complying with the ASC bivalve standard [80] that recommends comparing how long it takes to clear the body (CT) of water with the renewal time (RT), the ratio CT/RT should be >1. The same standard recommends that in locations where CT < RT, then the ratio CT/PPT should be >3. PPT is the number of days required for the replacement of phytoplankton biomass (PB) in the water body considering phytoplankton growth (PB/PPP; PPP = phytoplankton primary production). When the clearance rate is higher than the renewal rate plus the primary production, the body of water is depleted of phytoplankton and the bivalve stock will have less food available. Other authors suggest employing the regulation ratio defined as a fraction of the phytoplankton turnover rate RR = (1/CT)/(1/PT), where PT is the time it takes to renew the phytoplankton stock in an area [36,91]. We have applied the same formula for CC in the internal and external areas without considering primary production time. Renewal time in the internal areas is longer and it would be appropriate to include primary production time in the calculations. There are already mussel and oyster farms inside the embayments, therefore the chlorophyll-a concentration measured is affected already by the consumption of the bivalves from these farms. More complex models are needed to assess, with more detail, the carrying capacity inside the embayments. Renewal times are not uniform inside the embayments where hydrodynamic models are coupled with biological models such as the Regional Ocean Model System (ROMS) coupled to a biogeochemical nutrient-phytoplankton-zooplankton-detritus model (NPZD) such as in Dabrowski et al. [92]. Hydrodynamical modeling coupled to ecological modeling has been applied in other geographical areas to decide the most appropriate locations of the bivalve farms inside an embayment to avoid negative consequences for the ecosystem [93]. We have shown how small changes in renewal time affect the carrying capacity of the different water bodies in our study. The standard from the Aquaculture Stewardship Council [80] recommends to used retention time calculated as the number of days for tides to flush a volume of water equal to the volume of the area. In the Mediterranean Sea, tides have low amplitude, and so flushing time has to be calculated by other means.
The results on the modeled number of rafts in each area show that these external areas can hold a significant production of mussels in waters that are less impacted by warm events in summer. Mussel growth in these areas may be slower in comparison to the embayments due to lower chlorophyll-a concentration. The differences in seawater temperature measured during this study in summer were 1-2 °C for the same day and depth between internal and external areas; this small difference may be crucial to reduce mussel mortality in summer.

Conclusions
The proposed methodology can be reproduced using open-source tools such as SNAP software, enabling end-users to obtain their own maps of chlorophyll-a concentration from Sentinel-2 images. This methodology allows one to estimate chlorophyll-a concentration from Sentinel-2 with an accuracy higher than 70% in most cases. The best performant model was achieved with C2XC for atmospheric correction of S2-MSI L1C imagery. We have shown that the application of the DEB theory using multispectral remote sensing imagery allows one to rank the suitability of the different areas for shellfish aquaculture. It also allows one to determine the periods of the year with favorable and unfavorable conditions for shellfish growth. The results presented in this study show that the carrying capacity of the different areas is highly variable depending on the renewal times of the water, which at the same time are highly variable during the year. Nevertheless, the external areas located close to the Ebro delta embayments can hold a significant mussel production that would allow the reduction of mussel mortality during summer months due to the lower seawater temperature in these areas. This study shows the potential of the combined use of satellite remote sensing, in situ sampling, and carrying capacity models as a tool for helping the low trophic aquaculture industry to expand their activities to new suitable areas.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs14205235/s1, Figure S1. Secchi depth (m) at the different sampling stations. Note that the axis of sampling stations is reversed in the graphs showing chlorophyll-a concentration from Figure 3. Figure S2. Chlorophyll-a concentration, from Sentinel-2 images, averaged for each polygon and date outside the embayments. Figure S3. Matrix showing the results of the statistical analysis used to detect differences on the chlorophyll-a concentration, from Sentinel-2 images, averaged for each polygon and date outside the embayments. The squares colored in green correspond to polygons where statistical significant differences have been found, red color is for polygons where no differences have been detected, orange is for those polygons where no significant differences were found although they may exist. Figure S4. Chlorophyll-a concentration, from Sentinel-2 images, averaged for each of the 3 parts of each polygon outside the embayments, A-C are areas at different distances from the shore inside each polygon.