Deriving Water Quality Parameters Using Sentinel-2 Imagery: A Case Study in the Sado Estuary, Portugal

Abstract

Monitoring water quality parameters and their ecological effects in transitional waters is usually performed through in situ sampling programs. These are expensive and time-consuming, and often do not represent the total area of interest. Remote sensing techniques offer enormous advantages by providing cost-effective systematic observations of a large water system. This study evaluates the potential of water quality monitoring using Sentinel-2 observations for the period 2018–2020 for the Sado estuary (Portugal), through an algorithm intercomparison exercise and time-series analysis of different water quality parameters (i.e., colored dissolved organic matter (CDOM), chlorophyll-a (Chl-a), suspended particulate matter (SPM), and turbidity). Results suggest that Sentinel-2 is useful for monitoring these parameters in a highly dynamic system, however, with challenges in retrieving accurate data for some of the variables, such as Chl-a. Spatio-temporal variability results were consistent with historical data, presenting the highest values of CDOM, Chl-a, SPM and turbidity during Spring and Summer. This work is the first study providing annual and seasonal coverage with high spatial resolution (10 m) for the Sado estuary, being a key contribution for the definition of effective monitoring programs. Moreover, the potential of remote sensing methodologies for continuous water quality monitoring in transitional systems under the scope of the European Water Framework Directive is briefly discussed.

1. Introduction

Estuaries, defined as water bodies partially saline due to their proximity to coastal waters, but substantially influenced by the freshwater flows [1], are extremely important ecosystems for different life forms. These complex systems are characterized by high levels of primary and secondary productivity, supporting high abundance of organisms and providing key nursery services.

For centuries, estuaries have been considered as regions of extremely high relevance for society as many of the resources encountered on these ecosystems reach high economic and social values. Although human civilizations have historically benefited from estuaries (a majority of the world’s largest cities are located in proximity to estuaries [2]), only in the last decades we have become aware of the disturbing effects of human activities in these habitats [3]. Given the increasing anthropogenic pressure on such important ecosystems, the assessment of water quality is now essential to ensure the best possible water conditions. Therefore, the development of cost-effective tools for these assessments is crucial to promote continuous evaluations and maintain sustainable exploitation of these aquatic systems, keeping their natural balance. As such, it is critical to standardize water quality analysis and turn these processes accessible and understandable, not only for researchers, but also for environmental managers and policy makers, which are the ones who can put the results of this type of analysis into practice [4]. In this context, the European Union created a legislation aiming to balance the requirements of human needs and ecosystems conservation by a sustainable planning of human activities: The Water Framework Directive (WFD) (2000/60/EC). The WFD obliges all Member States to implement water management mechanisms, as well as to assess the ecological status of their water bodies by monitoring and classifying them according to regulated parameters [5]. Routine monitoring of the water quality status is an important step in order to comply with WFD requirements and objectives.

Water quality (WQ) can be measured in terms of biological, chemical, and physical indicators, such as turbidity, chlorophyll-a, harmful algae, temperature, metals, dissolved oxygen, nutrients (primarily phosphorus and nitrogen), and many other contaminants [2]. Here, we focus on core WQ indicators that can be derived directly from ocean colour satellite remote sensing: coloured dissolved organic matter absorption (a_CDOM), chlorophyll-a (Chl-a), suspended particulate matter (SPM) and turbidity. These four indicators are traditionally monitored through field sampling programs, which can be labour intensive, time consuming and expensive. Hence, it is not feasible for some of the member states to carry out continuous monitoring, an issue that needs to be overcome. Moreover, it is often necessary to assume that field samples, which are often limited both spatially and temporally, are representative of the total area of interest. In this context, remote sensing (RS) is a valuable tool for monitoring WQ parameters. The incrementing interest in understanding the potential of this technique is driven by the reduced costs and the high spatial and temporal resolution that allows obtaining information for large areas. Furthermore, the possibility to access historical data allows us to track changes and patterns of the different WQ parameters through time [6].

Nevertheless, RS has been far less successful in transitional waters than in other areas given the high optical complexity of the waters and the proximity to the land (adjacency effect). Retrieving reliable information from RS in transitional systems can be considered a complex task, with a relevant number of challenges involved. The first challenge is performing a good atmospheric correction (AC), a crucial step when using space-born imagery, as it removes the signal coming from the atmosphere from the total signal received by the satellite-sensor system, allowing to isolate the signal coming from just below the water surface. The AC in such mutable and complex environments requires different approaches from the ones used for land and open ocean applications. Moses et al. [7] indicated three main issues that make the AC challenging in optically complex waters: (i) the proximity to terrestrial sources of atmospheric pollution, that may result in optically heterogeneous atmosphere which is difficult to model; (ii) the adjacency effect, which is the contamination of the signal received at the sensor by the contribution of the adjacent land pixels; and (iii) the non-negligible reflectance of water in the near-infrared (NIR) region of the spectra, due to high sediment concentrations preventing the use of the AC schemes adopted for open oceans.

Once the water-leaving signal is retrieved, the complexity of the waters poses another challenge: to retrieve the different parameters contributing to the optical signal. A high range of optical variability can be found among estuaries. These waters can be a mixture of optically shallow and optically deep waters, with different optically active compounds (i.e., phytoplankton, CDOM, SPM) that do not co-vary and/or interact between them. Such aquatic systems can present short spatial and temporal gradients of clear to turbid waters, as well as oligotrophic to hypertrophic productive waters, often driven by tidal conditions. This optical variability challenges the application of standard algorithms for WQ monitoring. Therefore, the validation of such products by performing comparisons with in situ data is a crucial step for reliable analysis. Moreover, most of the bio-optical algorithms for WQ parameters retrieval are often calibrated regionally for the optical characteristics of the different sites. Thus, the performance of such algorithms outside their calibration region is also a critical topic.

In the past years, the number of WQ studies using space-born imagery over estuaries has been growing promisingly. Usually, this type of analysis in marine environments is made with ocean colour sensors, as they are specifically designed for marine purposes and have a good temporal (daily) and spatial resolution (300 m⁻¹ Km). However, for smaller regions, such as estuaries, the spatial resolution of this type of sensor could be inadequate. As an alternative, it is possible to use land sensors, as the MultiSpectral Instrument (MSI) on board the Sentinel-2 (S2) satellites, which have a finer spatial resolution (10, 20, or 60 m), even though they have a lower revisiting time (5 days considering 2-satellite constellation) [8]. As this sensor was designed primarily for land applications with bands in different wavelengths, it is important to compare different AC processors to find the most suitable one for water applications [5].

The use of Sentinel-2 MSI for WQ analysis has been showing interesting results. For example, [9] retrieved better Chl-a estimations from S2-MSI for inland and coastal waters distributed worldwide using a machine learning algorithm in comparison with state of the art ocean colour algorithms. They identified the atmospheric correction as the main hurdle in generating high-quality Chl-a products. Ansper and Alikas [5] retrieved Chl-a from S2-MSI in Estonian lakes and concluded that the sensor has potential for deriving WQ parameters from small areas, fulfilling European Union WFD monitoring and reporting requirements. In [10], the authors made use of the Sentinel-2 derived turbidity product to monitor changes in the water transparency in the Venice lagoon before and after the first COVID-19 lockdown. In Portugal, using Sentinel-2 data for this type of systems is revealing to be a good option [11,12], including for the Sado estuary [13]. However, these approaches seem to have some gaps in terms of image processing or in their framing in a water quality analysis. It is hard to find published studies based on the use of satellite imagery to discuss water quality management in the Portuguese territory, in support of water policy and legislation as the WFD. Brito et al. [14] and Cristina et al. [15] are two of the few examples and are applied to coastal waters in Western Iberia.

In the last decades, hundreds of RS publications have proposed solutions to overcome the challenges previously described and accurately quantify the biogeochemical parameters of WQ [8,16,17,18]. However, most of these analyses focus on the development of methods and do not use RS as a tool to better understand the dynamics of WQ in estuaries [8]. Additionally, the reviewed literature showed a lack of analysis that integrates several WQ parameters and not focusing only on one, such as chlorophyll-a or suspended solids as it is often observed.

