On the Potential of Robust Satellite Techniques Approach for SPM Monitoring in Coastal Waters : Implementation and Application over the Basilicata Ionian Coastal Waters Using MODIS-Aqua

Monitoring river plume dynamics and variations in complex coastal areas can provide useful information to prevent marine environmental damage. In this work, the Robust Satellite Techniques (RST) approach has been implemented and tested on historical series of Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) data to monitor, for the first time, Suspended Particulate Matter (SPM) anomalies associated to river plumes. To this aim, MODIS-Aqua Level 1A data were processed using an atmospheric correction adequate for coastal waters, and SPM daily maps were generated applying an algorithm adapted from literature. The RST approach was then applied to these maps to assess the anomalous presence of SPM. The study area involves the Basilicata region coastal waters (Ionian Sea, South of Italy). A long-time analysis (2003–2015) conducted for the month of December allows us to find that the maximum SPM concentration value was registered in December 2013, when an extreme hydrological event occurred. A short-time analysis was then carried out applying RST to monitor the dynamics of anomalous SPM concentrations. Finally, the most exposed areas, in terms of SPM concentration, were identified. The results obtained in this work showed the RST high potential when used in combination with standard SPM daily maps to better characterize and monitor coastal waters.


Introduction
Variations in Suspended Particulate Matter concentration (SPMc) due to natural (e.g., flood events) or man-made (e.g., dam building or marine dredging activities) causes can impact the coastal water quality and, therefore, marine ecosystems [1,2].The European Water Framework Directive (WFD) [3] establishes a framework for the protection of the transition waters, such as coastal waters and estuaries, implementing the monitoring activities of phytoplankton and benthic algae [4].Coastal areas, especially in presence of river mouths, are subject to large SPM variations in time and space [5], which can be considered as a stress factor for many pelagic and benthic organisms [5][6][7].
An increase in SPMc, which can be also produced by wind-wave resuspension in shallow waters [8], determines water turbidity variation with a direct effect on warming processes, influencing dissolved oxygen concentration [9].The SPM distribution varies in different coastal areas and over a large spatiotemporal domain, hence making the SPM monitoring unprofitable by using traditional field methods, which can provide accurate measurements only at a very small spatiotemporal scale.Remote sensing techniques have become a valuable alternative to in situ measurements [10], mainly using data acquired in the red and near-infrared (NIR) spectral regions of the electromagnetic spectrum, where SPM usually dominates the seawater spectral signature [10,11].In particular, [11] found that the wavelength associated to the maximum in remote sensing reflectance (R rs , normalized water-leaving radiance/mean extraterrestrial solar irradiance, (sr −1 )) increases from green to red up to NIR wavelengths with increasing SPMc.Focusing on Gironde estuary, [11] found that R rs grows between 400 and 700 nm for moderately high SPM values (35-250 g/m 3 ), while R rs tends to saturate at these wavelengths for extremely turbid waters (SPMc higher than 250 g/m 3 ) but still considerably increases in the spectral NIR region.
Many satellite radiometers, capable of acquiring information in the visible and NIR spectral regions, have been exploited so far for such a kind of application.Among them, MODIS is the one operational at a global scale with a good trade-off between spatial and temporal resolution.It is on-board Earth Observing System (EOS) Terra (since 2000) and Aqua (since 2002) satellites and has two spectral bands in the red (band 1, 620-670 nm) and NIR band 2, 841-876 nm) regions at a spatial resolution of 250 m, acquired twice per day, almost all over the world.Several authors have used MODIS single band algorithms [1,9,[12][13][14][15] as well as the combination of bands [16][17][18][19][20] for retrieving SPMc in different geographical areas characterized by different SPM features.Even if the proposed algorithms perform well for the specific region where they were calibrated, their accuracy is still poorly demonstrated over different geographic regions [21].Some efforts have been made in the past years to develop a more general algorithm valid for a wide range of geographical and geophysical settings as well as SPM concentrations.[10] proposed a semi analytical single band algorithm for SPMc retrieving, for different sensors and wavelengths, calibrated using in situ data (SPM range: from ~1 to ~110 g/m 3 ) from the Southern North Sea, finding an error of about 30% and 40% in calibration and validation respectively.More recently, [21] tested different algorithms (including [10]) for retrieving SPMc on a broad range of concentrations (from 0.15 to 2626 g/m 3 ) and for various coastal environments dominated by river discharge, resuspension or phytoplankton bloom in Europe, French Guiana, North Canada, Vietnam, and China.Among the algorithms tested, the one proposed by [10] resulted the best, assuming that the model parameters were adapted using the in situ measurements.These preliminary studies confirmed the uncertainties of such algorithms if applied to different SPMc ranges or geophysical characteristics, suggesting the need of further analyses for a better accuracy assessment.Another source of uncertainty for these algorithms concerns the implemented atmospheric correction procedures.Even if several atmospheric correction algorithms have been proposed for turbid coastal waters (e.g., [22][23][24]), their performance changes on the basis of the analyzed location and range of SPMc.Moreover, the SPMc values and variations are usually unknown for those areas where no previous studies are available.Therefore, considering the high SPM variability worldwide, it is quite difficult to understand the significance of an absolute SPMc value.The same difficulty, considering the SPM variability in the spatiotemporal domain, can occur also in the same study area.For instance a SPMc of 50 g/m 3 in an area normally characterized by high SPMc values (e.g., estuarine or coastal area under the influence of very turbid rivers, etc.) could have a low specific impact on the local environment.Instead, the same SPMc value, if registered in areas normally characterized by lower SPM concentrations, might have stronger effects and lead to serious environmental risks.
In this work, we investigate the potential of the Robust Satellite Techniques (RST) [25][26][27] approach to monitor SPMc variations in coastal area.Taking into account the signal historical behavior at local (i.e., pixel) level, RST allows the characterization of the unperturbed conditions (normally registered) of the analyzed areas independently from each specific site-setting.Once such unperturbed conditions are known, statistically anomalous SPMc variations can be detected, allowing us to infer their dynamics and spatiotemporal evolution.The Basilicata region (South Italy) coastal area, poorly investigated by ocean color remote sensing and never for SPM, has been selected as study area.Implementing RST on SPMc daily maps, anomalies associated with several river plumes in the short-and long-term periods are identified and analyzed.

Study Area
The Basilicata Region coastal waters (BRCW, Figure 1), located in the Ionian Sea (Mediterranean Sea), have been investigated in this work.The Basilicata region is situated in Southern Italy and its main rivers, all flowing into the Ionian Sea are, from N-E to S-W, the Bradano, Basento, Cavone, Agri, and Sinni.
Remote 2016, 8, 922 3 of 20 coastal area, poorly investigated by ocean color remote sensing and never for SPM, has been selected as study area.Implementing RST on SPMc daily maps, anomalies associated with several river plumes in the short-and long-term periods are identified and analyzed.

Study Area
The Basilicata Region coastal waters (BRCW, Figure 1), located in the Ionian Sea (Mediterranean Sea), have been investigated in this work.The Basilicata region is situated in Southern Italy and its main rivers, all flowing into the Ionian Sea are, from N-E to S-W, the Bradano, Basento, Cavone, Agri, and Sinni.The Bradano, Basento and Cavone river basins are characterized by low levels of precipitations and a small number of springs, implying lower flow rates during summer when compared to the Agri and Sinni basins characterized by a higher precipitation rate and a larger number of springs [28].The annual minimum flow rate of the Sinni, Agri, Bradano and Basento is 1.38, 0.5, 0.04, and 0.08 m 3 /s, respectively, whereas for the Cavone it is practically nil [28].In situ measurements of maximum flow rate are not available for these rivers.The maximum flow associated to a return-period of 30 years, computed according to [29], is 2635, 1408, 637, 2073 and 1282 m 3 /s for the Bradano, Basento, Cavone, Agri and Sinni, respectively [30].
Very few studies have been made using ocean color satellite data over BRCW (e.g., [31]), and none of them refers to SPM monitoring.Studies using in situ measured data to analyze the SPM dynamics in the BRCW are also very scarce.Among them, Buccolieri et al. [32] investigated the heavy metals distribution and accumulation in Taranto Gulf marine sediments, finding that Basilicata rivers are an important source of heavy metals for this area.Buccolieri et al. [32] speculated that the main reason of their accumulation in the East part of the gulf is due to the combined effects of the general anticlockwise marine currents and the river discharges (as source of contaminants) on the Basilicata region coast.
In recent years, several flood events affected the Basilicata region coastal area [33] due to the increase in the frequency and intensity of moderate-to-heavy rainfalls [34].Figure 2a shows the average daily river water levels at the five gauging stations closest to the river mouths (black dots in The Bradano, Basento and Cavone river basins are characterized by low levels of precipitations and a small number of springs, implying lower flow rates during summer when compared to the Agri and Sinni basins characterized by a higher precipitation rate and a larger number of springs [28].The annual minimum flow rate of the Sinni, Agri, Bradano and Basento is 1.38, 0.5, 0.04, and 0.08 m 3 /s, respectively, whereas for the Cavone it is practically nil [28].In situ measurements of maximum flow rate are not available for these rivers.The maximum flow associated to a return-period of 30 years, computed according to [29], is 2635, 1408, 637, 2073 and 1282 m 3 /s for the Bradano, Basento, Cavone, Agri and Sinni, respectively [30]. Very few studies have been made using ocean color satellite data over BRCW (e.g., [31]), and none of them refers to SPM monitoring.Studies using in situ measured data to analyze the SPM dynamics in the BRCW are also very scarce.Among them, Buccolieri et al. [32] investigated the heavy metals distribution and accumulation in Taranto Gulf marine sediments, finding that Basilicata rivers are an important source of heavy metals for this area.Buccolieri et al. [32] speculated that the main reason of their accumulation in the East part of the gulf is due to the combined effects of the general anticlockwise marine currents and the river discharges (as source of contaminants) on the Basilicata region coast.
In recent years, several flood events affected the Basilicata region coastal area [33] due to the increase in the frequency and intensity of moderate-to-heavy rainfalls [34].Figure 2a shows the average daily river water levels at the five gauging stations closest to the river mouths (black dots in Figure 1) for the 1 January 2011-30 June 2015 period (data were collected by the Civil Protection Department of Basilicata Region).All rivers show a strong seasonal pattern, characterized by maximum water levels during the autumn-winter period.Among the five rivers, the Sinni shows the lowest variability due to the presence of the Monte Cotugno dam on its path, which is able to control the river flow variation in case of extreme events (floods and droughts).During the almost five years of available river level data the maximum values were registered on 1 December 2013 for the Sinni (3.6 m) and Agri (4.19 m), on 2 December 2013 for the Basento (7.92 m) and Cavone (5.22 m).The two maximum water levels for the Bradano were registered on 2 March 2011 (6.04 m) and 2 December 2013 (5.39 m).These data suggested focusing on the extreme event occurred in the Basilicata region between 30 November and 7 December 2013, caused by a significant amount of rainfall due to the "Ciclone Nettuno" storm.A daily value of 130 mm was recorded on 1 December 2013 and a cumulative one of about 200 mm for the entire event [33] was found at the Bradano river gauging station.
Figure 1) for the 1 January 2011-30 June 2015 period (data were collected by the Civil Protection Department of Basilicata Region).All rivers show a strong seasonal pattern, characterized by maximum water levels during the autumn-winter period.Among the five rivers, the Sinni shows the lowest variability due to the presence of the Monte Cotugno dam on its path, which is able to control the river flow variation in case of extreme events (floods and droughts).During the almost five years of available river level data the maximum values were registered on 1 December 2013 for the Sinni (3.6 m) and Agri (4.19 m), on 2 December 2013 for the Basento (7.92 m) and Cavone (5.22 m).The two maximum water levels for the Bradano were registered on 2 March 2011 (6.04 m) and 2 December 2013 (5.39 m).These data suggested focusing on the extreme event occurred in the Basilicata region between 30 November and 7 December 2013, caused by a significant amount of rainfall due to the "Ciclone Nettuno" storm.A daily value of 130 mm was recorded on 1 December 2013 and a cumulative one of about 200 mm for the entire event [33] was found at the Bradano river gauging station.An important contribution to SPMc variations may also be caused by wind-wave resuspension in shallow waters [8].To better evaluate such an effect, wind data acquired by the Italian agency for environmental Protection and Research (ISPRA, [35]) were analyzed.In detail, considering that no wind station is available within the study area, we used the information acquired at the Taranto one, located just outside the study area (less than 2 km) in the North-East sector.The registered wind intensity and direction are shown in Figure 3 by the length and the direction of the green bars respectively, measured every 10 min in the 1-13 December 2013 period.An important contribution to SPMc variations may also be caused by wind-wave resuspension in shallow waters [8].To better evaluate such an effect, wind data acquired by the Italian agency for environmental Protection and Research (ISPRA, [35]) were analyzed.In detail, considering that no wind station is available within the study area, we used the information acquired at the Taranto one, located just outside the study area (less than 2 km) in the North-East sector.The registered wind intensity and direction are shown in Figure 3 by the length and the direction of the green bars respectively, measured every 10 min in the 1-13 December 2013 period.During the investigated period the highest wind speed value, equal to 13.9 m/s, was registered on 1 December 2013, while the averaged wind speed values measured during the event, i.e., on 4, 5 and 6 of December (analyzed in detail in Section 4) were 3.5 m/s, 1.8 m/s and 1.7 m/s, respectively.The value recorded on 13 December, i.e., on the first image used later for the falsification analysis (see Section 4.4), was 3.4 m/s.The absence of systematic in situ observations as well as of previous works on the relationship between SPM and wind in the studied area does not allow assessing whether these values can produce a sediment resuspension.Nevertheless, the comparable wind intensity recorded both during the days of the event (i.e., 4, 5 and 6 of December 2013) and after its end (i.e., 13 December 2013), seems to indicate that such intensities are not able to produce a detectable resuspension effect.Moreover, considering, as an example, that [36] found for their specific area of study (Barataria, SouthEast Louisiana, USA) that the resuspension effect was predominant over other SPM sources for wind speed higher than about 4 m/s, we can further speculate that the main contribution of SPM in BRCW for the specific investigated period is river input.A more general discussion about the relationship among SPM, river level and wind speed will be presented in Section 4.3.

Satellite Data Processing for SPM Daily Maps Retrieval
Data recorded by the MODIS sensor on board Aqua satellite have been analyzed in this study.Due to the presence of turbid waters in BRCW, the NASA's standard ocean color Level 2 (L2) product, obtained applying the NIR atmospheric correction [22,37], systematically presents a lot of negative Rrs values in the red and NIR spectral bands, especially near the coast, thus making these data unusable.For this reason, MODIS Level 1A data, directly downloaded from NASA's ocean color web site [38], were reprocessed using the SeaDAS software (version 7.3).After producing Level 1B and geolocation files, two atmospheric correction methods, NIR-SWIR [24] and MUMM [23], properly designed for coastal waters, were evaluated.The NIR-SWIR approach first computes a turbidity During the investigated period the highest wind speed value, equal to 13.9 m/s, was registered on 1 December 2013, while the averaged wind speed values measured during the event, i.e., on 4, 5 and 6 of December (analyzed in detail in Section 4) were 3.5 m/s, 1.8 m/s and 1.7 m/s, respectively.The value recorded on 13 December, i.e., on the first image used later for the falsification analysis (see Section 4.4), was 3.4 m/s.The absence of systematic in situ observations as well as of previous works on the relationship between SPM and wind in the studied area does not allow assessing whether these values can produce a sediment resuspension.Nevertheless, the comparable wind intensity recorded both during the days of the event (i.e., 4, 5 and 6 of December 2013) and after its end (i.e., 13 December 2013), seems to indicate that such intensities are not able to produce a detectable resuspension effect.Moreover, considering, as an example, that [36] found for their specific area of study (Barataria, SouthEast Louisiana, USA) that the resuspension effect was predominant over other SPM sources for wind speed higher than about 4 m/s, we can further speculate that the main contribution of SPM in BRCW for the specific investigated period is river input.A more general discussion about the relationship among SPM, river level and wind speed will be presented in Section 4.3.

Satellite Data Processing for SPM Daily Maps Retrieval
Data recorded by the MODIS sensor on board Aqua satellite have been analyzed in this study.Due to the presence of turbid waters in BRCW, the NASA's standard ocean color Level 2 (L2) product, obtained applying the NIR atmospheric correction [22,37], systematically presents a lot of negative R rs values in the red and NIR spectral bands, especially near the coast, thus making these data unusable.
For this reason, MODIS Level 1A data, directly downloaded from NASA's ocean color web site [38], were reprocessed using the SeaDAS software (version 7.3).After producing Level 1B and geolocation files, two atmospheric correction methods, NIR-SWIR [24] and MUMM [23], properly designed for coastal waters, were evaluated.The NIR-SWIR approach first computes a turbidity index used for selecting which bands, between the NIR or SWIR ones, have to be used for the atmospheric correction (aerosol contribution).The NIR approach is used for waters with a lower level of turbidity, while the SWIR one is used in case of highly turbid waters, as the seawater reflectance in the SWIR region is negligible independently from the SPMc value [39].Briefly, for waters where the turbidity index is lower than a certain value (i.e., 1.3) the standard NIR atmospheric correction is applied [24], while for pixels with a higher value the SWIR bands are used for the atmospheric correction.The MUMM atmospheric correction is, instead, based on the assumption of a constant ratio of both the aerosols and water leaving reflectances measured at 765 and 856 nm [23], the so-called similarity spectrum [40].This assumption, initially formulated for SeaWiFS, has been exported on MODIS sensor assuming the homogeneity of the 748:869 nm ratio for aerosol and water leaving reflectance.The MODIS 748 (band 15) and 869 (band 16) bands are both at 1000 m of spatial resolution.These two bands were automatically downscaled at 250 m spatial resolution by SeaDAS during the atmospheric correction process in order to apply such a correction to the two MODIS bands used for SPM computation in this work (band 1, 645 nm and band 2, 859 nm) at 250 m spatial resolution.
L2 satellite products obtained applying [23,24] for the subset shown in Figure 1b (Lat/Long projection, WGS84) were found to be very similar for most of the images analyzed over BRCW.However, the NIR-SWIR atmospheric correction returns a noisier product when compared to the MUMM one, producing several negative R rs values in the waters with a low level of turbidity.Therefore, the latter was chosen as the most suitable for atmospheric corrections in BRCW.As a first result, considering the validity range of the MUMM atmospheric correction method [41], we can speculate that BRCW is characterized by low-to-moderate SPMc values.
Subsequently we focused on the SPMc retrieval.Considering that no in situ measurements were available to define a local algorithm, we used a standard one.Nechad et al. [10] proposed a general semi-analytical algorithm for retrieving SPMc in turbid coastal waters and validated it within the 1-100 g/m 3 concentration range in Southern North Sea.The general formulation of this single band algorithm can be written as follows: where SPMc is the SPM concentration in g/m 3 , ρ w (λ) = π × R rs (λ) is the water leaving reflectance, A p (λ) (in g/m 3 ) and C p (λ) (dimensionless) are the algorithm coefficients, which depend on the spectral band considered.The use of MODIS-250 m red band (645 nm) reflectance in Equation (1) produces, especially when high river water levels occur, SPMc values higher than 200 g/m 3 in BRCW.These high values are in contrast with those obtained by previous studies [31,42], where the use of red band proved to be suitable for the retrieval of SPMc up to ~50 g/m 3 , while the use of the NIR band is more appropriate for higher SPMc values.In order to avoid such unphysical SPM concentrations achieved using only the red band, a combined approach, as previously proposed by [21], can be implemented.Red band R rs can be used for SPMc <~50 g/m 3 , while the NIR one is preferable for higher SPMc values.In detail, using the coefficients given by [10] for the first two MODIS bands, Equation (1) can be written as follows: The coefficients A p (λ) and C p (λ) of Equation ( 1), were specifically provided by [10] for λ = 645 nm in their Table 1.Whereas those for λ = 859 nm, not furnished, were instead obtained from a linear interpolation of the ones provided for λ = 857.5 nm and λ = 860 nm, respectively.Equation ( 2) is used for the slightly-to-moderately turbid waters (R rs (645) ≤ 0.03 sr −1 , corresponding to SPMc < ~50 g/m 3 ), while Equation ( 3) is applied for more turbid waters (R rs (645) ≥ 0.04 sr −1 , SPMc > ~150 g/m 3 ).SPMc corresponding to the 0.03 < R rs < 0.04 sr −1 range is computed using a weighted average, calculating the weights by applying a logarithmic function (see [21] for more details on the methods).Keeping in mind that [10] found, as mentioned above, an error of about 30% and 40% in calibration and validation respectively, it can be speculated that such an error affects at least the absolute SPMc maps computed over BRCW.
A summary of the pre-processing steps implemented to obtain the SPMc daily maps starting from the MODIS Level 1A data is shown in Figure 4a.Equation ( 2) is used for the slightly-to-moderately turbid waters (Rrs(645) ≤ 0.03 sr −1 , corresponding to SPMc < ~50 g/m³), while Equation ( 3) is applied for more turbid waters (Rrs(645) ≥ 0.04 sr −1 , SPMc > ~150 g/m³).SPMc corresponding to the 0.03 < Rrs < 0.04 sr −1 range is computed using a weighted average, calculating the weights by applying a logarithmic function (see [21] for more details on the methods).Keeping in mind that [10] found, as mentioned above, an error of about 30% and 40% in calibration and validation respectively, it can be speculated that such an error affects at least the absolute SPMc maps computed over BRCW.
A summary of the pre-processing steps implemented to obtain the SPMc daily maps starting from the MODIS Level 1A data is shown in Figure 4a.

The RST Approach for Identification of SPM Spatiotemporal Anomalies
RST is a general multi-temporal approach of data analysis, successfully implemented for investigating different natural and anthropic risks.A detailed description of the approach and a list of its wide range of applications can be found in [27], while in [43] its recent application to identify areas of near-surface chlorophyll-a anomalies in oligotrophic waters is described.Figure 4b shows the RST implementation steps for the specific purpose of this study.The approach starts collecting multi-year homogeneous (e.g., same investigated signal/parameter, same location, same month and acquisition time) series of satellite records.Then, the expected value and natural variability of the investigated signal are computed.Finally, a signal anomaly is statistically identified when the signal under investigation exceeds such an expected value taking also into account its natural variability.To this aim, a local variation index, named ALICE (Absolutely Local Index of Change of the Environment, namely ⊗V(x, y, t)), is computed [25][26][27] and expressed in its general formula as:

The RST Approach for Identification of SPM Spatiotemporal Anomalies
RST is a general multi-temporal approach of data analysis, successfully implemented for investigating different natural and anthropic risks.A detailed description of the approach and a list of its wide range of applications can be found in [27], while in [43] its recent application to identify areas of near-surface chlorophyll-a anomalies in oligotrophic waters is described.Figure 4b shows the RST implementation steps for the specific purpose of this study.The approach starts collecting multi-year homogeneous (e.g., same investigated signal/parameter, same location, same month and acquisition time) series of satellite records.Then, the expected value and natural variability of the investigated signal are computed.Finally, a signal anomaly is statistically identified when the signal under investigation exceeds such an expected value taking also into account its natural variability.To this aim, a local variation index, named ALICE (Absolutely Local Index of Change of the Environment, namely ⊗ V (x, y, t)), is computed [25][26][27] and expressed in its general formula as: In this equation V(x, y, t) is the investigated signal acquired at location (x, y) and at time t, selected according to the type of phenomenon to be studied.It can correspond to the signal measured in a single channel, a combination of spectral bands or another physical parameter (e.g., surface temperature, SPM or Chlorophyll-a concentration, etc.).V ref (x, y) and σ ref (x, y) represent the value of the investigated signal expected in unperturbed conditions and its natural variability named the "reference fields".The two reference fields are computed considering only cloud-free records belonging to the collected satellite imagery, acquired under similar observational conditions.V ref is usually expressed as the temporal mean µ V (x, y), but, for specific applications, it can be the minimum min V (x, y) or maximum max V (x, y) value historically recorded [25].σ ref (x, y) is usually expressed as the standard deviation.Considering the scope of the present study, the index in Equation ( 4) can be written as follows: where SPMc(x, y, t) is the SPM concentration at location (x, y) and time t, µ SPM (x, y) and σ SPM (x, y) are the interannual reference fields representing the unperturbed conditions and the normal variability at pixel level.Both reference fields were computed using time series of the daily MODIS-SPMc maps for the month of December for the years 2003-2015, obtained applying the approach suggested in the previous steps (Section 3.1, Figure 4a).For each pixel of the study area an iterative kσ procedure [25] was implemented to exclude outliers, to make the reference fields actually representative of unperturbed conditions.About 500 MODIS-Aqua images recorded in December were collected, including those only partially covering the BRCW, to produce the reference fields, which are shown in Figure 5a,b.While processing the images, for each of them, the pixels affected by clouds, sun glint, or large satellite zenith angle were discarded from further analysis.Figure 5c shows that for many pixels in the area more than 100 records (survived to the above mentioned discarding procedure) were used for computing µ SPM (x, y) and σ SPM (x, y), with a median value of ~92 useful record per pixel over the whole area.Such a number has been found to be, by an independent study applied on RST approach [44], high enough to be representative of the investigated long-term time-series signal.The identification of cloudy pixels has been performed using data acquired in the SWIR bands [45].After testing various thresholds, we used a value of 0.012 on the Rayleigh-corrected reflectance (ρ c = ρ TOA − ρ r ) 2130 nm SWIR band for detecting cloudy pixels.It is worth mentioning that a smaller number of images was used for computing the reference fields close to the coastline, compared to the other pixels of the study area.The sub-pixel effect due to the presence of a water-land mixed signal might produce a failure of atmospheric correction, especially for off-nadir acquisitions, making the SPM retrieval impossible in these areas.
Concerning the reference fields shown in Figure 5a,b, as expected, the highest SPMc values and variability were found near the coast, particularly close to the river mouths.
The index of Equation ( 5) should exhibit high positive values when SPMc is considerably higher than normal conditions, also considering the normal variability of the measured signal for the specific place and period of the year.The normal variability allows taking into account all the possible sources of noise (random or not, known or not) that can modify the measured signal at sea surface in coastal areas, such as river plumes, tidal effects, adjacency effects, re-projection issues and atmospheric correction problems.The possible variability in SPM absolute values due to computation uncertainties also provides an increase of the normal variability as well as a signal trend, while potential errors related to the use of an algorithm calibrated for another specific region are minimized by the inherent differential nature of RST approach.
Considering that the ALICE index is a Gaussian standardized variable (with µ~0 and σ~1, see Figure 5d), ⊗ SPM (x, y, t) increasing values can be associated to statistically anomalous events.⊗ SPM (x, y, t) values higher than 2, 3 and 5, indeed, indicate a decreasing occurrence probability equal to 2.27%, 0.15% and 2.8 × 10 −5 %, respectively, while higher ⊗ SPM (x, y, t) values are associated to an even lower occurrence probability.The index of Equation ( 5) should exhibit high positive values when SPMc is considerably higher than normal conditions, also considering the normal variability of the measured signal for the specific place and period of the year.The normal variability allows taking into account all the possible sources of noise (random or not, known or not) that can modify the measured signal at sea surface in coastal areas, such as river plumes, tidal effects, adjacency effects, re-projection issues and atmospheric correction problems.The possible variability in SPM absolute values due to computation uncertainties also provides an increase of the normal variability as well as a signal trend, while potential errors related to the use of an algorithm calibrated for another specific region are minimized by the inherent differential nature of RST approach.
Considering that the ALICE index is a Gaussian standardized variable (with μ~0 and σ~1, see Figure 5d), ⊗SPM(x, y, t) increasing values can be associated to statistically anomalous events.⊗SPM(x, y, t) values higher than 2, 3 and 5, indeed, indicate a decreasing occurrence probability equal to 2.27%, 0.15% and 2.8 × 10 −5 %, respectively, while higher ⊗SPM(x, y, t) values are associated to an even lower occurrence probability.

Interannual Analysis
Before applying RST to the event occurred in December 2013, the SPM interannual variability was studied for this month during the 2003-2015 period (Figure 6).To this aim, first the reference field outputs of the interannual analysis previously described (Figure 5a,b) were spatially averaged over BRCW.The spatial averaged values were found to be 1.70 g/m 3 and 0.97 g/m³, respectively, for

Interannual Analysis
Before applying RST to the event occurred in December 2013, the SPM interannual variability was studied for this month during the 2003-2015 period (Figure 6).To this aim, first the reference field outputs of the interannual analysis previously described (Figure 5a,b) were spatially averaged over BRCW.The spatial averaged values were found to be 1.70 g/m 3 and 0.97 g/m 3 , respectively, for the interannual mean (µ SPM ) and standard deviation (σ SPM ).µ SPM is represented by the red dashed line in Figure 6, while the dashed green one indicates the normal SPMc variability interval (i.e., µ SPM ± 1 × σ SPM = 2.67 g/m 3 ).For each December from 2003 to 2015, the mean SPMc as well as its standard deviation were also computed at monthly (M) scale.Both these values were spatially averaged over BRCW, obtaining the blue rhombus (µ SPM−M ) and the error bars (µ SPM−M ± 1 × σ SPM−M ), shown in Figure 6.
Three µ SPM−M main relative peaks have been identified in the 2003-2015 period: December 2004 (2.67 g/m 3 ), December 2008 (2.47 g/m 3 ) and December 2013 (2.75 g/m 3 ).Among these only the one corresponding to December 2013 shows a value higher than the normal variability interval (December 2004 is still within this range), indicating an anomalous behavior of this year with respect to the mean SPMc interannual value.Furthermore, for this year also the highest SPMc standard deviation variability (i.e., the error bar) was detected, suggesting strong spatial and temporal SPM dynamics.For the other years the behavior was similar and the higher µ SPM−M is, the higher σ SPM−M is.Indeed, years characterized by lower relative µ SPM−M values, such as 2010, 2012, 2014 and 2015, are also characterized by a lower σ SPM−M .A decreasing weak SPMc trend (with p-value = 0.025) has also been recognized (dashed black line in Figure 6), which, in any case has not significantly influenced signal variability and hence RST sensitivity.The standard deviation of trended and detrended series shown in Figure 5, indeed, is very similar (0.55 g/m 3 for the first and 0.49 g/m 3 for the latter).Such a trend can be associated to a reduction in the precipitation rate during the autumn-winter season in the Mediterranean basin as observed for the 1951-2010 period for the Bradano river catchment [46].Despite this general decreasing tendency, the occurrence of moderate-to-heavy rainfalls has increased in frequency and intensity over the last decade [34], generating for example the extreme December 2013 event.A climatological analysis which takes into account all the calendar months should be carried out in the future to better assess such behavior.
corresponding to December 2013 shows a value higher than the normal variability interval (December 2004 is still within this range), indicating an anomalous behavior of this year with respect to the mean SPMc interannual value.Furthermore, for this year also the highest SPMc standard deviation variability (i.e., the error bar) was detected, suggesting strong spatial and temporal SPM dynamics.For the other years the behavior was similar and the higher ̅ is, the higher is.Indeed, years characterized by lower relative μ values, such as 2010, 2012, 2014 and 2015, are also characterized by a lower .A decreasing weak SPMc trend (with p-value = 0.025) has also been recognized (dashed black line in Figure 6), which, in any case has not significantly influenced signal variability and hence RST sensitivity.The standard deviation of trended and detrended series shown in Figure 5, indeed, is very similar (0.55 g/m 3 for the first and 0.49 g/m 3 for the latter).Such a trend can be associated to a reduction in the precipitation rate during the autumnwinter season in the Mediterranean basin as observed for the 1951-2010 period for the Bradano river catchment [46].Despite this general decreasing tendency, the occurrence of moderate-to-heavy rainfalls has increased in frequency and intensity over the last decade [34], generating for example the extreme December 2013 event.A climatological analysis which takes into account all the calendar months should be carried out in the future to better assess such behavior.

Daily Analysis for the Event of December 2013
Here we analyze the results obtained for the extreme hydrological event affected the Basilicata region in December 2013.SPMc maps were computed by applying the combined approach described in Section 3.1, then the RST-based Equation ( 5) was applied to these data.It is worth noting that the first MODIS-Aqua cloud-free image over the BRCW related to the event was acquired on 4 December 2013 at 12:25 GMT, therefore with a three-day time lag since the maximum water level was registered at the Agri and Sinni stations, and a two-day time lag since that measured at the Basento, Bradano and Cavone (see Section 2 and Figure 2b).
Figure 7 shows the SPM maps obtained applying the method described in Section 3.1 to the MODIS-Aqua images respectively acquired on 4, 5 and 6 December 2013.In these maps the color bar describes the SPMc values, while in grey are depicted those pixels where no SPM retrieval was possible due to no-data value (failure of atmospheric correction due to sub-pixel effects close to the coast) or cloudy pixels.On 5 December, when the general currents have apparently changed, the area indicated by the green arrow in Figure 8a is disconnected from the other anomalous pixels, which move, in the case of the Sinni and Agri river plumes, to the North-East direction.On 6 December only very few anomalous pixels are still visible, extending in a rectilinear shape for more than 10 km from the coastline.
Although Figures 7 and 8 highlight similar SPMc features, the results reported in these two figures are supplementary.Figure 7a shows where high absolute SPMc values are located, but without any knowledge about the history of the area at pixel level, it is difficult to attribute a significance to these SPMc values.On the other hand, RST allows us to identify areas with SPMc values significantly higher than the value observed in normal conditions, taking also into account the signal variability (standard deviation).Therefore, the same SPMc absolute value can generate very different ⊗SPM(x, y, t) magnitudes.To better understand the differences between the anomalies obtained computing ⊗SPM(x, y, t) and the SPMc values, the signal profiles along the A-B transect in Figures 7c and 8c are reported in Figure 9.In detail, in this figure, the ⊗SPM(x, y, t) variation is depicted in red and the SPMc one in blue.The purple line provides the value of SPM temporal mean plus two times the standard deviation (i.e., μSPM(x, y) + 2 × σSPM(x, y)) for each pixel of the profile investigated.Moreover, the red dashed line indicates the ⊗SPM(x, y, t) threshold of 2. In the first part of the transect The plumes associated to the five main rivers in the study area are clearly visible, with SPMc values higher than 100 g/m 3 , especially very close to the river mouths on 4 and 5 December (Figure 7a,b).The plume characterized by the greatest extension and the highest SPMc values is that of the Basento river, while the smallest one, both in terms of extension and values, is the Agri plume.These results are well in agreement with the river water levels measured (Figure 2b).Comparing the images acquired in the above-mentioned three days, it is possible to confirm the very high SPM spatial and temporal dynamics, at least in the upper water layers.In fact, on 6 December, Figure 7c shows SPM values lower than 20 g/m 3 in almost all the BRCW, with the highest SPMc values in correspondence of the Bradano and Basento mouths.The distance from the coastline where significant SPMc values were found differs for each river.In fact, the plumes associated to the Sinni and Agri extend over a distance up to 20 km, while the Cavone, Basento and Bradano plumes cover a distance up to 10 km on 5 December (Figure 7b).On 4 December all river plumes extend from North-East to South-West (Figure 7a), while on 5 December, only one day after, the plume directions changed becoming almost perpendicular to the coastline for the Sinni, Agri and Bradano.The Basento plume direction, on the other hand, did not change significantly, while the Cavone plume was hardly distinguishable on 5 December (Figure 7b).The high spatial variability of plumes directions over BRCW suggest high superficial current variations in the study area.
Figure 8 shows the results obtained when applying ⊗ SPM (x, y, t) (Equation ( 5)) to the SPM maps (Figure 7) of 4, 5 and 6 December 2013.In Figure 8, ⊗ SPM (x, y, t) values higher than 2 have been depicted in different colors on the basis of their magnitude, while grey pixels correspond to no-data value, as in Figure 7.The maps show a significant amount of anomalous pixels with very high values all along the coast in particular for the firsts two days analyzed.On 4 December (Figure 8a), ~21% of the total area (712 km 2 for a total area of 3390 km 2 ) shows an ⊗ SPM (x, y, t) higher than 2 and ~18.6% of the total area has an ⊗ SPM (x, y, t) higher than 3. Bearing in mind that ⊗ SPM (x, y, t) values higher than 3 have an occurrence probability of 0.15%, it is worth mentioning that pixels showing such a value undergo significant anomalies compared to the normal conditions.Similar results were found for the image of 5 December, where about 19.6% of the BRCW pixels show ⊗ SPM (x, y, t) values higher than 3 (Figure 8b), a spatial extension greater than that found on 4 December.The situation, in terms of relative SPMc variation, seems to progressively return to almost normal conditions on 6 December, when areas characterized by ⊗ SPM (x, y, t) values higher than 3 are mainly located at the Basento and Bradano mouths and in the BRCW South-West part (Figure 8c), affecting less than 1% of the entire study area.Values of ⊗ SPM (x, y, t) higher than 2 on 6 December, instead, are mainly located in the South-West portion of the study area (Figure 8c) and seem to be related to the Sinni plume, extending over a distance of 20 km from the coastline.It is worth noting that the BRCW portion indicated by the green arrow in Figure 8a probably indicates both an SPM accumulation area caused by anticlockwise currents [32] on 4 and 5 December as well as a local resuspension effect related to the presence of the small seamount shown in Figure 1b (in the South-West sector of the image).This area is linked with the main anomalous areas (in correspondence of coastal waters in front of the main rivers) through a thin portion of anomalous pixels very close to the coast (blue arrow in Figure 8a).
Remote 2016, 8, 922 13 of 20 (up to 6.5 km, highlighted by the black dashed vertical line in Figure 9), the SPMc (blue line) is quite high but never higher than the (μSPM(x, y) + 2 × σSPM(x, y)) value (purple line) so that no anomaly is detected.In this first portion, SPMc and its variability (and therefore μSPM(x, y) and σSPM(x, y)) are normally higher than in the remaining part because of the proximity to the coastline and river mouths, which makes those areas usually affected by high mean and variability values.Moving on the right side of the black dashed line, even if SPMc is lower than in the first portion of the transect, it is high enough to be greater than (μSPM(x, y) + 2 × σSPM(x, y)) (sometimes also higher than 3 × σSPM(x, y)).Hence, statistically significant ⊗SPM(x, y, t) values are obtained, indicating that for the period and location analyzed anomalous SPMc is detected.7 and will be used in Figure 9.
On 5 December, when the general currents have apparently changed, the area indicated by the green arrow in Figure 8a is disconnected from the other anomalous pixels, which move, in the case of the Sinni and Agri river plumes, to the North-East direction.On 6 December only very few anomalous pixels are still visible, extending in a rectilinear shape for more than 10 km from the coastline.
Although Figures 7 and 8 highlight similar SPMc features, the results reported in these two figures are supplementary.Figure 7a shows where high absolute SPMc values are located, but without any knowledge about the history of the area at pixel level, it is difficult to attribute a significance to these SPMc values.On the other hand, RST allows us to identify areas with SPMc values significantly higher than the value observed in normal conditions, taking also into account the signal variability (standard deviation).Therefore, the same SPMc absolute value can generate very different ⊗ SPM (x, y, t) magnitudes.To better understand the differences between the anomalies obtained computing ⊗ SPM (x, y, t) and the SPMc values, the signal profiles along the A-B transect in Figures 7c and 8c are reported in Figure 9.In detail, in this figure, the ⊗ SPM (x, y, t) variation is depicted in red and the SPMc one in blue.The purple line provides the value of SPM temporal mean plus two times the standard deviation (i.e., µ SPM (x, y) + 2 × σ SPM (x, y)) for each pixel of the profile investigated.Moreover, the red dashed line indicates the ⊗ SPM (x, y, t) threshold of 2. In the first part of the transect (up to 6.5 km, highlighted by the black dashed vertical line in Figure 9), the SPMc (blue line) is quite high but never higher than the (µ SPM (x, y) + 2 × σ SPM (x, y)) value (purple line) so that no anomaly is detected.In this first portion, SPMc and its variability (and therefore µ SPM (x, y) and σ SPM (x, y)) are normally higher than in the remaining part because of the proximity to the coastline and river mouths, which makes those areas usually affected by high mean and variability values.Moving on the right side of the black dashed line, even if SPMc is lower than in the first portion of the transect, it is high enough to be greater than (µ SPM (x, y) + 2 × σ SPM (x, y)) (sometimes also higher than 3 × σ SPM (x, y)).Hence, statistically significant ⊗ SPM (x, y, t) values are obtained, indicating that for the period and location analyzed anomalous SPMc is detected.

Evaluation of Wind Effect on SPMc Variation
To better investigate the wind-wave resuspension effect in the SPMc variation over BRCW a long-term study was performed for the month of December in the 2011-2015 period.The monthly averaged wind data recorded at the Taranto station were compared with SPMc and river levels (derived by the data illustrated in Figure 2a) aggregated at the same temporal scale (Table 1).7c and 8c).The black vertical dashed line highlights the first occurrence of anomalous signal.

