Detecting the Effects of Extreme Events on Estuarine Suspended Particulate Matter Using Satellite Remote Sensing (Scheldt Estuary): Challenges and Opportunities

: Suspended Particulate Matter (SPM) plays an important role in controlling water quality, particularly in turbid estuaries. SPM may be impacted by changes in weather and climate, including potential changes in the frequency or intensity of extremes. Yet, the linkages between extreme events in wind and river discharge, and particularly the role these have on local dynamics and spatial patterns of estuarine SPM are, to date, largely unknown. This study investigates the effects that wind and river discharge have on SPM in a turbid estuary. It uses atmospherically corrected remotely sensed reﬂectances from Landsat-5, -7, and -8. From these data, we establish a 37-year-long time series, from 1984 to 2020, of satellite-derived SPM for the Scheldt Estuary. SPM was estimated using an algorithm applied to the near-infrared band and subsequently compared to in situ SPM data. Additionally, a time series of river discharge and wind speed were used to assess the frequency and severity of extreme events of wind speed and river discharge. In general, statistically signiﬁcant but weak relationships between SPM and river discharge and between SPM and wind speed were observed. SPM correlated with river discharge and wind in large parts of the estuary, indicating the important role of these drivers in the entire estuarine system. Our study demonstrates how synoptic satellite snapshots can be combined with in situ time series of drivers, such as river discharge and wind, to capture where these drivers relate to (and likely affect) SPM within an estuary. However, our study demonstrates an inability to capture SPM during windstorms both from in situ and satellite data. We discuss the challenges and limitations of assessing the effects of extreme events from satellite and in situ SPM. We recommend the deployment of complementary moored turbidimeters for continuous observations at two strategic locations, as indicated by our spatial study, to improve the ability to capture the effects of extreme events in both space and time.