In this study, we show the potential of WQ monitoring using Sentinel-2 MSI observations in the Sado estuary (Portugal), through an algorithm intercomparison exercise and time-series analysis of different WQ parameters (CDOM, Chlorophyll-a, SPM, and turbidity). This analysis was developed with the purpose of demonstrating that satellite RS can be useful for ecological assessment of estuarine systems through a water quality evaluation, studying the feasibility of using Sentinel-2 for studies in transitional waters, even if it was originally designed for land applications. Moreover, we also briefly discuss the potential and benefits of using RS methodologies to implement continuous WQ monitoring plans in transitional systems under the scope of the WFD to overcome the issue of the lack of continuous monitoring due to the high cost of in situ samplings.

2. Materials and Methods

2.1. Study Area

The Sado estuary is a transitional system located near Setúbal, Portugal. It has a total area of 212.4 km² [19], being, therefore, one of the biggest estuaries in Europe and the second largest national estuary (Figure 1). Given its high productivity, biodiversity and aesthetic value, it has been part of a natural reserve since 1980 [20]. Besides all the characteristics that make the Sado estuary a unique estuarine system, it is also a privileged place for aquaculture. The Sado oyster (Magallana angulata) is the main example of this activity in the region and is currently one of the main icons of the Setúbal region. In the context of WFD, nationally it is classified as a transitional water of typology A2- mesotidal well-mixed estuary [21].

The Sado estuary is characterized by a complex morphology and can be divided into two regions: (i) the outermost area (A in Figure 1), the Setubal bay, with 5 km wide, 20 km long and an average depth of 10m; and (ii) the innermost region (B in Figure 1), that presents a mixture of navigation channels and a vast intertidal zone of reduced depth (average depth of 10 m), relevantly influenced by the tide [22]. The connection to the ocean is made through a deep and narrow channel of approximately 1.5 km long and maximum depth of 50 m [23].

The Sado river flow can be considered low, varying on average from 0.7 m³/s in the summer to 60.0 m³/s in the winter [19,20]. Therefore, there is a strong influence of the tidal effect on the water circulation within the estuary [21,22]. The inflow and outflow of water is generally made along the two navigation channels (North and South channels, where stations 6 and 8 are located, respectively), according to the tide, with more intense currents in the South channel [24,25].