Evaluation of Wind Effect on SPMc Variation
To better investigate the wind-wave resuspension effect in the SPMc variation over BRCW a long-term study was performed for the month of December in the 2011-2015 period.The monthly averaged wind data recorded at the Taranto station were compared with SPMc and river levels (derived by the data illustrated in Figure 2a) aggregated at the same temporal scale (Table 1).A clear agreement is observed in Table 1 between the variation of SPM concentrations and river levels, while no correlation seems to be present between SPM and wind.In particular, for the 5-year period investigated, the minimum wind speed value was observed in December 2013 (2.58 m/s), in correspondence with the maximum SPMc (2.75 g/m 3 ) and river level (2.13 m) values.Besides, the maximum wind speed value (4.61 m/s) was observed in December 2015, when a minimum in SPMc (1.20 g/m 3 ) and river level (1.03 m) were registered.Based on these data it can be speculated that, if present, the wind induced resuspension phenomena have only a second-order effect on the SPMc variations for CWBR.These SPMc variations, once again, are much more related to river level increasing caused by rainfall events such as the "Ciclone Nettuno" storm, confirming the previous assumption on the wind speed effect provided in Section 2.

Confutation Analysis
The analysis for a day of low SPMc was also conducted to assess RST performances in case of normal river discharges.The scope of this analysis is to show that no "false positive" value is detected in case of low/normal SPMc.This assessment was carried out on 13 December 2013, when the rivers outflowing in BRCW were no longer affected by the flood episode of the early December 2013 (Figure 2b).In detail, the mean rivers level was 2.22 m, well below those recorded on the three days already analyzed in Section 4.2 (5.21 m for 4 December 2013, 4.43 m for 5 December 2013 and 3.73 m for 6 December 2013).SPM concentrations for this day are lower than 10 g/m 3 in almost all the BRCW, with the highest values (the SPMc values do not exceed 50 g/m 3 ) obtained at the Bradano mouth (Figure 10a).Results obtained applying RST for the same day show no significant anomaly-⊗ SPM (x, y, t) < 2 -(Figure 10b).Neither are anomalies detected at the Bradano mouth where the highest SPMc values are registered.This is because the computed SPMc is in the normal range of values and variability for that specific area and period of the year.As the river levels recorded on 13 December are not affected by high river water level (Figure 2b), this result confirms once again the ability and the robustness of the RST approach to detect only the pixels affected by anomalous SPMc values.
To assess if the result just shown can be related to a specific contingency, we analyzed other two MODIS images acquired on 30 December 2012 and 13 December 2014, hence in different years from 2013 (Figure 11).The mean river levels were 1.20 m for the first image, and 1.02 m for the latter, indicating again the presence of normal river discharge over the BRCW.For both the acquisitions, SPMc was everywhere lower than 10 g/m 3 (Figure 11a,c), and no significant anomaly has been detected by RST (Figure 11b,d), confirming the previous results.10a).Results obtained applying RST for the same day show no significant anomaly-⊗SPM(x, y, t) < 2 -(Figure 10b).Neither are anomalies detected at the Bradano mouth where the highest SPMc values are registered.This is because the computed SPMc is in the normal range of values and variability for that specific area and period of the year.As the river levels recorded on 13 December are not affected by high river water level (Figure 2b), this result confirms once again the ability and the robustness of the RST approach to detect only the pixels affected by anomalous SPMc values.To assess if the result just shown can be related to a specific contingency, we analyzed other two MODIS images acquired on 30 December 2012 and 13 December 2014, hence in different years from 2013 (Figure 11).The mean river levels were 1.20 m for the first image, and 1.02 m for the latter, indicating again the presence of normal river discharge over the BRCW.For both the acquisitions, SPMc was everywhere lower than 10 g/m³ (Figure 11a,c), and no significant anomaly has been detected by RST (Figure 11b,d), confirming the previous results.