Introduction
Suspended particulate matter (SPM) is a critical variable in determining water quality [1] and the functioning of aquatic systems. It plays an important role in the fate of pollutants and pathogens in estuaries [2,3], modulates light availability for primary productivity [4,5], benthic biota [6], and contributes to the obstruction of waterways [7].
Within estuaries, SPM is expected to be a response to complex interactions among tides, river discharge, and winds. While in microtidal estuarine environments, river discharge and winds are known as the main forcing mechanisms of SPM dynamics; in meso-to macrotidal environments, the tidal regime and the asymmetry between ebb and flood currents are commonly the main driving mechanisms of changes in SPM [8]. For example, the association between ebbing tidal currents and river discharge increases the flow of SPM to the coastal zone, while flooding tidal currents lead to the increased residence time of freshwater and modulate SPM entrainment in turbid maximum zones, or TMZ [9]. The difference between these processes in a tidal cycle indicates residual sediment transport. Although tidal currents and river discharge typically control the distribution and fate of SPM, wind-driven waves can induce erosion and resuspension of sediments [10,11]. Adding to the complexity of the hydrodynamical relationship with SPM in tidal estuaries, the combined effect of tides, rivers, and wind-driven circulation typically creates a constantly changing bathymetry in soft-bottom estuaries [12].
Estuaries are among the areas most susceptible to climate change and extreme meteorological events. Examples include cyclones, storm surges, and extreme rainfall, which are likely to increase in the decades to come [13]. These extreme events impact several national security goods and services due to higher waves, floods, and severe wind conditions affecting navigability and enhancing changes in SPM [14,15]. In this scenario, a long-lasting effect of remobilized SPM may, for example, lead to a change in the trophic state [16].
Therefore, quantifying and understanding the role of SPM and extricating what affects SPM is particularly pressing for estuaries. Most of the studies addressing these questions rely on hydrodynamical modeling (e.g., references [17][18][19]) and water sampling (e.g., [20,21]); however, a growing number of studies have developed a foundation for the use and application of satellite remote sensing to understand SPM dynamics in turbid and tidal-dominated estuarine environments (e.g., references [10,[22][23][24]).
The main objective of this study is to explore the potential of long-term monitoring by high spatial resolution satellite-borne sensors to quantify the effects of extreme events on the spatial variability in surface SPM in a complex estuarine environment. We assess the relationships of SPM with wind and river discharge to identify the main spatial patterns and to attribute the contribution of each environmental driver in a spatial context. Our goal will be achieved by performing a long-term spatiotemporal analysis (from 1984 to 2020) of remote sensing data (Landsat-5, -7, and -8).

Study Site
The Western Scheldt Estuary is an open estuary on the southwestern part of the Dutch coast. The estuary is characterized as a funnel-shaped multiple-channel system, separated by elongated shallow areas and tidal flats. The Western Scheldt, with its multiple-channel system, forms one estuarine system with the adjacent one-channel Sea Scheldt in Belgium, together referred to as Scheldt Estuary ( Figure 1). The Scheldt Estuary functions as an essential area for navigation, recreation, and natural habitat for, among others, birds, fish, and saltmarsh vegetation [25].
Occasional capital dredging and regular maintenance dredging are carried out to secure safe navigation to the port of Antwerp for increasingly larger ship vessels [21]. Dredged material is redistributed either in the deeper areas, in secondary channels, or used to counteract the erosion of shallow subtidal zones near tidal flats [25,26]. Sand mining has also played a role in sediment transport in the Scheldt [7]. The sediment consists of mainly fine inorganic sediments [27,28] ranging from sandy-clay mud in shallower areas, to sand in the deepest channel [29]. Cohesive organo-mineral particles are observed in the TMZ (Turbidity Maximum Zone) [30]. At least four TMZs were reported in the Scheldt (see references [8,21,28,30,31]), with the most prominent TMZ located in the mid estuary (between about 58 and 110 km from the estuary mouth [8]).
The sediment dynamics in the Scheldt are influenced by tides, river discharge, and winds [20,23]. These forcings divide the estuary into three parts governing SPM transport and distribution patterns [20]: (i) the lower estuary (0-45 km), where tides and wind-driven waves are dominant forcing mechanisms, and a net landward SPM transport occurs; (ii) mid estuary (45-75 km), where tidal energy dominates and SPM from the lower and upper estuary converges, and (iii) the upper estuary (75-95 km), where tides and river discharge are governing the dynamics and river discharge may be the prevalent forcing mechanism Remote Sens. 2023, 15, 670 3 of 17 in winter months. In the upper estuary, a net seaward SPM transport takes place, and the main sources of sediment to the estuary are the eroded tidal flats and shoals and SPM carried by river discharge. and distribution patterns [20]: (i) the lower estuary (0-45 km), where tides and winddriven waves are dominant forcing mechanisms, and a net landward SPM transport occurs; (ii) mid estuary (45-75 km), where tidal energy dominates and SPM from the lower and upper estuary converges, and (iii) the upper estuary (75-95 km), where tides and river discharge are governing the dynamics and river discharge may be the prevalent forcing mechanism in winter months. In the upper estuary, a net seaward SPM transport takes place, and the main sources of sediment to the estuary are the eroded tidal flats and shoals and SPM carried by river discharge. ). The yellow hatched area depicts the region of TMZ occurrence (within this study's boundaries). Regions in light yellow represent no bathymetric data available for the Sea Scheldt (BEL). Sites I and II approximately indicate the sites of suggested installation/deployment of high temporal resolution turbidity sensors discussed in Section 4.4.
The Scheldt Estuary is under the influence of a semi-diurnal meso to macrotidal regime, with tides reaching up to 4.85 m on spring tides and 2.81 m during neap tides [20,27]. The tidal range increases regularly in the upstream direction, with a maximum range registered as 5.25 m. Another important aspect of the tidal dynamics in the Scheldt is its asymmetry. This is possibly the main factor influencing sediment exchange between the estuary and the adjacent coast [32] during low river discharge [33].
River discharge varies from 2 m 3 s −1 to 600 m 3 s −1 (respectively in summer and winter months [34]) with an annual average between 100 and 200 m 3 s −1 [7]. During peaks in river discharge, the Scheldt changes from mostly a well-mixed to a partially stratified estuary, which leads to increased retention of river sediments from the increased flow of sediment- The Scheldt Estuary is under the influence of a semi-diurnal meso to macrotidal regime, with tides reaching up to 4.85 m on spring tides and 2.81 m during neap tides [20,27]. The tidal range increases regularly in the upstream direction, with a maximum range registered as 5.25 m. Another important aspect of the tidal dynamics in the Scheldt is its asymmetry. This is possibly the main factor influencing sediment exchange between the estuary and the adjacent coast [32] during low river discharge [33].
River discharge varies from 2 m 3 s −1 to 600 m 3 s −1 (respectively in summer and winter months [34]) with an annual average between 100 and 200 m 3 s −1 [7]. During peaks in river discharge, the Scheldt changes from mostly a well-mixed to a partially stratified estuary, which leads to increased retention of river sediments from the increased flow of sediment-rich river discharge [28]. Meanwhile, because the flow of marine water upstream is lessened, so is the supply of marine sediment upstream and the intensity of the TMZ.
Wind climate in the Scheldt is characterized by stronger and more frequent south-west (SW) winds and by north-east (NE) winds [35]. The results in reference [35] show that during persistent SW winds, it is likely that the Scheldt Estuary imports marine mud from the Belgian-Dutch inner shelf, which can lead to increased SPM in the lower estuary. Wind-driven waves can also play a significant role in sediment transport in the estuary [20]. This mechanism is likely to be effective in the estuarine mouth (i.e., lower estuary), whereas locally generated wind-driven waves are expected in the middle estuary [36]. The effect of such waves is assumed to be local and dependent on the wind direction. During severe windstorms, wave height can reach peaks up to 3 m [36].

In Situ Data
In situ data comprise long time series of SPM data and the drivers of estuarine dynamics: water level as a proxy for tidal stage, river discharge, and wind speed and direction (see overview in Table S1).
SPM data provided by Rijkswaterstaat (hereafter RWS, https://waterinfo.rws.nl/#! /nav/index, accessed on 7 April 2021) consists of stations distributed along the Scheldt Estuary domain and covering different temporal ranges (refer to Table S1). The four longterm monitoring stations Vlissingen Boei SSVH, Hansweert Geul, Terneuzen boei 20, and Schaar van Ouden Doel are further in this document referred to as RWS Vlissingen, RWS Hansweert Geul, RWS Terneuzen, and RWS Schaar van Ouden Doel, respectively. The second set of SPM data is provided by [37] and distributed by the R for Science project (https://www.rforscience.com/data/westerschelde/, accessed on 16 March 2021) of NIOZ Royal Netherlands Institute for Sea Research. This dataset consists of seventeen stations sampled along the central axis of the Scheldt Estuary from 1995 to 2004.
Among the drivers of estuarine dynamics, information on water level and river discharge is also available from the RWS database. Water level data (https://waterinfo. rws.nl/#!/nav/index/, accessed on 23 November 2021) were continuously recorded every 10 min from 1984 to 2020 at Hansweert Geul station. We used the estimated time difference from these water level records in [23] to determine the tidal stage (i.e., low tide, high tide, ebb, and flood) associated with the SPM data. For example, high water level at RWS Vlissingen occurs 56 min prior to higher water level at RWS Hansweert Geul (reference station), 37 min prior at RWS Terneuzen, and 34 min after at RWS Schaar van Ouden Doel. The same time differences were applied to the four different regions, centered at each RWS in situ station. Daily river discharge was recorded at station Schaar van Ouden Doel from 1984 to 2020. Extreme events in river discharge (Q) were classified as extremely high when discharge was equal to or higher than the 95th percentile of the entire time series. Wind data were obtained from Royal Netherlands Meteorological Institute (KNMI, https://www-p.knmi.nl/nederland-nu/klimatologie/uurgegevens, accessed on 16 October 2020) from Vlissingen (station 310) with a 1 h temporal resolution and covering the period between 1984 and 2020. To obtain daily, monthly, and yearly wind data, wind speed (U) and direction cannot be directly averaged. First, wind components (u and v) must be calculated; then, each component can be averaged to the needed timestep, and subsequently converted back to wind speed and direction. The stronger the wind speeds, the higher the discrepancy between vector and scalar averages. KNMI defines extremely high wind (i.e., windstorms) when wind speed is equal to or higher than 20.8 m·s −1 (https: //www.knmi.nl/kennis-en-datacentrum/uitleg/storm, accessed on 16 October 2020).

Remote Sensing Data and Validation of SPM Products
A time series of all available, full resolution, level 1 data from Landsat-5 (1984Landsat-5 ( -2011, Landsat-7 (1999Landsat-7 ( -2003, and Landsat-8 (2013-2020) were processed using the Royal Belgium Institute of Natural Sciences (RBINS) processor ACOLITE (https://odnature.naturalsciences. be/remsem/software-and-data/acolite, accessed on 20 December 2022; version 20210114.0), using the default dark spectrum fitting scheme [38,39] to obtain atmospherically corrected remotely sensed reflectance (Rrs) for the Scheldt Estuary (swath inside the area 3.19 • -4.42 • E and 51.15 • -51.49 • N, refer to Figure 1). Every satellite scene was mapped to a cylindrical equidistant projection at a spatial resolution of 30 m. All scenes were then manually inspected with the help of true-color composites, and cloud-compromised scenes were removed, leaving 178 good-quality scenes (114, 22, 42 scenes for Landsat-5, -7, and -8, respectively). SPM was then estimated using the semi-analytical algorithm described in reference [40], Equation (1): where ρ w (ρ w = Rrs · π) is the water-leaving reflectance, and coefficients A and C are hyperspectral narrow-band coefficients originally provided by reference [40]. Here, we test and compare the application of A and C coefficients as originally provided and bandweighted coefficients from the original narrow-band coefficients (provided in Table S2). The choice of algorithm is motivated by the fact that it was developed from an in situ dataset collected in the southern North Sea. Satellite-derived SPM was calculated from both the red and the near-infrared (NIR) bands to check for the best suitable satellite band and prevent underestimation of SPM due to possible saturation of Rrs signal at high SPM [41]. Following the good-practices criteria proposed in [3], the median and standard deviation of satellite-derived SPM were extracted by applying 3 × 3-pixel boxes centered at the location of each in situ sampling site. Additionally, two performance assessments were carried out to analyze the algorithm retrievals at both red and NIR bands (see Section 3.1). Firstly, median and standard deviation parameter from SPM data retrieved from the satellite was compared with in situ SPM data collected within 30 min and within 60 min of each satellite scene (hereafter called match-up analysis). The shorter time window may provide more accurate estimates, which can be significant given the strong hydrodynamical forcings the Scheldt Estuary is under, while a longer time window increases the number of data match ups [3,42]. For validation metrics, a type II linear regression [43] least squares fit (MATLAB script by [44]) was applied [3]. Additional statistical descriptors such as the mean absolute error, MAE (g·m −3 ), root mean square error, RMSE (g·m −3 ), and mean absolute percentage error, MAPE (%), were estimated.
Secondly, the frequency distributions of full time series of both satellite-derived and in situ SPM were compared at four long-term established stations in the Scheldt Estuary: Vlissingen SSVH, Terneuzen Boei 20, Hansweert Geul, and Schaar van Ouden Doel (hereon RWS Vlissingen, RWS Terneuzen, RWS Hansweert Geul, and RWS Schaar van Ouden Doel, respectively). Location of stations is shown in panels b, g, k, and l in Figure 1. Such a comparison tests if the independent sets of SPM data yield similar frequency distribution (i.e., shape, center, spread, skewness), as this indicates that the satellite and in situ instruments observe similar patterns [3,45]. Differences may be due to calibration inaccuracies (possibly causing a bias), ineffective atmospheric correction, or differences in measurement statistics due to, e.g., cloud cover, cloud shade, or time of observation.

Data Analysis
Climatological yearly cycles were extracted from the monthly river discharge and wind speed time series, and a least square curve was fitted to each parameter. Extreme events in daily river discharge and hourly winds were analyzed for their frequency, duration, and magnitude from January to December of each year. Theil-Sen regression (MATLAB script by [46]) was carried out to establish the significance of temporal trends, applying a significance level of 5%.
Maps were then generated with SPM derived from satellite data using the NIR-based algorithm. From these SPM maps, we extracted median and standard deviation from the region surrounding each of the four long-term monitoring in situ stations (see polygons in Figure 1). The resulting satellite-derived time series of SPM were compared to time series of river discharge and wind speed using Theil-Sen regressions. Spatial correlations of SPM with wind speed and river discharge were examined, applying Kendall-tau correlations at a 5% significance level (with statistical descriptors Kendall-tau correlation coefficient (r), probability (p), and the number of observations (n). Estimated significant correlation coefficients were classified as weak if r ≤ 0.3 and moderate if 0.3 < r ≤ 0.6 [47].

Validation of SPM Products
SPM match ups and frequency distributions convey qualitative and quantitative information in different ways, both visually and statistically. The match-up analysis covered a large extent of the Scheldt Estuary (black and dark-yellow dots in Figure 1) and a considerable range of in situ SPM (4 to 147 g·m −3 ) despite the limited number of cloud-free satellite scenes matching in situ measurements. Match ups totaled n = 27, of which 18 come from Landsat-5, 7 from Landsat-7, and 2 from Landsat-8 (refer to Table S3).
A first validation analysis was carried out as a comparison of in situ SPM against the match ups for SPM derived from narrow-band and from band-weight coefficients. Results show larger differences in the NIR band with lower uncertainties from the band-weighted coefficients (i.e., improved MAPE by about 14%), and technically similar uncertainties were observed in the red band (less than 1% difference). From these results, the remainder of this study is carried out by applying band-weighted coefficients.
The Comparison of frequency distributions of all in situ SPM (at four long-term monitoring stations) to the satellite-derived SPM from the red and the NIR bands (at the location represented by the four in situ stations) allowed for more samples than for the match ups. For all four long-term monitoring stations, satellite-derived SPM distributions from the red band are underestimated by about 50% (Figure 3, distributions in dark yellow). On the other hand, satellite-derived SPM from the NIR band yielded a reasonable maximum underestimation of 25%, and distributions follow a similar shape and skewness (Figure 3, distributions in white). It corroborates the assumption that the NIR band is the most useful satellite band for applications in the Scheldt Estuary under moderate to high SPM concentrations, compared to shorter wavelengths (as suggested in [48]). Thus, the results presented further are carried out using the NIR band. Additionally, the frequency distribu- Comparison of frequency distributions of all in situ SPM (at four long-term monitoring stations) to the satellite-derived SPM from the red and the NIR bands (at the location represented by the four in situ stations) allowed for more samples than for the match ups. For all four long-term monitoring stations, satellite-derived SPM distributions from the red band are underestimated by about 50% (Figure 3, distributions in dark yellow). On the other hand, satellite-derived SPM from the NIR band yielded a reasonable maximum underestimation of 25%, and distributions follow a similar shape and skewness ( Figure 3, distributions in white). It corroborates the assumption that the NIR band is the most useful satellite band for Remote Sens. 2023, 15, 670 7 of 17 applications in the Scheldt Estuary under moderate to high SPM concentrations, compared to shorter wavelengths (as suggested in [48]). Thus, the results presented further are carried out using the NIR band. Additionally, the frequency distributions suggested that SPM concentrations increase upstream by a large factor (e.g., medians show that SPM at RWS Schaar van Ouden Doel is about 100% higher than at RWS Vlissingen).

Description and Context of Estuarine Drivers during SPM Data Acquisition
With strong interannual variability and regular seasonality, river discharge is generally low in the spring and end of summer and high in winter months (Figure 4b), while wind speed is low in late winter and spring with highs in late fall and early winter (Figure 4d). Both river discharge and wind speed peak in the winter (see Figure S4 for a qualitative interannual comparison with SPM).

Description and Context of Estuarine Drivers during SPM Data Acquisition
With strong interannual variability and regular seasonality, river discharge is generally low in the spring and end of summer and high in winter months (Figure 4b), while wind speed is low in late winter and spring with highs in late fall and early winter (Figure 4d). Both river discharge and wind speed peak in the winter (see Figure S1 for a qualitative interannual comparison with SPM). Analysis of the in situ data for extreme events of river discharge revealed no significant trends in the magnitude (Figure 5a) or in the frequency of extreme river discharge ( Figure 5b); there was a significant negative trend in time in the duration of extreme river discharge (Figure 3c; r = −0.59, p << 0.001). No significant long-term trends were detected in the magnitude (Figure 5d), frequency (Figure 5e), and duration ( Figure  5f) of extreme wind speeds (i.e., windstorms); however, a seemingly decadal oscillation of magnitude of windstorms, with a relatively weak period between 1998 and 2013, is observed (Figure 5d). Neither satellite scenes nor in situ SPM were uniformly acquired at tidal stages or seasons (Figure 6a,b). Interestingly, the heterogeneity in SPM sampling is larger for in situ data than for satellite observations of tidal stages (i.e., while ebbing and flooding tides are overrepresented from satellite scenes, there is no clear pattern in tidal bias observed from in situ samples and the heterogeneity increases in the upstream stations). Analysis of the in situ data for extreme events of river discharge revealed no significant trends in the magnitude (Figure 5a) or in the frequency of extreme river discharge ( Figure 5b); there was a significant negative trend in time in the duration of extreme river discharge (Figure 3c; r = −0.59, p < 0.001). No significant long-term trends were detected in the magnitude (Figure 5d), frequency (Figure 5e), and duration (Figure 5f) of extreme wind speeds (i.e., windstorms); however, a seemingly decadal oscillation of magnitude of windstorms, with a relatively weak period between 1998 and 2013, is observed (Figure 5d).  Analysis of the in situ data for extreme events of river discharge revealed no significant trends in the magnitude (Figure 5a) or in the frequency of extreme river discharge ( Figure 5b); there was a significant negative trend in time in the duration of extreme river discharge (Figure 3c; r = −0.59, p << 0.001). No significant long-term trends were detected in the magnitude (Figure 5d), frequency (Figure 5e), and duration ( Figure  5f) of extreme wind speeds (i.e., windstorms); however, a seemingly decadal oscillation of magnitude of windstorms, with a relatively weak period between 1998 and 2013, is observed (Figure 5d). Neither satellite scenes nor in situ SPM were uniformly acquired at tidal stages or seasons (Figure 6a,b). Interestingly, the heterogeneity in SPM sampling is larger for in situ data than for satellite observations of tidal stages (i.e., while ebbing and flooding tides are overrepresented from satellite scenes, there is no clear pattern in tidal bias observed from in situ samples and the heterogeneity increases in the upstream stations). Neither satellite scenes nor in situ SPM were uniformly acquired at tidal stages or seasons (Figure 6a,b). Interestingly, the heterogeneity in SPM sampling is larger for in situ data than for satellite observations of tidal stages (i.e., while ebbing and flooding tides are overrepresented from satellite scenes, there is no clear pattern in tidal bias observed from in situ samples and the heterogeneity increases in the upstream stations). Seasonally, satellite observations are underrepresented by about 50% in fall and winter compared to summer and spring ( Figures S5.A and S5.B) for both river discharge and wind ( Figure 7) due to cloudiness. Of all satellite observations, high extreme events were only recorded during the winter season for river discharge (Figure 7a), while high extreme events in wind speed (i.e., windstorms) were not observed within satellite overpass for the period (Figure 7b). In situ SPM was not sampled during windstorms. On the other hand, in situ SPM was sampled during extreme river discharge events in spring, fall, and winter ( Figures S5.A and S5.B). SPM was also investigated for peaks under predominant wind directions, for which NE and SW winds were observed to be the most frequent and stronger, in accordance with the literature (refer to Figure S6).

Spatial Relationships of River Discharge and Wind Speed with SPM
We observed a statistically significant daily correlation of SPM with river discharge (Kendall-tau correlation r » 0.1-0.3; Figure 8a) and a weaker correlation with daily wind speed (Kendall-tau correlation r » 0.1-0.2; Figure 8b). The TMZ-dominated area presented no significant relationship with either river discharge or wind speed for the period. Seasonally, satellite observations are underrepresented by about 50% in fall and winter compared to summer and spring ( Figure S2A,B) for both river discharge and wind (Figure 7) due to cloudiness. Of all satellite observations, high extreme events were only recorded during the winter season for river discharge (Figure 7a), while high extreme events in wind speed (i.e., windstorms) were not observed within satellite overpass for the period (Figure 7b). In situ SPM was not sampled during windstorms. On the other hand, in situ SPM was sampled during extreme river discharge events in spring, fall, and winter ( Figure S2A,B). SPM was also investigated for peaks under predominant wind directions, for which NE and SW winds were observed to be the most frequent and stronger, in accordance with the literature (refer to Figure S3). Seasonally, satellite observations are underrepresented by about 50% in fall and winter compared to summer and spring ( Figures S5.A and S5.B) for both river discharge and wind ( Figure 7) due to cloudiness. Of all satellite observations, high extreme events were only recorded during the winter season for river discharge (Figure 7a), while high extreme events in wind speed (i.e., windstorms) were not observed within satellite overpass for the period (Figure 7b). In situ SPM was not sampled during windstorms. On the other hand, in situ SPM was sampled during extreme river discharge events in spring, fall, and winter ( Figures S5.A and S5.B). SPM was also investigated for peaks under predominant wind directions, for which NE and SW winds were observed to be the most frequent and stronger, in accordance with the literature (refer to Figure S6).

Spatial Relationships of River Discharge and Wind Speed with SPM
We observed a statistically significant daily correlation of SPM with river discharge (Kendall-tau correlation r » 0.1-0.3; Figure 8a) and a weaker correlation with daily wind speed (Kendall-tau correlation r » 0.1-0.2; Figure 8b). The TMZ-dominated area presented no significant relationship with either river discharge or wind speed for the period.

Spatial Relationships of River Discharge and Wind Speed with SPM
We observed a statistically significant daily correlation of SPM with river discharge (Kendall-tau correlation r » 0.1-0.3; Figure 8a) and a weaker correlation with daily wind speed (Kendall-tau correlation r » 0.1-0.2; Figure 8b). The TMZ-dominated area presented no significant relationship with either river discharge or wind speed for the period.
Furthermore, we investigated in situ SPM in the long-term stations in comparison with satellite-derived SPM at each region surrounding the in situ stations (Figures 9 and 10). All relationships are statistically significant, and estimated intercepts are comparable (within about ±15% difference) between each satellite region (Figure 9a-d) and respective in situ stations (Figure 9e-h). Furthermore, we investigated in situ SPM in the long-term stations in comparison with satellite-derived SPM at each region surrounding the in situ stations (Figures 9 and  10). All relationships are statistically significant, and estimated intercepts are comparable (within about ±15% difference) between each satellite region (Figure 9a-d) and respective in situ stations (Figure 9e-h).  Figure 1). Panels (e)-(h) display the long-term SPM in situ stations. Shaded areas represent the region of extremely high river discharge. Tidal stages are color-coded, and the size of dots is proportional to wind speed on the day of data acquisition (see legend). Theil-Sen regressions are statistically significant (i.e., p < 0.05).

Figure 9.
Relationship between river discharge (Q) and SPM. Panels (a-d) represent the median satellite-derived SPM in the delimited region (refer to Figure 1). Panels (e-h) display the long-term SPM in situ stations. Shaded areas represent the region of extremely high river discharge. Tidal stages are color-coded, and the size of dots is proportional to wind speed on the day of data acquisition (see legend). Theil-Sen regressions are statistically significant (i.e., p < 0.05). Figure 9. Relationship between river discharge (Q) and SPM. Panels (a)-(d) represent the median satellite-derived SPM in the delimited region (refer to Figure 1). Panels (e)-(h) display the long-term SPM in situ stations. Shaded areas represent the region of extremely high river discharge. Tidal stages are color-coded, and the size of dots is proportional to wind speed on the day of data acquisition (see legend). Theil-Sen regressions are statistically significant (i.e., p < 0.05).

Figure 10.
Relationship between wind speed (U) and SPM. Panels (a)-(d) represent the median satellite-derived SPM in the delimited region (refer to Figure 1). Panels (e)-(h) display the long-term SPM in situ stations. Shaded areas represent the region of extremely high wind speed. Tidal stages are color-coded, and the size of dots is proportional to river discharge on the day of data acquisition (see legend). Theil-Sen regressions are statistically significant (i.e., p < 0.05).
While river discharge relationships with SPM are observed on the high extreme event range (grey regions in Figure 9), there was no observation from satellite nor in situ SPM carried out during extremely high wind speed settings ( Figure 10). The lack of SPM observations during windstorm conditions prevents our analysis from assessing possible responses in SPM in such circumstances. Additionally, the statistically significant relationships between SPM and wind speed are weak (and weaker than the relationships with river discharge).  Figure 1). Panels (e-h) display the long-term SPM in situ stations. Shaded areas represent the region of extremely high wind speed. Tidal stages are color-coded, and the size of dots is proportional to river discharge on the day of data acquisition (see legend). Theil-Sen regressions are statistically significant (i.e., p < 0.05).
While river discharge relationships with SPM are observed on the high extreme event range (grey regions in Figure 9), there was no observation from satellite nor in situ SPM carried out during extremely high wind speed settings ( Figure 10). The lack of SPM observations during windstorm conditions prevents our analysis from assessing possible responses in SPM in such circumstances. Additionally, the statistically significant relationships between SPM and wind speed are weak (and weaker than the relationships with river discharge).

Assessing Estuarine SPM Concentrations from Landsat Satellite Sensors
We assessed the effects of river discharge and wind on the variability in Suspended Particulate Matter (SPM) using remote sensing and in situ data. The unique 37-year data record of Landsat-5, -7, and -8 sensors, coupled with improved data processing, allows for investigating these relationships in a spatially explicit context.
The overall agreement between in situ SPM and satellite-estimated SPM concentrations confirms the accuracy of satellite observation and the used NIR-band algorithm [40] to derive SPM. In particular, the band-weighted coefficients used in the present paper give a reasonably better performance for SPM in the NIR band than the original narrow-band coefficients. Band-weighting is a more accurate technique for (broad-band) satellite-derived water quality products. This highest accuracy makes band-weighted coefficients desirable over narrow-band coefficients (e.g., reference [39]) for best algorithm performance [49]. However, band-weighting may be somewhat limited by the SRF (Spectral Response Function) spectral limits in the NIR band of Landsat-5 (760-900 nm), Landsat-7 (770-9000 nm), and Landsat-8 (829-900 nm) because their ranges are not fully covered by the spectral resolution provided in the original narrow-band coefficients (520-885 nm). This source of uncertainty was addressed in Supplementary Material S1. Our results, therefore, provide a cautious note and a recommendation on applying band-weighted coefficients for Landsat-5, -7, and -8 given the spectral data limitations.
For both NIR and red-based SPM estimates, other relevant sources of uncertainties may result from (1) space and time differences between in situ sampling and satellite overpass when sampling in patchy waters or highly dynamic environments, i.e., we observed that a time interval of 60 min might be a large interval for a system such as the Scheldt; and (2) (specific) inherent optical properties, namely absorption and backscattering of particles may vary, i.e., due to the natural or human-induced variability in particle size (e.g., formation of flocs), or particle composition (e.g., seasonal plankton bloom; refer to Supplementary Material I in reference [38]). Additionally, results showed that in the turbid waters of the Scheldt Estuary, the performance of the red band was less optimal than the NIR, which may also be attributed, to a degree, to the saturation of reflectance and the contribution of Colored Dissolved Organic Matter (CDOM) to the total SPM estimate which reference [40] observed could underestimate SPM by about 22% in the red band and by 2% in the NIR.
Finally, from the frequency histograms, the magnitudes and spatial patterns of satellitederived SPM from the NIR band were consistent with long-term in situ observations at four stations along the Scheldt Estuary. These are the reasons why we believe applying NIR bands is the most adequate choice for SPM estimates, especially under extreme event conditions.

SPM Variability from Interaction with Estuarine Forcing Mechanisms and Its Limitations
We have found that, in general, the relationship of SPM was stronger with river discharge, and weaker with wind speed (by 50%). In addition, while river discharge presents a similar spatial relationship with SPM along the estuary, winds present more zones of statistically non-significant relationships. A relatively stronger correlation was observed west of RWS Vlissingen station and in the RWS Hansweert Geul region. Our analysis also suggests large SPM variability under lower river discharge (i.e., Q < 250 m 3 ·s −1 ) and a linearly increasing trend in SPM concentrations under higher river discharge (i.e., Q > 250 m 3 ·s −1 , including extremely high river discharge episodes), while large SPM variability is observed at all wind speeds (U < 15 m·s −1 ). Wind direction may also relate to SPM variability: we observed that overall stronger wind speeds and peaks of SPM are more frequent during NE and S-SW wind conditions.
We speculate that the large SPM variability and consequently relatively low correlation coefficients observed from the interaction with river discharge and wind speeds may be explained by tidal influence and delayed effects on SPM. In the Scheldt, tides were observed to be significant across the estuary. From the studies by references [23,27], we speculate that the high variability in SPM may partially be explained by sediment transport from tidal flats (likely observed at the end of low tide and flooding); entrainment and advection processes, especially in the mid and upper estuary (i.e., presence of the TMZ); and high tidal current velocities responsible for mixing processes (in high tide and ebbing). In our study, we considered only the instantaneous effects of river discharge and wind over SPM concentrations, while a delayed effect (e.g., long settling time) of high SPM could have been captured when SPM data were sampled.
Analyses based on satellite-derived or in situ SPM are proven to have their inherent limitations. Reference [23] concluded that remote sensing by sun-synchronous sensors may overrepresent SPM at certain tidal stages, while SPM is underrepresented seasonally (i.e., winter months) due to cloud cover. From in situ SPM, we have observed a similar misrepresentation of SPM. When SPM in situ sampling is carried out on the same day (i.e., within hours among RWS Vlissingen, RWS Terneuzen, RWS Hansweert Geul, and RWS Schaar van Ouden Doel stations), changes in a tidal stage are likely observed, as an interval of just a few hours is long enough for a change in tidal stage, current direction, and SPM dynamics. Finally, in situ sampling is seasonally less biased than satellite remote sensing, yet it also underrepresents extreme wind events.

Effects of Extreme Events on SPM
Extreme events in river discharge and wind speed (windstorms) were expected to play a role in the complex SPM dynamics, especially-if not exclusively-in winter months, as river discharge is at its maximum flow, likely carrying a higher input of inland sediments in suspension, and when wind-driven waves are expected to regulate resuspension of bottom sediment and prevent SPM from settling. Our study, however, shows the limitation of data collection of SPM in extreme wind events: neither long-term high spatial resolution satellite imagery (limited overpass, likely cloud conditions during storms) nor in situ data (monitoring ships will not sail during storms) were available to capture SPM under extremely high wind conditions. Therefore, our regressions concerning wind speed and SPM do not consider extremely high wind events preventing a more detailed analysis of windstorms. In contrast, SPM observations (both in situ and satellite-derived estimates) were possible during extremely high river discharge conditions, although in a low number of samples. In that sense, extrapolation of the system's response during extreme river discharge events surpassing the range of data analyzed here is discouraged.
Our regression analysis did not reveal significant linear changes in the magnitude, frequency, and duration of windstorms revealing the stochastic nature of extreme wind events. However, we qualitatively observed decadal peaks in the magnitude of windstorms. We suggest for future studies that a longer time series may be used to investigate the occurrence of any non-linear temporal patterns of extreme wind events in the Scheldt. Extreme river discharge has experienced stepwise declines over the past decades, indicating a reduction in sediment discharge, i.e., a decline in the winter magnitude of river discharge but also a significant decline in the duration of extreme events. We speculate that this might be associated with engineering works particularly in the Sea Scheldt, such as the so-called SigmaPlan (https://www.sigmaplan.be, accessed on 21 December 2022) in Belgium, which includes the construction of natural flood-control areas in the river valleys to capture excess river water to manage floodwater levels.
Numerous studies have described the relationship and effects of estuarine drivers of SPM dynamics in multiple coastal sites, e.g., references [16,50,51], most of which, such as the Scheldt Estuary, shelter large population counts and have undergone a series of humanmade adaptations to prevent the unwanted effects of extreme events. However, due to the particularities of each system, the responses to extreme events are diverse. For example, the Mekong River Delta (Vietnam) [50] was described with positive feedback between river discharge and the amount of SPM during the monsoon season due to stronger river discharge and intensified wind-induced resuspension. In contrast, reference [51] observed a sharp loss of turbidity in the occurrence of flood events (river discharge above the 75th percentile) in the macrotidal Gironde Estuary (France), and reference [16] observed lower turbidity in the microtidal estuary of Río de la Plata (Uruguay) caused by periods of intensified river discharge due to excessive rain. However positive (i.e., increase in turbidity with extreme events) or negative (i.e., turbidity dilution) the trend with SPM, it is rational to expect that extreme events represent changes in the local SPM dynamics and that these changes likely also affect other aspects of the estuary (e.g., changes in water quality and sediment-related biota [21], or shoreline recession of tidal flats [52]). We also speculate that the effects of extreme events in the diverse biogeochemical aspects of the estuary are not uniform along the estuary, e.g., reference [53], thus affecting more the sites and regions that are not used to large variability in SPM. In such sites, the biota may be less resilient to environmental changes, or a temporary loss of species richness and species composition may be observed. For example, the dynamics may change more drastically with an above-the-records extreme event in the region surrounding Hansweert Geul station than in the region of Schaar van Ouden Doel station. The latter region is more susceptible to (frequent) large SPM variability and/or species biologically adapted to low light conditions and increased turbulent conditions in the water column, for example.

Transferable Lessons and Opportunities to Monitor SPM in Extreme Events
The effects of extreme events on SPM patterns (and patterns of other optically active parameters such as chlorophyll or colored organic matter) using remote sensing techniques may also be assessed on other estuarine and coastal systems, regardless of the main drivers of dynamics. Large coastal systems may resort to low spatial and high temporal resolution satellite sensors such as MODIS or VIIRS to study short-term or event-specific variability (e.g., reference [54]), or even to geostationary satellite sensors (e.g., reference [10]). Studies of smaller coastal water bodies still rely on a high spatial resolution characterized by a lower temporal resolution such as Sentinel-2 or Landsat-5, -7, -8 (e.g., the present study, [24]). Sensors with higher temporal resolution may be exploited to observe the state of the estuary during extreme events, or to make a comparison of the situation directly before and after the occurrence of events, as long as atmospheric conditions allow. In addition, the combined application of multiple satellite sensors could be more effective and fill gaps in our current observation of estuarine dynamics.
Our synergetic use of satellite-derived SPM and in situ data of river discharge and wind show exactly where the drivers of SPM variability have the most effect within the estuary. Satellite data providing unique spatial observations could then be complemented with moored stations at strategic locations (e.g., derived from the correlation maps of SPM and their drivers) equipped with optical (e.g., turbidimeters) and/or acoustic (e.g., Acoustic Doppler Current Profiler-ADCP) calibrated instruments to represent the highest temporal and continuous data sampling available-regardless of weather or sea state (e.g., references [51,55]).
For instance, based on our results and the successful examples of high-resolution turbidity sensors (e.g., the MAGEST program in the Gironde estuary (https://magest. oasu.u-bordeaux.fr/, accessed on 21 December 2022) [51,56]), we reason that the two best suitable locations for those sensors would likely be site I, West of the estuary mouth, where both spatial correlations of river discharge and wind speed were stronger and which is likely a region to be most affected by windstorms, and site II at the upper-most section of the estuary, where no spatial correlation was statistically significant, a possibility due to the almost permanent effect of the TMZ in the area (see location in Figure 1 for the suggested sites I and II, respectively). We speculate that with a data acquisition on the order of minutes, those high temporal resolution instruments, in strategic locations (sites I and II), are most likely able to capture the complex SPM variability that in situ data or satellite sensors are only partially able to explain. Additional opportunities include providing valuable insights towards assessing sediment budget/flow, migration, and expansion of the TMZ, or feeding ecological and hydrodynamical models.

Conclusions
This study contributes to the understanding of the spatial relationship between SPM and river discharge and wind speed. From the results of this study, we conclude that (1) the duration of extreme river discharge significantly decreased in the past 37 years (1984-2020), (2) satellite data at low temporal resolution is likely unable to capture extreme wind or river discharge events, while in situ SPM data are not sampled during extreme wind events for safety reasons, (3) SPM derived from the red band (or wavelengths < 800 nm) is less optimal to capture the peaks in SPM concentration that are expected during extreme events, and (4) synoptic satellite (Landsat) time series can be used in combination with time series of drivers, such as river discharge and wind, to capture where these forcing mechanisms affect SPM within an estuary as a tool to establish best locations for moored sensor stations.
Additional temporal point data of moored turbidimeters at strategic locations can complement these insights by providing data under particularly extreme events. Such a synergistic approach is valid both for retrieving and attributing the spatial and long-term temporal variability in SPM in estuaries around the world.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs15030670/s1. Table S1: Overview of in situ data. Table S2: Coefficients A and C of SPM algorithm. Original narrow-band and band-weighted coefficients (A and C; Equation (1) in main text) in the used SPM algorithm for each satellite sensor. Acronyms for satellite sensors are: TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper), and OLI (Operational Land Imager). Table S3: Match-up information data. In situ data is provided by RWS and by NIOZ. Refer to main text for full SPM data description. L5 is Landsat-5, L7 is Landsat-7, and L8 is Landsat-8. Figure S1: Interannual variability of SPM from in situ RWS datasets (grey scatter) and satellite-derived SPM (yellow scatter) relative to each region delimited in Figure 1 (main text). Overall yearly median is represented by black line. Figure S2A: Distribution of total number of satellite scenes and in situ SPM per classes of boreal season. Figure S2B: Distribution of total number of in-situ SPM data (from in situ data RWS) per classes of boreal season (a), (c), (e), (g), daily river discharge and (b), (d), (f), (h), daily wind speed. Shared areas mark the threshold for extreme events. Figure S3: Relationship between wind direction and SPM. Panels (a), (c), (e), (g) represent the median satellite-derived SPM in the delimited region (refer to Figure 1  Funding: This research received no external funding.

Data Availability Statement:
The datasets generated/analyzed for this study will be found in the DANS (Data Archiving and Networked Services) repository. Additionally, data processing scripts can be made available upon request to the corresponding author.