The estuary presents relatively low sun elevation during winter and higher wind intensities during spring-summer (hourly average with a maximum of ~8 m/s in 2018 and ~6m/s in 2019, Source: https://snirh.apambiente.pt (accessed on 26 November 2020)). These factors along with the tidal variation, the low water column depth and, consequently, the observed bottom resuspension, makes accurate ocean colour remote sensing challenging in this system.

2.2. In Situ Data

For the present analysis, 8 sampling stations distributed along the estuary, shown in Figure 1, were considered. The estuary was divided into two areas according to its morphologic characteristics and the variation of its physicochemical properties, as described earlier, thus stations #6, #7, and #8 were located in the outermost area (A) and stations #1 to #5 in the innermost region of the estuary (area B) (Figure 1). All sampling stations were located in optically deep waters (verified by the in situ measurements of Secchi depth). Monthly campaigns were performed between March 2018 and November 2019 within the framework of the AQUASado project. All campaigns were planned to start the water samples collection at similar tidal conditions (high water), to avoid the influence of the tide on the temporal variability of the WQ parameters, and to coincide with the S2-MSI satellite passage over the Sado estuary, when possible.

Surface water samples were collected and taken to the laboratory for quantification of a_CDOM, Chl-a, SPM, and turbidity. It is relevant to note that the total number of samples for the different parameters may vary according to the quality control results and data availability.

Moreover, the absorption coefficients of phytoplankton (a_phy) and non-algal particles (a_nap) were also obtained following the method described in [26]. These data were only used for the characterization of the study area and will be further analysed in a different work. Water temperature and salinity were also measured for site characterization purposes. The temperature was measured using an EXO2 Multiparameter Sonde (YSI) and the salinity through a high precision salinometer (Guildline Autosal 8400B) [27].

2.2.1. CDOM Absorption Coefficient

The methodology described in [26] was used to obtain the CDOM absorption coefficient (a_CDOM (λ), m⁻¹). Water samples were filtered in triplicate through a 0.2 µm pore size polycarbonate filter (Whatman). The absorbance readings of the filtrate were obtained at the range of 300–800 nm through a dual-beam spectrophotometer (Shimadzu 2600). From the average of the absorbance readings and the subtraction of the reference sample, the values obtained at each λ were multiplied by the light extinction coefficient (2.303) and divided by the cuvette length (0.01 m) [26].

2.2.2. Chlorophyll-a

Water samples were filtered using GF/F filters with 25 mm diameter and a pore of 0.7 μm (Whatman). The filters were then wrapped in aluminium and stored at −80 °C to be later analysed through high performance liquid chromatography (HPLC) (C18 column) to quantify the different phytoplankton pigments following the method described in [28] and adapted by [29]. For the extraction, 3 mL of a solution composed of 95% cold methanol, buffered with 2% of ammonium acetate and containing 0.35 mgL⁻¹ of trans-β-Apo-8′-carotenal (as internal standard) were used. After placing the extraction solution in the sample test tube, the filter was grinded with a glass rod. Then, the tube was placed in the freezer (−20 °C) for 30 min after which it was placed in the ultrasound for 5 min. The test tube was then put back in the freezer (−20 °C) for another 30 min and, just after, centrifuged for 10–15 min at 4000 RPM (4 °C). Lastly, the sample was filtered and run in the HPLC machine.

2.2.3. Suspended Particulate Matter

SPM, defined as the dry mass of particles per unit volume of water sample (mg L⁻¹), was determined following the standard gravimetric method described in [30]. A known volume of water sample was filtered through pre-ashed and pre-weighted Whatman GF/F filters (nominal pore size 0.7 µm and 47 mm diameter), with volume depending on the amount of material present in the water. At the end of filtration, the filters, including the outer edge of the filters, and filtration apparatus were rinsed with MilliQ water to remove any salt. The dry mass of particles collected on the filters was then measured with a precision balance.

2.2.4. Turbidity

The turbidity was determined with a portable infrared turbidimeter (Lovibond TB 210 IR) which quantifies the turbidity in nephelometric turbidity units (NTU) by measuring the light (λ = 860 nm) scattered at 90°, with an accuracy of ±2.5% of reading. The instrument calibration was checked before each sampling campaign using the calibration standards provided with the instrument (<0.1, 20, 200, and 800 NTU) and track the instrument stability. The turbidity was always recorded in triplicate for the same water sample and averaged.

2.3. Sentinel-2 MSI Data

Sentinel-2 (A/B) MSI Level 1C images were used in the present study. The S2-MSI granules acquired between March 2018 and March 2020 were downloaded from Sentinel Scientific Data Hub (https://scihub.copernicus.eu/ (accessed on 26 November 2020)). Cloud free L1C images with coincident field campaign dates were used to compare satellite-derived data with in situ measurements or match-ups (processing chain shown in Figure 2). The images were re-sampled to 10 m as part of the pre-processing and the average of 3 × 3 pixels window centred on the sampling station was used for the in situ comparison and to perform the time-series.

To perform the matchups, the available images were filtered considering the time difference between the satellite passage and the in situ measurements. The time window chosen was ±2 h. The number of total match-ups depended on the AC processor used. To ensure high quality of satellite data, different quality flags were considered (

Table 1. Flags used from AC-C2RCC and Polymer.

AC Processor	Flag	Meaning
C2RCC	Rtosa_OOS	The input spectrum to the atmospheric correction neural net was out of the scope of the training range and the inversion is likely to be wrong.
	Rtosa_OOR	The input spectrum to the atmospheric correction neural net out of training range
	Rhow_OOS	The Rhow input spectrum to the Inherent Optical Properties (IOP) neural net is probably not within the training range of the neural net and the inversion is likely to be wrong.
	Rhow_OOR	One of the inputs to the IOP retrieval neural net is out of training range
	Cloud risk	High downwelling transmission indicates cloudy conditions
Polymer	!bitmask & 1023 = 0	Invalid pixel

). These quality flags are part of the L2 processing and, when raised, they indicate that something might have gone wrong during the data processing and that the retrieved product should be carefully analysed. Moreover, the IdePix (v2.2) was used on L1C images for pixel type identification.

Given the lack of in situ radiometric measurements for the direct test of AC, three AC processors were evaluated by testing a combination of AC + bio-optical algorithms (Figure 2). The AC processors tested were Acolite v20190326 (RBINS) [31], C2RCC v2.0 [32] and Polymer v4.12 [33]. The selection of the aforementioned AC algorithms was based on their free availability, and on the correcting effects that could be applied, such as sun glint and adjacency of land pixels, common sources of errors in transitional waters.

Acolite is an AC processor developed for coastal and inland waters. By default, it performs atmospheric correction using the dark spectrum fitting approach, but it can be configured to use an exponential extrapolation approach [34].

The C2RCC processor relies on a database of radiative transfer simulations of water-leaving reflectance and related top-Of-atmosphere radiances (satellite signal) [32]. The inversion of water signal and satellite signal is performed by neural networks and is fully described in [32].

Polymer AC is a coupled water-atmosphere model that uses a spectra-matching optimization approach to obtain a best-fitting combination of water and atmosphere unknowns, being one of its strengths to retrieve water reflectance in the presence of sun glint. This AC does not use exclusively the signal in the NIR, but it also uses all the available spectral bands in the visible [33].

It is important to note that these processors are all under active development but they are publicly available and are considered mature and useful to compare [35]. Each processor was used with its default settings, as this should be the best option for general use, without a priori knowledge of the water body or atmospheric conditions.

Once the image processing was complete, the atmospherically corrected outputs, i.e., water reflectance ρ_w or remote sensing reflectance R_rs (R_rs = ρ_w/π), were used to test different algorithms for WQ parameters retrieval (

Table 2. Bio-optical algorithms tested for the different water quality parameters.

Algorithm	Equation	Type of Algorithm	MSI Bands Involved (nm)	Reference
acdom(m−1)
TS443	$acdom\left(443\right)=0.00129+0.6543\times \left(\frac{Rrs670}{Rrs490}\right)$	Empirical	490, 665	[36]
TS412	$acdom\left(412\right)=0.00411+2\times \left(\frac{Rrs670}{Rrs490}\right)$	Empirical	490, 665	[36]
MAs	$acdom$ (443) = ln[($\left(\frac{Rrs490}{Rrs551}\right)$− 0.4247)/2.453]/(−13.586)	Empirical	490, 560	[37]
MAM	$acdom$ (443) = ln[($\left(\frac{Rrs490}{Rrs551}\right)$− 0.4363)/2.221]/(−13.126)	Empirical	490, 560	[37]
CZ	$acdom$ (440) = 0.2987x–1.369, x = Rrs B1/Rrs B3	Empirical	490, 665	[38]
CH	$acdom$ (440) = 28.966e−2.015x, x = Rrs 560 nm/Rrs 665 nm	Empirical	560, 665	[39]
KU	$acdom\left(420\right)=5.20\times {\left(\frac{Rrs565}{Rrs660}\right)}^{-2.76}$	Empirical	560, 665	[40]
Chlorophyll-a
C2RCC	-	Neural network		[32]
OC3	$OC3=lo{g}_{10}\left[\left(\left[R_{rs}\left(442\right),R_{rs}\left(492\right)\right] \right)\times R_{rs}{\left(560\right)}^{-1}\right)$ $y=0.3308-2.684x+1.5990x^2+0.5525x^3-1.4876x^4$ $Chl{a}_{OC3}={10}^y$	Empirical	443, 490, 560	[17]
2-Band	$Chl{a}_{2Band}={\left(35.75\times \frac{R_{rs}\left(708\right)}{R_{rs}\left(665\right)}-19.3\right)}^{1.124}$	Semi-analytical	665, 705	[41]
Gons et al. (2005) (GS)	$b_b\left(783\right)=1.56\times {\rho }_w\left(783\right)/\left[0.082-0.6\times {\rho }_w\left(783\right)\right]$ $a_{phy}\left(665\right)=\left(0.70+b_b\right)\times {\rho }_w\left(705\right)/{\rho }_w\left(665\right)-0.4-b_b^p$ $Chl{a}_{Gons}=a_{phy}\left(665\right)/a_{phy}^{*}\left(665\right)$	Semi-analytical	665, 705, 783	[42]
Suspended Particulate Matter
C2RCC	-	Neural network		[32]
Nechad	$SPM=\frac{A^{\lambda }{\rho }_w\left(\lambda \right)}{1-\left(\frac{{\rho }_w\left(\lambda \right)}{C_{\lambda }}\right)}+B_{\lambda }$	Semi-analytical, single band	665, 705, 740, 783, 865	[43]
Siswanto	$lo{g}_{10}\left(SPM\right)=0.649+25.623X_1-0.646X_2$	Empirical	560, 665, 490	[44]
Turbidity
Dogliotti	$Turb=\frac{A^{\lambda }{\rho }_w\left(\lambda \right)}{1-\left(\frac{{\rho }_w\left(\lambda \right)}{C_{\lambda }}\right)}$	Semi-analytical, single band	665, 865	[45]
Nechad	$Turb=\frac{A^{\lambda }{\rho }_w\left(\lambda \right)}{1-\left(\frac{{\rho }_w\left(\lambda \right)}{C_{\lambda }}\right)}+B_{\lambda }$	Semi-analytical, single band	665, 705, 740, 783, 865	[46]

), with the aim of selecting the best processing chain for the study area. Several algorithms have been developed to retrieve WQ parameters from S2-MSI data. However, the majority of them are often calibrated regionally for the optical characteristics of the different sites. Thus, the performance of such algorithms outside their calibration region is a key topic. To address this matter, an algorithm intercomparison exercise was performed for the different parameters under evaluation. The selection of the bio-optical algorithms to test was conducted considering: (i) the similarity between the original ranges of the in situ parameters used for their calibration and the range of in situ parameters observed in this study; (ii) their good performance outside their calibration region based on the literature review; and (iii) their applicability to S2-MSI data.

2.4. Statistical Indicators

To assess the agreement between the in situ observations and the satellite data, a statistical evaluation was performed using: linear regression parameters including coefficient of determination (R²), slope and intercept and uncertainty estimates including: root mean square deviation (RMSE), bias error (BIAS, δ), unbiased root mean squared error (URMSE, ∆), absolute percentage difference (APD), and relative percentage difference (RPD). The RMSE, the bias, and the APD were used to evaluate the accuracy of the match-ups and the RPD was used as an indicator of the systematic error (relative bias). These parameters can be defined as follows:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{ \left(sa{t}_i-in sit{u}_i\right)}^2}{N}}$$

(1)

$$BIAS=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\left(sa{t}_i-in sit{u}_i\right)$$

(2)

$$URMS =\sqrt{\frac{1}{N}{\displaystyle \sum }_{1=1}^N{\left[\left(sa{t}_i-{\mu }_{sat}\right)-\left(in sit{u}_i-{\mu }_i\right)\right]}^2}$$

(3)

$$APD(\%)=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\frac{\left|sa{t}_i-in sit{u}_i\right|}{in sit{u}_i}\times 100$$

(4)

$$RPD(\%)=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\frac{sa{t}_i-in sit{u}_i}{in sit{u}_i}\times 100$$

(5)

where Sat is the satellite-derived data and in situ the observations collected in the estuary. N is the total number of samples, i is the sample index and ${\mu }_{sat}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$N$}\right.{\displaystyle \sum }_{i=1}^{N}sa{t}_i$– analogous definition applies for µ_{in situ}.

Statistical analyses of Chl-a are presented in logarithmic scale to account for the log-normal distribution of the bio-optical parameter [47]. According to [48] the RPD and APD were calculated without log-transforming the Chl-a values.

Additionally, summary diagrams such as Taylor and Target plots, referenced in [49] and [50], were used for a better visualization and inter-comparison of results. The Taylor diagram is based on the scoring of the normalized standard deviation (STD_sat/STD_{in situ}, [51]) and correlation coefficient (R). The Target diagram summarizes the scoring of the URMS and BIAS instead.