Identification of the Most Critical Areas in Terms of SPMc Values
The RST-based analysis performed at daily scale to study the December 2013 event was extended to the whole December 2003-2015 dataset to identify the BRCW areas most seriously affected by anomalous SPMc values.To this aim ⊗SPM(x, y, t) was computed all over the time period considered, and the ⊗SPM(x, y, t) > 3 occurrence frequency was computed for each pixel of the study area.Results are reported in Figure 12 in terms of percentage of occurrence with respect to the total images considered for the reference field computation at pixel level (Figure 5c).Besides some spurious and isolated pixels far from the coast, the higher frequency of occurrence was found near the coast, with a pattern very similar to that obtained studying the December 2013 event.Orange-red pixels in Figure 12 highlight those pixels where [frequency⊗SPM(x, y, t) < 3] is higher than 20%, which are mainly located in front of the Bradano, Basento and Cavone river mouths.Another area of persistence, located South of the Sinni mouth (black arrow in Figure 12), does not seem to be directly connected to any river mouth and is likely to correspond to an SPM accumulation due to the currents or resuspension phenomena.The dispersion direction of each river plume is also observable-the Bradano, Basento and Cavone plumes follow the North-East/South-West direction, while the Agri and Sinni plumes are hardly distinguishable.The highest ⊗SPM(x, y, t) > 3 occurrence frequency was  The map reported in Figure 12, referring to the whole 2003-2015 period, can be very helpful to the regional agencies for environmental monitoring because it indicates the most critical areas in terms of anomalous SPMc values.Those areas have been subject to a water transparency reduction as well as other phenomena, such as pollutant transport that could have affected the marine environment status.In the light of this, among the possible users, aquaculture companies can find this map really useful in order to exclude the most vulnerable areas from their potential production facilities.Such information, obtained by applying RST approach to SPMc maps, would have been