2.5. Time-Series Analysis

The statistical indicators described in Section 2.4 were used to select the best processing chain from the match-ups intercomparison exercise and perform a spatio-temporal analysis of the different WQ parameters, with the main objective of evaluating the capabilities of S2-MSI imagery for continuous monitoring of WQ variability in estuarine systems.

The best performing combination of AC + bio-optical algorithms was applied to all cloud-free L1C images sensed between March 2018 and March 2020. As done for the match-ups exercise, quality flags were considered to ensure the highest quality of satellite data.

Given the strong variability of the in situ data and the high dynamic that characterizes the Sado estuary, the time-series analysis was performed separately for the region A (outermost area) and the region B (innermost area) (Figure 1). The time-series performed with RS products and in situ observations considered the average of the pixels corresponding to the sampling stations included in each of the areas. Moreover, the seasonal variability (astronomical seasons) of the WQ parameters in the whole estuary and its adjacent area, is provided through high resolution maps.

3. Results

3.1. Study Area Characterization

Based on the monthly surface data collected in situ between March 2018 and February 2019 analysed in the present study (see Section 2.2), the estuary presents clear seasonal patterns of the physical and chemical properties (

Table 3. Annual and seasonal averages of the temperature (T), salinity, chlorophyll-a (Chl-a), suspended particulate matter (SPM), colored dissolved organic matter absorption (a_CDOM), and turbidity in the Sado estuary (based on the monthly surface data collected in situ between March 2018 and February 2019). A and B are the outermost and innermost regions of the estuary, respectively.

		T (°C)		Salinity		Chl-a (mg/m3)		SPM (mg/L)		aCDOM 443 nm (m−1)		Turbidity (NTU)
		Mean	STD	Mean	STD	Mean	STD	Mean	STD	Mean	STD	Mean	STD
Winter	A	14.3	0.9	35.89	0.40	0.90	0.26	1.14	0.46	0.21	0.05	1.10	0.33
	B	12.9	2.4	32.45	1.95	1.45	0.92	2.61	1.13	0.38	0.07	2.14	0.81
	Estuary	13.4	2.1	33.76	2.29	1.25	0.78	2.19	1.19	0.36	0.10	1.84	0.85
Spring	A	15.5	1.2	35.02	0.89	2.09	1.33	1.71	0.75	0.16	0.09	-	-
	B	17.4	2.0	29.13	5.13	3.54	5.85	8.40	7.74	0.37	0.12	-	-
	Estuary	16.7	1.9	30.91	5.04	3.03	4.79	6.09	7.00	0.26	0.15	-	-
Summer	A	18.7	1.0	36.02	0.27	1.66	0.35	2.06	1.06	0.15	0.12	1.30	0.86
	B	23.3	0.9	35.23	1.42	4.21	2.41	8.37	4.42	0.47	0.15	6.58	4.49
	Estuary	17.7	4.0	35.46	1.25	3.34	3.66	5.74	4.94	0.31	0.20	4.28	3.86
Autumn	A	17.7	2.2	35.91	0.15	1.47	1.29	1.80	0.46	0.14	0.09	1.31	0.49
	B	20.6	4.3	34.14	2.83	1.59	0.53	3.91	1.69	0.42	0.11	3.03	1.38
	Estuary	18.9	3.3	34.73	2.42	1.55	0.81	3.21	1.72	0.30	0.18	2.46	1.41
Year	A	16.1	2.3	35.57	0.78	1.70	1.29	1.55	0.73	0.16	0.10	1.11	0.44
	B	17.4	4.2	32.31	4.30	3.40	3.75	5.99	4.95	0.41	0.13	4.59	3.92
	Estuary	17.0	3.7	33.36	3.87	3.05	3.43	4.97	4.74	0.29	0.17	3.81	3.75

) with a clear tendency to find higher values of a_CDOM, Chl-a, SPM, and turbidity during spring and summer. Furthermore, the estuary should not be considered spatially homogeneous. The innermost region of the estuary (B) shows higher values of the optically active components and a greater susceptibility to seasonal variations of the weather, given the greater thermal amplitude observed.

The Sado estuary is optically dominated by CDOM, as observed in Figure 3, with higher values of a_CDOM compared in respect to phytoplankton and non-algal particle absorption coefficients, a_phy and a_nap respectively. Figure 2 shows the dominance of a_CDOM mainly in stations located in area B, showing smaller contributions of a_nap and a_phy. It should be reminded that a_nap and a_phy are only here presented to characterise the study area and will not be further analyzed in the present study. In region A, the contribution of a_phy was slightly higher. However, a_CDOM remains dominant in most of the stations located in this area. The important contribution of the a_CDOM (443 nm) even in stations located in area A, allows to classify the whole estuary as optically complex with Case 2 waters [52], regardless of the studied period and the tidal conditions.

3.2. S2-MSI Match-Ups

For this study, the combination of AC and bio-optical algorithms were evaluated against in situ values of the different WQ parameters. An algorithm intercomparison exercise was performed for the available match-ups with the objective of selecting the best processing chain for a spatio-temporal analysis of the different WQ parameters under investigation. Taylor and Target diagrams were used for the selection of the best performing processing chain.

Even if no in situ radiometric measurements were available for a thorough comparison, the ${\rho }_w$ spectra derived from the matchups of the three AC processors were compared separating the two regions of the estuary (Figure 4).

Figure 4 shows that Acolite retrievals are significantly higher across the whole spectrum compared to the other ACs. In spite of this, similar shapes were obtained by the three processors. From a visual inspection, it can be noted a difference between the spectral shapes of areas A and B. In particular, for area B, ${\rho }_w$ in the blue, is lower compared to the green region of the spectra, i.e., the green peak is more pronounced than in area A. This could be partially due to higher importance of CDOM (a_CDOM) and sediments (a_nap), which agrees with in situ measurements (Figure 3). Moreover, it can be observed that reflectance in the NIR region of the spectra is higher in area B compared to A, which matches the higher SPM found in the inner part of the estuary.

Regarding a_CDOM, AC-C2RCC and Polymer showed the best performance (mean APD = 79% and APD = 68%, respectively) when combined with the bio-optical algorithms, compared to the results obtained using Acolite (mean APD = 164%). In the Taylor diagram (Figure 5), the TS443 and KU algorithms combined with Polymer (pTS443 and pKU) were closer to the reference (represented by a star). Considering the AC-C2RCC processor, the KU, CZ, and TS443 were the best combination performing algorithms. In the Target plot, the TS443 algorithm showed the results closest to the reference (zero), both combined with Polymer and AC-C2RCC. All other processing chains (except AC-C2RCC) underestimated the a_CDOM. Among the algorithms that showed the best performance, the TS443 was the one with the slope values closest to 1 (0.9 with AC-C2RCC and 0.8 with Polymer), and the lowest RMSE (0.151 and 0.106, respectively for AC-C2RCC and Polymer) compared to the other algorithms (

Table 4. Summary of the results obtained in the application of the different algorithms to images processed with C2RCC, Polymer and Acolite, for the quantification of CDOM, chlorophyll-a, SPM, and turbidity. The table contemplates the data of the average of 3x3 pixels centred in the pixel of the sampling stations. The best performing processing chain of each parameter are highlighted in bold.

	AC Processor	Algorithm	Equation	R2	RMSE	URMS	BIAS	APD (%)	RPD (%)	N
CDOM	AC-C2RCC	C2RCC (443 nm)	Y = 0.4x + 0.22	0.201	0.170	0.167	0.032	99.019	70.360	21
		TS (443 nm)	Y = 0.9x + 0.05	0.471	0.151	0.149	0.027	63.571	40.000	21
		TS (412 nm)	Y = 1.6x + 0.18	0.563	0.701	0.447	0.540	143.856	134.141	21
		MAS (443 nm)	Y = 0.1x + 0.01	0.427	0.319	0.142	−0.285	78.199	−78.199	20
		MAM (443 nm)	Y = 0.3x + 0.1	0.325	0.182	0.124	−0.133	63.744	−8.546	20
		CZ (443 nm)	Y = 0.4x − 0.004	0.462	0.226	0.117	−0.194	57.911	−51.400	21
		CH (443 nm)	Y = 0.6x − 0.08	0.294	0.265	0.169	−0.204	74.368	−68.336	21
		KU (420 nm)	Y = 0.7x − 0.05	0.446	0.283	0.199	−0.201	50.334	−40.553	21
	Polymer	TS (443 nm)	Y = 0.8x + 0.06	0.608	0.106	0.105	−0.013	46.612	19.630	21
		TS (412 nm)	Y = 1.4x + 0.21	0.651	0.516	0.303	0.418	106.388	101.800	21
		MAS (443 nm)	Y = 0.1x + 0.01	0.607	0.342	0.152	−0.306	85.400	−85.400	20
		MAM (443 nm)	Y = 0.2x + 0.09	0.515	0.229	0.133	−0.187	68.307	−27.708	21
		CZ (443 nm)	Y = 0.3x + 0.003	0.582	0.245	0.114	−0.217	62.152	−60.788	21
		CH (443 nm)	Y = 0.7x − 0.07	0.376	0.220	0.159	−0.152	66.332	−55.218	21
		KU (420 nm)	Y = 0.8x − 0.04	0.529	0.231	0.187	−0.135	41.981	−31.224	21
	Acolite	TS (443 nm)	Y = 0.4x + 0.3	0.561	0.142	0.106	0.095	117.027	114.796	19
		TS (412 nm)	Y = 0.8x + 0.89	0.672	0.787	0.165	0.770	258.640	258.640	19
		MAS (443 nm)	Y = 0.1x + 0.09	0.620	0.261	0.139	−0.222	74.266	−35.156	19
		MAM (443 nm)	Y = 0.1x + 0.09	0.620	0.262	0.138	−0.222	73.605	−36.073	19
		CZ (443 nm)	Y = 0.2x + 0.1	0.565	0.213	0.123	−0.174	70.461	−18.879	19
		CH (443 nm)	Y = 0.6x + 0.31	0.153	0.578	0.291	0.499	306.715	304.302	19
		KU (420 nm)	Y = 0.6x + 0.8	0.182	0.673	0.310	0.598	244.380	243.929	19
Chlorophyll-a	AC-C2RCC	Chla-C2RCC	Y = 4.9x − 1.33	0.223	0.714	0.672	0.242	281.060	243.552	13
		OC3	Y = 11.0x − 2.49	0.259	0.885	0.519	0.717	749.362	733.932	13
		GS	Y = 0.4x − 0.07	0.625	0.565	0.205	−0.527	66.985	−66.985	13
		2-Band	Y = 4.4x + 1.6	0.608	0.765	0.126	0.754	493.038	493.038	13
	Polymer	Polymer	Y = 2.5x + 0.67	0.365	0.503	0.228	0.448	221.230	216.124	13
		OC3	Y = 2.8x + 1.27	0.294	0.603	0.260	0.543	309.703	304.136	13
		GS	-	-	-	-	-	-	-	-
		2-Band	-	-	-	-	-	-	-	-
	Acolite	OC3	Y = 2.6x + 0.63	0.217	0.5151	0.229	0.461	224.311	224.311	11
		GS	Y = −0.2x + 1.14	0.0004	0.327	0.306	−0.115	61.556	−1.586	11
		2-Band	Y = −1.4x + 13.07	0.0003	1.044	0.233	1.018	1116.360	1116.360	11
SPM	AC-C2RCC	SPM-C2RCC	Y = 2.2x − 0.59	0.467	25.696	3.927	3.025	107.691	88.159	23
		N (665 nm)	Y = 0.8x + 1.82	0.443	1.688	1.344	1.020	56.509	46.041	23
		N (705 nm)	Y = 0.8x + 1.03	0.470	1.343	1.307	0.3111	35.319	17.718	23
		N (740 nm)	Y = 0.8x + 0.91	0.493	1.259	1.248	0.170	31.814	12.884	23
		N (783 nm)	Y = 0.8x + 1.21	0.498	1.443	1.292	0.642	43.561	30.144	23
		N (865 nm)	Y = 0.5x + 2.11	0.503	1.237	1.079	0.606	47.695	37.523	23
		SI	Y = 0.7x + 0.62	0.409	1.416	1.373	−0.343	32.926	−8.814	23
	Polymer	N (665 nm)	Y = 0.7x + 2.12	0.536	1.611	1.102	1.175	57.332	53.049	23
		N (705 nm)	Y = 0.7x + 0.85	0.568	1.055	1.018	−0.274	26.082	−2.703	23
		N (740 nm)	Y = 0.4x + 2.17	0.229	1.427	1.421	0.129	42.072	18.463	23
		N (783 nm)	Y = 0.8x + 2.94	0.475	2.521	1.279	2.173	95.546	94.803	23
		N (865 nm)	Y = 0.7x + 0.94	0.253	1.593	1.540	0.409	57.568	6.599	21
		SI	Y = 0.5x + 1.05	0.589	1.198	1.002	−0.658	21.122	−14.695	23
	Acolite	N (665 nm)	Y = 1.1x + 3.57	0.489	1.256	1.198	3.849	153.309	153.309	20
		N (705 nm)	Y = 1.2x + 3.28	0.516	1.391	1.230	3.838	152.629	152.629	20
		N (740 nm)	Y = 2.0x + 10.27	0.196	2.409	1.513s	13.444	526.817	526.817	20
		N (783 nm)	Y = 1.6x + 11.81	0.133	2.829	1.759	13.804	554.024	554.024	20
		N (865 nm)	Y = 2.2x + 15.72	0.503	1.496	1.330	19.432	780.628	780.628	20
		SI	Y = 0.7x + 1.75	0.548	0.928	0.928	0.644	35.952	31.502	20
Turbidity	AC-C2RCC	N (665 nm)	Y = 0.5x + 1.00	0.762	1.563	1.508	−0.412	28.999	−0.165	15
		N (705 nm)	Y = 0.5x + 0.77	0.805	1.581	1.450	−0.630	29.688	−11.269	15
		N (740 nm)	Y = 0.6x + 0.38	0.837	1.706	1.393	−0.985	37.970	−31.209	15
		N (783 nm)	Y = 0.6x + 0.43	0.841	1.479	1.280	−0.742	33.400	−22.474	15
		N (865 nm)	Y = 0.4x + 0.23	0.838	2.344	1.717	− 1.597	53.064	−53.064	15
		DO	Y = 0.5x + 0.56	0.793	1.790	1.522	−0.943	37.657	−26.332	15
	Polymer	N (665 nm)	Y = 0.5x + 1.43	0.692	1.714	1.699	−0.220	31.803	15.192	15
		N (705 nm)	Y = 0.4x + 0.58	0.761	2.093	1.717	−1.198	44.136	−36.130	15
		N (740 nm)	Y = 0.3x + 1.05	0.806	2.245	1.966	−1.085	38.869	−14.188	15
		N (783 nm)	Y = 0.4x + 2.12	0.462	2.029	2.009	0.290	69.471	53.834	15
		N (865 nm)	Y = 0.4x + 0.44	0.656	2.716	1.957	−1.883	51.924	−49.191	11
		DO	Y = 0.4x + 0.95	0.715	1.893	1.725	−0.781	29.646	−13.270	15
	Acolite	N (665 nm)	Y = 0.2x + 1.48	0.747	2.577	2.370	−1.013	47.825	3.754	13
		N (705 nm)	Y = 0.2 + 1.55	0.744	2.570	2.385	−0.957	50.436	9.282	13
		N (740 nm)	Y = 0.2x + 4.14	0.282	2.983	2.464	1.680	182.207	173.672	13
		N (783 nm)	Y = 0.2x + 4.38	0.258	3.137	2.485	1.915	199.765	193.217	13
		N (865 nm)	Y = 0.2x + 5.86	0.186	4.319	2.567	3.474	296.190	292.246	13
		DO	Y = 0.1x + 0.97	0.340	2.858	2.418	−1.742	53.563	−30.192	13

).