Identification of the Most Critical Areas in Terms of SPMc Values
The RST-based analysis performed at daily scale to study the December 2013 event was extended to the whole December 2003-2015 dataset to identify the BRCW areas most seriously affected by anomalous SPMc values.To this aim ⊗ SPM (x, y, t) was computed all over the time period considered, and the ⊗ SPM (x, y, t) > 3 occurrence frequency was computed for each pixel of the study area.Results are reported in Figure 12 in terms of percentage of occurrence with respect to the total images considered for the reference field computation at pixel level (Figure 5c).Besides some spurious and isolated pixels far from the coast, the higher frequency of occurrence was found near the coast, with a pattern very similar to that obtained studying the December 2013 event.Orange-red pixels in Figure 12 highlight those pixels where [frequency⊗ SPM (x, y, t) < 3] is higher than 20%, which are mainly located in front of the Bradano, Basento and Cavone river mouths.Another area of persistence, located South of the Sinni mouth (black arrow in Figure 12), does not seem to be directly connected to any river mouth and is likely to correspond to an SPM accumulation due to the currents or resuspension phenomena.The dispersion direction of each river plume is also observable-the Bradano, Basento and Cavone plumes follow the North-East/South-West direction, while the Agri and Sinni plumes are hardly distinguishable.The highest ⊗ SPM (x, y, t) > 3 occurrence frequency was found in front of the Bradano river mouth (47%).On the other hand, an area with a low frequency of SPMc anomaly occurrence is located in front of the Agri mouth, in agreement with the river levels measured.
The map reported in Figure 12, referring to the whole 2003-2015 period, can be very helpful to the regional agencies for environmental monitoring because it indicates the most critical areas in terms of anomalous SPMc values.Those areas have been subject to a water transparency reduction as well as other phenomena, such as pollutant transport that could have affected the marine environment status.In the light of this, among the possible users, aquaculture companies can find this map really useful in order to exclude the most vulnerable areas from their potential production facilities.Such information, obtained by applying RST approach to SPMc maps, would have been difficult to highlight by analyzing the absolute SPMc at daily or monthly scales.It is worth mentioning that such an analysis should be extended to all the other calendar months to have a clearer view of the whole area.

Conclusions and Prospective
Complex coastal areas characterized by the co-existence of several river mouths can be affected by a large spatiotemporal SPMc variability.Inferring the role played by each river in these

Conclusions and Prospective
Complex coastal areas characterized by the co-existence of several river mouths can be affected by a large spatiotemporal SPMc variability.Inferring the role played by each river in these ecosystems, in terms of SPMc relative variability, can furnish a better understating of the water quality status.In this work, ocean color satellite data were used for monitoring SPM concentration variation in the Basilicata region coastal waters (BRCW), located in the Gulf of Taranto (Ionian Sea).It is the first time that this area, characterized by the presence of 5 river mouths on about 35 km of coastline, has been object of this kind of analysis.To better characterize the BRCW behavior, in terms of SPMc variability, long-and short-time analyses have been carried out.First, data acquired by MODIS-Aqua at 250 m of spatial resolution in the red and NIR spectral bands have been used to generate SPMc daily maps by adapting an algorithm from literature [10].After that, 13 years of SPMc maps acquired in the month of December from 2003 to 2015 have been analyzed following the Robust Satellite Techniques methodology prescriptions.Such a long-term study showed on one hand that the SPMc has got a weak decreasing trend, on the other that December 2013 shows the highest value in correspondence of the highest river levels recorded in the area.In fact, the "Ciclone Nettuno" storm affected a large portion of the Basilicata region from the end of November to the first days of December 2013, causing a large flood event and consequently an increase in SPMc.A negligible SPMc resuspension contribution was confirmed during the investigated period by a specific analysis on wind data affecting the area.The RST short-term temporal analysis focused on this event has shown the potential of this methodology in integrating information contained in the daily SPMc maps.Evident anomalies have been found in correspondence of SPMc significant statistical increase.RST helps to better understand the meaning of the recorded absolute SPMc values.The same SPMc value, indeed, can be representative of very different (and stressful) environmental conditions depending on the areas investigated.As shown in Figure 9, an SPMc value up to 10 g/m 3 found close to the coastline where higher mean SPMc values and variability are usually present, cannot generate any signal anomaly, while an SPMc of 5 g/m 3 can produce anomalies where concentrations measured under unperturbed conditions are typically lower.The reliability of the proposed approach has been further assessed by a falsification analysis carried out in a few days without significant river water levels.Results achieved confirmed the RST robustness, avoiding the detection of "false positive" values in normal conditions.Finally, RST has been used for defining the areas most prone to an SPMc anomalous increasing over BRCW for the month of December in the 2003-2015 period.The most exposed areas have been found in front of the river mouths, in particular for the Bradano, Basento and Cavone rivers.Very few persistent anomalous pixels have been found associated to the Sinni and Agri mouths.In conclusion, the RST implemented on SPM data has proved to be a useful tool to better understand the behavior of coastal waters and improve monitoring activities.A major limitation of such an approach applied to ocean color satellite data is related to the presence of clouds, which obviously can mask the coastal area preventing the method application.Clearly, the spatial resolution of the satellite data used also affects the capability to provide information for the areas closest to the coastline.The performed analysis allow us to assess, as stated above, the presence of a weak trend in the investigated time series that has not significantly influenced RST sensitivity in this specific case.It is clear that, in presence of more significant trends, their effect should be removed before implementing RST, otherwise they will reduce the reliability of the achieved results.
These first results suggest several future activities, such as the extension of the analysis over the whole time period of available MODIS-Aqua data, on monthly or seasonal scales.This will help the understanding of annual and interannual SPMc variations over the BRCW as well as a clearer view of the most critical areas.The RST implementation on MODIS-Terra data will guarantee two acquisitions in a very short time period (i.e., within about 3 h) to timely detect quick SPMc variations.Furthermore, the same approach could be easily applied to other satellite data such as VIIRS-SNPP and in the future also to OLCI-Sentinel 3. Finally, owing to the intrinsic RST properties, it could be easily exported and implemented in other geographical areas, also analyzing other marine environments such as estuaries.