For Chl-a, combining the results obtained with the Taylor diagram and the target plot, the bio-optical algorithm/AC processor combination that compared best with the in situ data was GS with AC-C2RCC (cGS) reflectances. For this processing chain, the highest determination coefficient (R² = 0.63) was obtained, while the remaining statistical parameters (RMSE, bias, APD and RPD) were the closest to zero (Figure 6). It should be noted that the application of 2-Band algorithm [41] to AC-C2RCC images also produced good results (the second best), superior to the ones of the standard Chla-C2RCC algorithm and even of OC3 (Taylor plot, Figure 5). On the other hand, the standard algorithm from Polymer (pPolymer) showed better results than the standard algorithm from C2RCC, but the first processor failed in the quantification of Chl-a using GS and 2-Band algorithms. In general, the satellite-derived data tended to overestimate the in situ data (except the combination with the GS algorithm). The retrieval of Chl-a remains relatively inconsistent, with large differences between satellite and in situ observations (high APD values).

For the SPM and turbidity products, better results were obtained when NIR bands were used. Considering the rather low measured SPM and turbidity values, signal in the NIR was expected to be low and a better correlation with the red bands could be expected. From the S2 extracted spectra (Figure 4) this is clear for region A, but less evident in region B. However, a poor relationship was found using the red bands for region A (not shown), which could be partially due to the overlap of this band with the location of the Chl-a absorption peak which might be affecting the retrieval of SPM and turbidity at this wavelength. Considering this and given the low number of match-ups available for region A, the dataset was not separated by region and thus, the general performance of the algorithm was used to select the best performing AC+algorithm. It should be noted though that for the clearer waters in region A the correlation is low (blue signs in Figure 6), thus results using the NIR bands should be considered with caution for this area.

For the SPM product, the best results were obtained from Polymer in combination with Nechad algorithm [44] at 705 nm (pN705) (R² = 0.57, APD = 26%). However, this algorithm decreased its performance when applied to images processed with the AC-C2RCC. On the other hand, with AC-C2RCC, good results were also obtained, using Nechad algorithm for 740 nm (cN740) (R² = 0.49, APD=31%) (Figure 6). Interestingly, the same algorithm (for the same wavelength) showed a worse performance when applied to images processed with Polymer (Table 4). For SPM, Acolite + bio-optical algorithms were the combinations that presented the worst performance (APD always > 150% except from aS, where APD = 31%) and the only ones in which there was a tendency to overestimate in situ values, suggesting that ${\rho }_w$ is overestimated compared to other AC (Figure 4). From the set of results, it was the application of Nechad for the 665 and 705 nm wavelengths that showed a greater consistency in the results between the different processors, although they did not have the overall best performance when compared with the in situ dataset.

The highest agreement between the satellite and the in situ turbidity data was achieved by applying Nechad at the wavelength of 783 nm to the images processed with AC-C2RCC (cN783, R² = 0.84, APD = 33%) (Figure 6). Overall, this AC processor presented good results for the set of algorithms and bands tested. Additionally, AC-C2RCC showed a tendency to underestimate satellite-derived turbidity data in comparison to the in situ measurements (negative RPD for all the algorithms). The application of the Dogliotti algorithm [46] generated results similar to those of Nechad for AC-C2RCC images, as expected. With Acolite, the application of Nechad to bands 740, 783 and 865 nm showed a relevant decrease in the quality of the match-ups results, when compared with the lower wavelengths.

Overall, two out of the three tested AC processors algorithms, AC-C2RCC and Polymer plus bio-optical, seemed to provide results more comparable with in situ data for the selected study area (Table 4). The combinations (AC+bio-optical algorithms) tested with these two processors presented similar results for all the parameters and can be considered more suitable for the study area than the combinations using Acolite (mean APD = 254%), which showed a general overestimation of the in situ values (mean Bias_ACOLITE = 2.78; mean Bias_C2RCC = 0.04; mean Bias_POLYMER = −0.07). Even though Polymer combinations have outperformed AC-C2RCC combinations in specific cases, such as in the quantification of CDOM or SPM, it did not present satisfactory results for Chl-a. Considering that AC-C2RCC combinations showed good capabilities for all these parameters, the best performing AC-C2RCC processor plus bio-optical algorithms were selected for the time-series analysis for being the best combination for the set of parameters analysed. In accordance, only the best match-ups results of the application of the algorithms to AC-C2RCC images will be presented in this section (Figure 6). The full dataset of statistical results is available in Table 4.

Since the combinations using the AC-C2RCC were chosen as overall best performing, the algorithms selected to perform the time-series analysis were the following: cT443 (a_CDOM), cGS (Chl-a), cN740 (SPM), and cN783 (turbidity). The table contemplates the data of the average of 3 × 3 pixels centred in the pixel of the sampling stations. The best performing processing chain of eaxh parameter are highlighted in bold.

Overall, for the set of parameters tested, encouraging results were obtained in the satellite data intercomparisons process, which indicates the usefulness of MSI images to monitor water quality parameters in the Sado estuary.

3.3. Spatio-Temporal Analysis: Sado Estuary Case Study

To demonstrate and evaluate the usefulness of a consistent Sentinel-2-derived WQ data record, we present here a spatio-temporal analysis of the four WQ parameters under investigation for the Sado estuary. To do so, the best performing processing schemes chosen during the intercomparison part of the present study were selected and applied to S2-MSI time-series, i.e., AC-C2RCC-TS443 for a_CDOM, AC-C2RCC-GS for Chl-a, AC-C2RCC-N740 for SPM, and AC-C2RCC-N783 for turbidity (Figure 6). Moreover, the time-series analysis was performed separately for the two different regions of the estuary, given the high optical variability that characterizes this transitional system.

Figure 7 and Figure 8 show the time-series of the different parameters considering the spatially averaged observations. From a visual inspection, it is hard to observe a clear seasonal variability of any of the selected parameters in area A. Moreover, region B presents a greater amplitude (over a short time, due to tidal variability) of values in respect to region A, which remains relatively constant. This might lead to bigger uncertainties associated with satellite-derived products over region B due to higher optical variability and greater interaction between the optically active compounds. Therefore, a worse performance of the bio-optical algorithms is expected. Moreover, the greater amplitude of values will lead to difficulties in detecting WQ anomalies without a consistent climatological record of the area of interest. On the other hand, S2-MSI time-series may allow us to better understand the short-term variability which would be highly laborious through in situ observation.

For a better visualization of the seasonal patterns, high spatial resolution seasonal maps of the full estuarine system are provided in Figure 9. For all the set of parameters, the highest values were found in Spring and Summer, mainly in the inner region of the estuary, following the tendency shown by the field measurements.

Regarding the a_CDOM (443 nm), the seasonal maps show higher values during Spring within the whole estuary. However, the seasonal average of in situ a_CDOM (Table 3) shows that the highest values were obtained in Summer in the innermost region of estuary, and in winter in region A. In this region, the outermost area, the greatest differences between the seasonal averages of satellite-data and in situ (Table 1) were found in summer (0.14 m⁻¹) and in autumn (0.11 m⁻¹), while spring and winter presented the smallest variations (0.03 and 0.05 m⁻¹ respectively). In region B, the greatest differences were in spring (0.28 m⁻¹) and in winter (0.20 m⁻¹), while the smallest in summer (0.11 m⁻¹) and during autumn (0.16 m⁻¹). Even if the satellite seems to be able to follow the in situ trend in both areas, the satellite-derived values slightly overestimate the field observations most of the time. This can be clearly seen in the time-series plot, especially in region B between April and October 2019.

For chlorophyll-a, it was possible to observe a significant difference (mean APD = 391%) in terms of absolute values between the in situ data and the satellite-derived ones (Figure 7). This deviation between datasets seemed to be more evident during Spring and Summer, both in region A and region B. It is possible to observe that the satellite-derived data were not able to match the highest chlorophyll in situ values found during the sampling period, presenting a significant underestimation in comparison to the in situ observations. During Autumn and Winter, when field Chl-a measurements were lower, a greater similarity between the datasets was found. In general, region B presented slightly higher Chl-a values than region A, as would be expected, given the in situ variability observed. As such, the MSI sensor was sensitive to the spatial variability of the chlorophyll-a in the Sado estuary. This spatial variability of Chl-a becomes more evident when looking at the entire estuary and its surrounding region (Figure 9). It should be noted that in Summer, it was not possible to obtain valid values for the region in the coastal ocean (Figure 9) due to sun glint occurrences that led to an overcorrection in the atmospheric correction part of the processing, which, combined with low concentrations of chlorophyll, resulted in invalid Chl-a values. During winter, it was also possible to observe high chlorophyll values in the interior of the estuary and in the region of the coastal ocean (higher values than during spring).

Regarding SPM, Figure 7 and Figure 8 show that the derived product can clearly capture the temporal variability and potential anomalies in water conditions. As expected, higher values were observed in the inner region of the estuary at all seasons, being spring and summer the seasons with the highest levels of suspended solids within the estuary. As well as for a_CDOM and turbidity, the SPM retrieved recurring to the S2-MSI presented a slight overestimation of the in situ observations.

As already mentioned in the match-ups section, the turbidity product was the best represented WQ parameter using S2-MSI data. This can be easily observed also from the time-series plots, where the satellite-derived turbidity followed the in situ observations with high similarity, showing the great potential of this satellite for routine monitoring of this important parameter. Again, the highest turbidity levels were reached in spring and summer, with higher concentrations in the inner region of the estuary. It is important to note that despite the match-ups stations being optically deep at the time of the sampling, for the time-series it was not possible to know it in advance. Moreover, in this study, it was not possible to assess if the shallower regions of the estuary were influenced by the bottom signal therefore Figure 9 should be interpreted with caution.

4. Discussion

4.1. Algorithms Intercomparison Exercise

The principal objective of this study was to evaluate the potential use of high spatial resolution satellite data to estimate water quality (WQ) parameters in a highly dynamic system. For this purpose, different processing chains (AC plus WQ algorithms) have been selected and applied to S2-MSI imagery in order to evaluate the usefulness of this satellite for water quality monitoring in the Sado estuary.

According to the results, the statistics indicate consistency between the satellite-derived values and in situ observations of turbidity, indicating that this parameter is the best WQ parameter to be retrieved via S2-MSI data. However, the agreement is also relatively strong for SPM and a_CDOM for the study area, while the Chl-a product showed some limitations and higher errors associated (Table 4).

Even if no field radiometric measurements were available to directly assess the performance of the AC processors, AC-C2RCC and Polymer combined with bio-optical algorithms were the combinations that provided the most comparable results to the in situ data. Similar results have also been discussed by Warren et al. [35], where different AC processors have been tested on S2-MSI data showing that AC-C2RCC and Polymer outperformed Acolite for applications in coastal waters. Pereira-Sandoval et al. [53] also showed similar results, pointing out the good performance of these two ACs in Spanish inland waters and suggesting the application of a water type classification to improve the performance of the processors.

Polymer and AC-C2RCC were the processor that presented the smallest reflectance values between the processors tested in this study. An underestimation from the Polymer processor was observed in other studies. Steinmetz et al. in [33,54], the authors pointed towards a miss-estimation of the wind speed from the afore-mentioned processor that might lead to an over-correction of the model, resulting in a significant reduction of the total signal, which could explain the invalid Chl-a values obtained using this processor. Again, the lack of in situ radiometric measurements did not permit to deeply explore the performance of the ACs in this study area. Even if the Polymer processor, in combination with bio-optical algorithms, outperformed AC-C2RCC in some cases, the neural network methodology was selected as it provided the best overall combinations performance and was used to perform the time-series analysis.

In complex waters, the CDOM generally increases absorption (reducing reflectance) in the blue bands while the SPM causes an increase in reflectance in the green and red bands due to scattering [53,55]. In the present study, region B of the estuary presented greater contribution of SPM and turbidity and, in fact, the water-leaving reflectance of all AC processors was relatively higher in the green-red bands compared to region A (Figure 4). However, very similar values were observed between the two areas of the estuary in the blue region of the spectra, where also higher values were expected. In conditions of high SPM concentration, the a_CDOM algorithms that use the blue/green-red band ratios tend to underestimate the measured a_CDOM. These factors may have affected the performance of regionally calibrated algorithms (i.e., MA_S, MA_M, CZ, CH, and KU) that, in general, underestimate the a_CDOM retrieval in conditions of higher SPM concentrations. In cases where the a_CDOM does not co-vary with Chl-a, as in the present study, this can significantly affect the bands in the blue to green part of the spectrum due the overlap on the reflected signal [56,57,58,59]. Thus, it is necessary to know the variation of the optically active components present in the water to find alternatives that provide greater accuracy to the retrieval of WQ parameters by satellite remote sensing.

Regarding chlorophyll-a, in the study conducted by Pahlevan et al. [9], the 2-band algorithm originated negative retrievals for low concentrations of Chl-a. The authors considered that this was somehow expected as 2-band algorithms are usually designed to be applied in highly eutrophic waters, a category in which the Sado estuary does not fit. Pereira-Sandoval et al. [53] tested several AC processors to S2-MSI data and observed that, for inland waters with low concentrations of chlorophyll-a, Polymer showed slightly worse results in the NIR bands than C2RCC. Possibly, the existence of low Chl-a values in the Sado estuary, combined with a possible overcorrection performed by Polymer processor, could be the cause of the absence of valid Chl-a values when the GS and 2-band were applied to the data derived from this processor. Interestingly, Chl-a processed by Plymouth Marine Laboratory (PML) using the Calimnos processing chain [60] also did not produce Chl-a matchups with our in situ data for the Gons algorithm [42] coupled with Polymer AC, demonstrating consistency with our results (data not shown). Thus, deriving reliable Chl-a data from S2 is still challenging and further developments are needed. In this context, the on-going H2020 CERTO project (www.certo-project.org (accessed on 26 November 2020)) is likely to provide relevant outputs. This project aims to produce improved satellite data that should be harmonized across different types of complex waters, i.e., from lakes to the sea, with a particular emphasis in transitional waters. Considering its great potential, future work will include testing CERTO’s products, which will be delivered during the project.

Regarding SPM and turbidity retrieval, given the strong relationship between the two parameters (R² = 0.92 in our study area) their relation is worth to be mentioned. Clear differences can be noticed between the turbidity and SPM S2-MSI retrieval, with turbidity presenting overall better correlations with in situ data. This might be caused by different factors. Firstly, turbidity parameter is less subject to in situ measurements errors than SPM. Röttgers et al. [61] demonstrated that significant uncertainties can be associated with most of the SPM measurements. Their results showed that these errors were mainly associated with salt retention and loss of filter material during washing and drying procedures. Moreover, turbidity is an optical parameter highly related to backscattering [45] and thus to reflectance. To explain the poorer retrieval of SPM using a global algorithm with respect to turbidity, the differences between the SPM and turbidity retrieval algorithms should be considered as suggested by [45]. In fact, since turbidity (a measure of side-scattering) is an inherent optical property (IOP), it is not necessary to consider the potential variability of mass specific optical properties, while SPM algorithms do. On the other hand, SPM retrieval algorithms will also be sensitive to particle size, density, and refractive index, which can be important sources of regional variability for retrieval of SPM concentration.

Furthermore, in Figure 6 it is possible to see that the selected algorithms for SPM and turbidity (which use NIR bands) perform better in region B (more turbid waters) compared to region A (clearer waters). Considering the strong spatial gradients, the non-linearity of SPM optics properties in estuarine environment [62] and the interaction with different optically active constituents such as phytoplankton cells, it is key to develop algorithms that are calibrated regionally to accurately quantify SPM from satellite data. Future work will include further investigation on the performance of single-band algorithms considering the different optical properties that characterize the different regions of the estuary. The work will focus on the development of a switching band algorithm that takes into account the different optical properties of region A and region B.

In general, the results of the present study suggest that for satellite-derived a_CDOM, Chl-a and SPM, regional approaches should be adopted for Portuguese transitional waters. Regarding turbidity, the statistical results are encouraging, with reasonably accurate turbidity retrieval using a global algorithm, at least for Region B. As a whole, for the Sado estuary, there is still the need for a more robust database accompanied by in situ and orbital radiometric measurements from remote sensors.

4.2. Spatio-Temporal Analysis