Figure 1 .
Figure 1.(a) Study area localization (red box) in the Mediterranean Sea; (b) magnification of the Basilicata Region Coastal Waters (BRCW) within the red box of (a), with the main rivers of the region and the bathymetry from 10 to 100 m.The black dots on the rivers indicate the position of the gauging stations.

Figure 1 .
Figure 1.(a) Study area localization (red box) in the Mediterranean Sea; (b) magnification of the Basilicata Region Coastal Waters (BRCW) within the red box of (a), with the main rivers of the region and the bathymetry from 10 to 100 m.The black dots on the rivers indicate the position of the gauging stations.

Figure 2 .
Figure 2. (a) Water levels measured at the five gauging stations closest to the Ionian Sea for the 01 January 2011-30 June 2015 period.For the Cavone river data are missing from 1 March 2011 to 1 January 2012, while for the Agri river the measurements started from 1 January 2013.(b) Magnification on the 01 November 2013-31 December 2013 period, within the black dashed rectangle in (a).

Figure 2 .
Figure 2. (a) Water levels measured at the five gauging stations closest to the Ionian Sea for the 1 January 2011-30 June 2015 period.For the Cavone river data are missing from 1 March 2011 to 1 January 2012, while for the Agri river the measurements started from 1 January 2013.(b) Magnification on the 1 November 2013-31 December 2013 period, within the black dashed rectangle in (a).

Figure 3 .
Figure 3.The registered wind intensity (m/s) and direction measured every 10 min in the 1-13 December 2013 period at the Taranto station (data from [35]).

Figure 3 .
Figure 3.The registered wind intensity (m/s) and direction measured every 10 min in the 1-13 December 2013 period at the Taranto station (data from [35]).

Figure 4 .
Figure 4. Flowchart of the various steps of elaboration.(a) Pre-processing of the MODIS L1A data to obtain SPM daily maps; (b) main steps of RST implementation.

Figure 4 .
Figure 4. Flowchart of the various steps of elaboration.(a) Pre-processing of the MODIS L1A data to obtain SPM daily maps; (b) main steps of RST implementation.

Figure 5 .
Figure 5. December 2003-2015 RST-based interannual reference fields for BRCW: (a) μSPM(x, y); (b) σSPM(x, y); (c) number of images used for each pixel to compute μSPM(x, y) and σSPM(x, y); (d) histogram of ALICE (Equation (5)) values computed for the whole time series in the 10 × 10 pixels white box reported in (c) with the theoretical Gaussian curve superimposed (blue line in (d)).

Figure 5 .
Figure 5. December 2003-2015 RST-based interannual reference fields for BRCW: (a) µ SPM (x, y); (b) σ SPM (x, y); (c) number of images used for each pixel to compute µ SPM (x, y) and σ SPM (x, y); (d) histogram of ALICE (Equation (5)) values computed for the whole time series in the 10 × 10 pixels white box reported in (c) with the theoretical Gaussian curve superimposed (blue line in (d)).