All the analysed parameters showed strong spatial variability within the estuary, as expected. This indicates that to analyze these parameters using in situ measurements, multiple sampling points are needed to fully represent the study area, which is not always feasible. In turn, the seasonal variability of the WQ parameters seemed to be well captured from the S2-MSI satellite products, evidenced in the time series by overlapping the in situ with the satellite-derived data.

For all parameters, higher values were obtained in the interior area of the estuary and during Spring and Summer, in concordance with previous studies [18,63], mainly due to phytoplankton blooms and bottom resuspension caused by stronger wind regimes in the referred period of the year. In the period here studied, the highest SPM values were also found in spring and summer, even if precipitation is stronger in autumn and winter, when higher terrigenous inputs from the rivers are expected. The higher values of SPM (which are historically found in spring) are related to the phytoplankton spring blooms and bottom resuspension due to stronger wind regimes, and not to terrigenous inputs of the SPM [64].The a_CDOM (443) parameter did not correlate or co-vary with the Chl-a, but it had an inverse relationship of R² = −0.57 (p < 0.05) with salinity (not shown here). This suggests that most of the CDOM’s origin is not local (e.g., from the phytoplankton derived components), but rather from bottom resuspension processes due to stronger wind regime within the estuary during spring/summer and from increased rain and riverine input mainly during winter [64,65].

As already mentioned, Chl-a was the parameter that presented the greatest challenges in achieving accurate satellite-derived data. However, it was possible to obtain useful results on the patterns of temporal and spatial variability of the phytoplankton proxy using Sentinel−2 MSI in the Sado estuary. The highest concentrations were found in spring and summer, as expected, even if high values of Chl-a were also observed during winter in the inner region of the estuary. It is also interesting to note that the region of the coastal ocean presented unexpected variation patterns, showing higher values of Chl-a during winter with respect to spring. Furthermore, during summer, the same region did not present any valid Chl-a data, suggesting that the selected algorithms may not be suitable for Chl-a retrieval for the clearer waters in the coastal ocean which present different optical properties, more similar to the open ocean.

Overall, the time-series analysis of all the studied parameters was found to be in agreement with the historical climatological data available for the area and, therefore, the use of the S2-MSI sensor can be considered useful for water quality analysis over the Sado estuary using the selected best-performing processing chain. The high spatial resolution maps provided for the study area are a relevant achievement since no such information was available until the present work, although improvements regarding the bottom effect should be pursued in future work.

4.3. Applicability of Satellite RS Products to WFD

The EU WFD requires that all Member States achieve and/or maintain good ecological status of their water bodies, including estuaries. Under the directive, the assessment of WQ is based on the deviation from a predefined reference condition. However, as discussed by [66], this is one of the most complex, and still under discussion, points of its implementation. For the Sado estuary, in fact, the phytoplankton Chl-a is one of the few water parameters that have a well-defined classification system to assess the ecological status of the water body. According to the range of values observed during the period under analysis, the concentrations never exceeded the reference value for the highest salinity class (6.67 mg m⁻³ for salinities higher than 25 [67,68]) being the observed Chl-a concentration (both in situ and satellite-derived) always <5 mg m⁻³, and therefore indicating that the Sado estuary has currently no clear sign of eutrophication problems. However, it should be reminded that the in situ data considered was always collected at high tide conditions. Higher values might be encountered during low tides.

The lack of historical data for the other physicochemical parameters in the Sado estuary increased the level of difficulty in applying RS WQ products in the context of ecological quality assessment, under the WFD. The definition of reference values is crucial to further explore the potential of RS in the context of the WFD, for which S2-MSI satellites have a great potential. The provided time-series of the WQ parameters obtained through remote sensing observations are crucial to define baseline levels and allowed to confirm the values used to assess the normal variation patterns (i.e., SPM < 20 mgL⁻¹, data not shown). This will help decision-makers to detect potential anomalies and support further studies on ecological aspects. Moreover, these high-resolution maps provide key information to implement effective and meaningful monitoring programs for WQ assessments under the WFD.

Considering the great importance of CDOM as a WQ indicator, we consider that more efforts should be put to retrieve more accurate CDOM absorption values from the S2-MSI data. As discussed previously, a_CDOM is strongly related to salinity [69] and therefore could be used by policymakers as an indicator of fresh-water input within the estuary or to infer drought periods.

One major benefit of using remote sensing data to support the monitoring of surface water status is the capability of providing systematic observations of the whole water-system. The high revisit frequency (5 days) allows an easier detection of small-scale changes and a more synoptic view of the water body. However, a relevant drawback of optical remote sensing data is that cloud coverage and sun glint effects may seriously affect its usefulness [70]. In the present work, only an average of 55% of the total images available (over the 2 years period under investigation) were cloud-free and focused on periods with good atmospheric conditions. This means that image availability is non-homogeneous around the year and can make routine monitoring challenging.

As an example of the applicability of S2-MSI images for water quality assessment, Figure 10A shows the RGB image acquired on 19 February 2020 by the satellite S2A, that captured the dredging activities conducted in the Sado estuary. In the RGB image (top), the sediment plume produced by the dredge is clearly visible (red box). Using the selected algorithm for SPM retrieval (Figure 10B), we observed that a value of 35 mgL⁻¹ of suspended solids was reached in the pixels corresponding to the sediment plume, compared to an average of 2 mgL⁻¹ of SPM in the surrounding pixels. Considering that concentration of SPM in the same area never exceeded 8 mgL⁻¹ during the studied period, it can be concluded that, although locally, the dredging activities performed in the area had a severe impact on the water clarity of that specific area of the estuary. Unfortunately, it was not possible to follow the evolution of impacts during the dredging period due to cloud-coverage of Sentinel-2 images.

Finally, we consider that the approach followed in this research may be extended to other Portuguese transitional systems and European countries in support of the WFD considering the low costs of the data production with respect to continuous in situ monitoring programs.

5. Conclusions

Over the past few decades, satellite remote sensing techniques have increasingly been used to support environmental monitoring and to assess different types of water bodies, mainly because of their cost-effectiveness and spatio-temporal resolution [71]. The aim of the present study was to demonstrate the capabilities and knowledge gaps of S2-MSI data for continuous monitoring of WQ parameters in a highly dynamic system and evaluate its usefulness under the WFD.

In fact, the Sentinel-2 satellites showed their encouraging potential for routine monitoring of WQ parameters in the Sado estuary. Their medium temporal resolution (5 days at mid latitudes) and, therefore, the higher frequency of satellite data regarding in situ observations, lead to an easier detection of small-scale changes and allow a synoptic view and analysis of the system.

In the context of the WFD, this study demonstrates the high potential of S2-MSI in water quality mapping and monitoring. However, a regional validation and calibration of the bio-optical algorithms is suggested. The turbidity product showed high consistency with in situ observations, proving the great capabilities of the MSI sensor to monitor this important water quality parameter. However, for the key parameter Chl-a further research is needed, with a more complete set of match-ups. Moreover, in the context of the WFD, the need to define reference values for variables such a_CDOM, SPM, and turbidity for the Sado estuary (and other Portuguese estuaries) may be overcome through the provided remote sensing data (maps and time series) which are fundamental to be used as a baseline information for the characterization of the normal conditions of the Sado estuary.

Improvements in the management of water resources can be achieved only through the delivery of continuous and sustained WQ data, to allow forecasts at local and regional scales and achieve an effective decision-making process. In this context, remote sensing is an extraordinary tool that allows a better understanding of how water bodies are changing due to environmental perturbations, anthropogenic pressure and climate change.

Future work will focus on a more direct link between S2-MSI derived products and WFD for the determination of the ecological status of Portuguese transitional systems. In addition, to further evaluate the best processing chain for estuarine systems as the Sado estuary, the collection of in situ radiometric data is key, in order to carry out robust satellite data validation and calibration exercises. It would also be interesting to use remote satellite detection to monitor the water quality during the course of specific engineering operations, such as dredging activities, that, as observed, may cause relevant changes in the estuarine system.

The present work was the first exercise to provide inter-annual and seasonal variability of several WQ parameters, with high spatial resolution (10 m) for the Sado estuary, being fundamental for the establishment of effective monitoring in the study area.