Figure 6 .
Figure 6.Trend of SPMc monthly mean spatially averaged over BRCW for the months of December from 2003 to 2015 ( ̅ , blue rhombus) compared with the interannual mean ( ̅ , dashed red line) and interannual mean ± standard deviation ( ̅ ± , dashed green line).The vertical error bar associated at each year represents the annual standard deviation spatially averaged over BRCW.The dashed black line indicates the trend line of the monthly SPMc mean for December.

Figure 6 .
Figure 6.Trend of SPMc monthly mean spatially averaged over BRCW for the months of December from 2003 to 2015 (µ SPM−M , blue rhombus) compared with the interannual mean (µ SPM , dashed red line) and interannual mean ± standard deviation (µ SPM ± σ SPM , dashed green line).The vertical error bar associated at each year represents the annual standard deviation spatially averaged over BRCW.The dashed black line indicates the trend line of the monthly SPMc mean for December.

Figure 7 .
Figure 7. SPMc maps over the BRCW obtained with the approach described in Section 3.1 on (a) 4; (b) 5 and (c) 6 December 2013.The black A-B segment in (c) will be used in Figure 9.

Figure 7 .
Figure 7. SPMc maps over the BRCW obtained with the approach described in Section 3.1 on (a) 4; (b) 5 and (c) 6 December 2013.The black A-B segment in (c) will be used in Figure 9.

Figure 8 .
Figure 8. ⊗SPM(x, y, t) maps obtained applying the RST approach for the MODIS-Aqua images acquired on: (a) 4; (b) 5 and (c) 6 December 2013.The black A-B segment in (c) is the same reported in Figure7and will be used in Figure9.

Figure 8 .
Figure 8. ⊗ SPM (x, y, t) maps obtained applying the RST approach for the MODIS-Aqua images acquired on: (a) 4; (b) 5 and (c) 6 December 2013.The black A-B segment in (c) is the same reported in Figure7and will be used in Figure9.

Figure 12 .
Figure 12.Frequency of occurrence of ⊗ SPM (x, y, t) > 3 for BRCW in the analyzed period.

Table 1 .
Monthly average SPM, river level and wind values for the month of December in the 2011-2015 period.

Table 1 .
Monthly average SPM, river level and wind values for the month of December in the 2011-2015 period.