COMBINED FLOODING AND WATER QUALITY MONITORING DURING SHORT EXTREME EVENTS USING SENTINEL 2: THE CASE STUDY OF GLORIA STORM IN EBRO DELTA

: Short extreme events have significant impact on landscape and ecosystems in low-lying and exposed areas such as deltaic systems. In this context, this paper proposes a combined methodology for the mapping and monitoring of the flooding and water quality dynamics of coastal areas under extreme storms from Sentinel 2 imagery. The proposed methodology has been applied in a coastal bay of the Ebro Delta (Catalonia, NE Spain) to evaluate jointly the impact of Gloria storm (January 2020) in land-flooding and water quality. The experimental results show that the Gloria storm had a strong morphological impact and altered the water quality (chl-a) dynamics. The results show a recovery in terms of water quality after some weeks but in contrast the coastal morphology did not show the same degree of resilience. This paper is the first step of an overall goal that is to set the bases in a long term, for a workflow for rapid response and continuous monitoring of storm effects in coastal areas and/or highly valuable ecosystems such as the Ebro Delta.


INTRODUCTION
The nature of coastal areas makes them highly dynamic and susceptible to extreme variations under anomalous weather circumstances. Particularly, low-lying areas such as deltaic systems are highly vulnerable to extreme storms. The impact of the storms, and associated hazards (e.g. wave-driven, sea-level rise hazards) induces morphodynamic processes such as beach erosion or overwash  and flooding episodes that may affect large extensions where multiple human uses coexist (e.g. urban areas, agri-food industry, ecosystem services). In addition to the physical impacts in the terrestrial ecosystem, significant rainfall and strong winds related to extreme storms increase run-off and turbulence in coastal waters, thus with an expected impact on water quality and aquatic ecosystems (Hernández et al., 2020). A better understanding of the relation between storm events and landscape/ecosystem development is needed. This is a key step for improving territorial and emergency management strategies on storm-induced hazards assessment, preparedness, response, and environmental resilience monitoring.
The Ebro Delta is an excellent case of study because supports important terrestrial and aquatic ecosystems, as well as economic activities sensitive to extreme storms. It is a subsident area with an unprotected coastline enclosing two coastal bays with high ecological and socio-economical value. The importance of storm hazards in the Ebro Delta in terms of magnitude and recurrency has been increasing since the 90s with most hazardous storms characterized by low-pressure systems off the Ebro Delta coast and eastern wave storms, usually accompanied by storm surges and heavy rainfall. Such events lead to a compound flooding that * Corresponding author affects most of the Ebro Delta surface and induces morphological disturbance in the coastal front and bays (Jiménez et al., 2011). Furthermore, increased discharge of the Ebro River, land run-off, and sediment suspension and transport in coastal waters suppose a direct impact on the aquatic ecosystems and related economic activities such as fisheries or aquaculture (Berdalet et al., 2020).
This, together with the expected exacerbation of these trends under a climate change scenario has led to several studies assessing the effect of these events in the area. Some studies have been focused on reporting and documenting the events with social media and media pictures, visual interpretation of airborne imagery, and reporting of the consequences (Jiménez et al., 2020, González et al., 2020. Alternatively, most scientific research efforts have been focused on the investigation of flooding risk and coverage, and morphodynamics through numerical model simulations under storm and climate change scenarios, and do not include reports on water quality dynamics (Amores et al., 2020, Garcia et al., 2013, Alvarado-Aguilar et al., 2012. All this information is a highly valuable source of knowledge for assessing the consequences of storm-induced hazards in the delta contributing to improving risk, land, and emergency management strategies. However, there is still a gap for improvement in terms of observation-driven methods for the monitoring of storminduced hazards, complementing and/or enhancing the current approaches in the area. Mainly, on the assessment of flooding events and associated morphodynamic changes under highintensity storms over large areas, and the monitoring of water quality parameters (e.g. Total Suspended Solids (TSS) or chlorophyll-a concentration (chl-a)), when access to in situ measurements may be limited. From this perspective, spaceborne multispectral remote sensing is a unique source of data as it provides a synoptic view over large areas and frequent systematic observations. Among remote sensing platforms, the Sentinel-2 multispectral imagery satellites (S2) operated by the European Spatial Agency (ESA) represent a breakthrough in terms of flooding and water quality monitoring. The S2 MultiSpectral Imagery (MSI) covers the VIS-NIR and SWIR regions of the spectrum at 10, 20, or 60 m spatial resolution (depending on the spectral band) with a frequency of 5 days, improving temporal and spatial resolutions of prior similar platforms such as Landsat. Multispectral data or optical imagery can be used for both flood and water quality mapping (Caballero et al., 2020, Sobel et al., 2020, Goffi et al., 2020, Cavallo et al., 2021.
Detection of flooding and land surface water change based on multispectral satellite imagery generally relies on the use of water-related indicators proposed as combinations in the VIS, NIR, or SWIR spectral regions in form of spectral index (Pekel et al., 2014). Overall, combinations of VIS/SWIR and NIR/SWIR spectral bands provide for a clear diagnostic of water surfaces, showing positive values for flood conditions and negative values for soil (both wet and dry) (Boschetti et al., 2014). However, as the selected index response varies as a function of the study area, water characteristics, environmental and atmospheric conditions, the combination of multi-source data is recommended for improving classification accuracy (Acharya et al., 2018). Multiparameter classification methods aim to reinforce the characterization by exploiting redundancy or complementarity provided by multiple spectral features (Goffi et al., 2020). Therefore, the integration of different band combinations has been reported to improve the accuracy of land cover mapping (including flooding and classification of dry and mixed surfaces) through different segmentation techniques (e.g. Cavallo et al., 2021, Tavus et al., 2020. Extensive research can be found related to the retrieval of water quality parameters such as TSS or chl-a from remote sensing reflectances (Rrs). The common approaches rely on empirical and/or bio-optical algorithms that exploit spectral signatures of different water constituents, along with different band combinations depending on the water optical properties (Moore, et al., 2014, Caballero, et al., 2020. All of them use a set of insitu values to calibrate and validate the different models. The performance of the different algorithms depends mainly on the type of water and the accuracy of the atmospheric correction, but also on the dynamic range of concentrations of water quality parameters within the area of study. Some research has been done in relation to the monitoring of the impact of short extreme events on water quality. (Sobel et al., 2020) used S2 imagery to evaluate the impact of several Hurricanes in the Gulf of Mexico in terms of TSS concentration. Several studies suggest that water quality recovers relatively fast from severe events (Chen et al., 2017, Huang et al., 2011 while others show long-term impacts (Beaver et al., 2013, Wetz andYoskiwitz, 2013).
This study presents a methodology for the monitoring of the flooding and water quality (i.e. chl-a) dynamics from S2 imagery under extreme storms in coastal areas. In particular, the proposed methodology has been applied in the bays of the Ebro Delta (Catalonia, NE Spain) to jointly evaluate the impact of landflooding and water quality in Alfacs Bay during and after the Gloria storm (January 2020). The study aims to demonstrate how S2 imagery together with appropriate modelling can be used as a tool for assessing the ecosystem response of a coastal area in the aftermath of an extreme climatic event, through an integrated approach. Moreover, this study presents the first steps of an overall goal that is to set the bases in a long term, for a workflow for rapid response and continuous monitoring of storm emergencies in low deltaic and highly vulnerable areas such as the Ebro Delta.

Study area
The Ebro Delta (Catalonia, NE Spain), with a total extension of 32.500 ha is one of the largest wetlands in the Western Mediterranean region (Figure 1). It supports highly dynamic and productive ecosystems with great ecological and economical value. Several protected areas (e.g. RAMSAR sites, UNESCO biosphere reserve, the Natura 2000 network) coexist with a number of anthropic activities. Agriculture is the core of the economy with ~ 65 % devoted to rice production (Genua-Olmedo et al., 2016). Inside Ebro Delta bays, the production of bivalves constitutes also a major economic activity, although fisheries and exploitation of saltworks are also part of the diverse economic activities in the area ( Figure 1).

Figure 1. Summarized depiction of the study area
From a topographic point of view, Ebro Delta is a low-lying coastal area with sandy, silty and sandy, and fine soils. However, an abrupt decrease in Ebro River sediments discharge and the accumulative subsidence have reversed the growing trend of the Ebro Delta surface during the last century. The Ebro Delta currently presents a wave-dominated coast with strong reshaping processes (Sayol and Marcos, 2018), highly exposed to storms and extreme events, subjected to massive flooding, beach erosion, and overwash episodes (Jiménez et al., 2020). From the perspective of water quality, it highlights the importance of the semi-enclosed bays formed by the northern and southern hemideltas ( Figure 1). Both coastal bays have a relatively small surface, shallow waters (0-7 m), and narrow mouths connecting them to the Mediterranean Sea. The bays harbour important populations of seagrass and benthic organisms, as well as fisheries and shellfish farming activities, which are sensitive to water quality deterioration due to storm episodes.
The results presented in this paper are focused on Alfacs Bay (Figure 1), the southern coastal bay where many of the threats driven by storms in the Ebro Delta are observed. Alfacs bay concentrates aquaculture, an extremely vulnerable piece of coast heavily affected by erosion (Trabucador barrier), and La Banya spit, a special Protection Area for birds including a dune system that is part of Ebro Delta Natural Park.

Gloria storm
Gloria storm was a marine storm affecting the Ebro Delta region in January 2020 (19.01.20 -23.01.20). It has been classified as one of the most intense among the events in the region during the last decades with immediate economic losses of several millions of euros. During the event, the situation was dominated by a southern deep low-pressure system that generated sea level surface elevation (up to 1 m) and strong easterly winds which induced wind waves (Amores et.al, 2020) with high significant height (maximum Hs > 7m) and heavy rainfall on land (maximum cumulative precipitation > 130 mm). The Gloria storm had devasting effects on urban areas, agriculture fields, the natural environment, and different crucial infrastructures such as roads, becoming the most hazardous event in Ebro Delta in the last decades.
NDWI  (1) NDWI XU = 3 − 11 3 + 11 (2) These indexes were used later as input in an unsupervised kmeans clustering to generate 4 classes for each image. Then, the mean and standard deviation of pixels within each class and SI (4 classes × day -1 × number of dates; N = 44) were extracted and reclassified into two classes by means, again, of unsupervised kmeans clustering. The NDWIXU cluster centers were then used for defining water (positive NDWIXU ) and non-water classes (negative NDWIXU) resulting in binary maps. This method aimed to first differentiate 4 classes which may be water or non-water and then solve the uncertainty of single images classification by time-series analysis of different classes accounting for different types of water and surface state. It was implemented using a customized script developed for R software v3.6 (R Core Team, 2017).

Chl-a mapping
The chl-a maps during the Gloria storm are part of a longer time series (October 1 st 2019, to September 30 th 2021) in the area of Ebro Delta (Figure 1 top). This subsection summarizes the methodology used to generate this time series of chl-a maps.
All available S2 L1C orbit R051 images during the two-year period of interest were downloaded from the Copernicus Services Data Hub (https://cophub.copernicus.eu/dhus/#/home). After image-by-image visual checking for clouds and shadows over the Ebro Delta coastal waters, the valid images (Nimages = 210) were pre-processed. This step included the resampling of the bands to 20 m pixel, the subset to the extent of the area of interest, and the atmospheric correction specific for coastal waters, all implemented through the Graph Processing Tool of SNAP v8.0. The atmospheric correction was done with the latest version of the Case 2 Regional Coast Colour processor (Doerffer, et al., 2007, Brockman et al., 2016, here referred as to C2XC. It was applied on all valid S2 images, with a fine-tuning of daily auxiliary data (pressure, ozone, altitude, salinity, and water temperature). For land/water segmentation, the valid pixel expression was set as a threshold on band B11 of S2 L1C images. The threshold was defined independently for each image with the histogram-based triangle method. The output of C2XC included Rrs at VIS-NIR S2 bands.
Field measurements (N = 146) were conducted along 26 dates, at different locations including waters in Alfacs bay, Fangar bay, and the surrounding exterior Mediterranean sea. The measurements corresponded to integrated water column samples, taken within a window of ±3 hours of S2 pass. For chl-a determination, water samples were filtered through 0.4-0.6 µm GF/F glass fiber filters, extracted using standard methods (Shoaf and Lium, 1976), and calculated with Jeffrey and Humphrey (1975) equations.
Then, the match-up exercise between field-measured chl-a data and C2XC Rrs was done using a quality-check procedure based on 3x3 pixel-window centered at field measurements. The quality check included the removal of pixels flagged inherently by C2XC, an additional outlier filtering using boxplot analysis and, the removal of negative Rrs pixels. Only pixel windows with more than 4 pixels remaining were accepted for posterior analysis. After this step and the rejection of cloud-affected S2 images, the number of valid field chl-a data was 73, ranging between 0.1 -8 mg/m3, along 14 different dates.
For chl-a estimation, a set of 11 band combinations (2 to 4 spectral bands) in form of spectral indices (SI) were computed for all valid match-ups. These SI included bands B1 (centered at 443 nm) to B6 (centered at 740 nm) and exploit specific chl-a absorption peaks in Blue and Red spectral regions, and/or use Red Edge information for mitigating the effect of absorption by non-algal particles and yellow substances. For instance, common ocean colour algorithms based on Blue to Green ratios (O'Reilly and Werdell, 2019), the three-band model of Gitelson et al. (2011), and the Normalized Difference Chlorophyll Index (Mishra and Mishra, 2012) were included among the tested SI. Chl-a was modelled with all of them. A 70 % of chl-a data was used for model calibration (cal) and 30 % for validation (val). The cal/val datasets were generated randomly. Models were developed for raw and log-transformed (when possible) Rrs and chl-a data (i.e. raw-raw, log-raw, log-log). Linear, linear piecewise (1 breakpoint), polynomial (2 nd , 3 rd , and 4 th order), logarithmic, power, and exponential models were tested. All the process was repeated through 100 iterations, with varying cal/val datasets. For each model (Band Combination × Type of Fitting × Iteration), performance was evaluated by means of the following statistics: the Mean Average Error (MAE), the Root Mean Squared Error (RMSE), the Average Percentage Difference (APD), the BIAS (eq.2-5) and the Pearson's r (r). The statistics were computed including calibration and validation datasets.
To generate the chl-a maps, all flagged pixels (C2XC flagging) or pixels with negative Rrs in any of the spectral bands of each valid S2 image were removed. A mask was applied for avoiding most shallow waters and the presence of seagrass, which make the retrieval of chl-a more uncertain. Then, the model with the best performance was applied to all valid images and pixels. Those with negative chl-a value (residual) were set to 0 mg/m 3 and, pixels with chl-a > 30 mg/m 3 were removed (noData) given the higher uncertainty in so high concentrations. Pixels corresponding to the mussel rafts were masked out with a rafts' shapefile ( Figure 1).
Where is the estimated value and the observed (i.e. field measured) value.

Chl-a estimates
Overall, models that performed better were related to Rrs bands combinations including B2 (centered at 490 nm) or B3, and B4 (centered at 665 nm) or B5 (centered at 705 nm). The best model for chl-a estimation was found with a simple B5/B2 band ratio fitted with a 2 nd polynomial function (eq. 6).  (Figure 2), overestimation of retrieved chl-a concentrations in scenarios with great concentrations of TSS is expected. Empirical algorithms, as the one used, can perform well only inside their range of derivation and for the area for which they were derived. They are more limited in their ability to discriminate between non-unique signals from parameters that may be covariant with chl-a, for instance, TSS (Matthews et al. 2012), which increases particularly during storm events. This could not be validated because no chl-a nor TSS data were available in such cases, but it was observed in RGB composites that after strong winds or storm events, the water of the bays turns darker brownish or whitish, suggesting higher TSS (e.g. in Figure 3 and 4B). In these cases, the selected model will tend to overestimate chl-a since the reflectance on B5 increases with great TSS, overlapping the absorption of chl-a in the Blue (Soriano et al., 2019) In these cases, other empirical or semi-analytical algorithms may be more accurate. However, it is important to highlight that the proposed methodology for chl-a retrieval is easy to replicate, and straightforward using several software tools such as SNAP, which may be important for a wide part of the scientific community related with ecosystems monitoring and modelling response but with few or no remote sensing skills.
Nevertheless, further research should move towards a multialgorithm combination, suitable across different types of water optical properties. The development of a blending algorithm requires accurately defining the optical water type (Moore et al., 2014) accounting for the spectral shape, magnitude, and distinctive Rrs spectral features. Furthermore, a greater amount of field data is needed, covering all different scenarios, which were not represented in the current chl-a dataset.

Flooding mapping
Systematic validation of flooding maps was not conducted, albeit a visual comparison of the resulting maps showed good agreement with what can be interpreted from the 10 m RGB composites of S2 in Alfacs bay (e.g. in Figure 3). The method allowed us to assess the evolution of the main morphological phenomena that occurred during the Gloria storm in the bay, the breaching of the Trabucador barrier. From visual interpretation, the main limitations of the method arise from; i) small features (e.g. shellfish rafts, borders between dry and flooded surfaces) which may be misclassified due to the pixel mixing with the flooded surrounding environment; ii) increasing uncertainty in presence of emerged or submerged vegetation and continuous succession of small emerged and submerged patches (Ozesmi and Bauer, 2002).
Nevertheless, it is important to remark that the spectral distribution before and after an extreme even could be quite different. Additional research is needed for further validating the accuracy of the proposed methodology. For instance, the minimum number of images to be used must be defined. Also, the 4 initial k-means classes may be a constraint in more heterogeneous environments, needing to explore the increasing of clusters number. Regarding the improvement in the detection of small features, suitability for resampling to 10 m instead of 20 m pixel can be studied, albeit it will strongly depend on the native resolution of S2 bands. In addition, information on the sea level simultaneously to the S2 pass must be considered in deeper analysis. Although that, the proposed method does not depend on large time series, avoids the difficulty of finding stable thresholds as in decision trees or rule-based approaches, is fast and unsupervised, and requires simple computation, allowing for rapid mapping in the aftermath of a storm, which is a key pillar in the management of extreme events.

Chl-a and morphology response to Gloria storm in Alfacs bay
For assessing the dynamics of water quality and flooding in the Alfacs Bay, three status maps ( Figure 4) and three 'change maps ( Figure 5) were generated within a period of one month encompassing the Gloria storm (21.01.20).
The Gloria Storm caused breaching of the Trabucador barrier, without showing signs of recovery after 20 days ( Figure 4B). The thin line that seems to reconnect the delta with La Banya spit ( Figure 4C and 5B) was the ongoing reconstruction of the road for in/out of the trucks transporting salt from the saltworks (Figure 1). La Banya spit remained disconnected from the delta during all the period shown ( Figure 4C and 5C), preventing the passage of vehicles and the transport of goods to and from the saltworks. It can be observed in the change maps that the Trabucador barrier breached from the center and receded back into the bay. The setback was evident in the extremes of the barrier as observed in Figures 5A and 5C corresponding to the difference between the 15 th and 20 th days after the storm and the image of 5 days before the Gloria storm.
In La Banya spit (Figure 1), the most important effects of Gloria storm were observed on the eastern side, at the dune fields. In this area, the barchan dunes in the shoreline side showed stability, with no flooding detected ( Figure 5). The inner dunes, however, showed a more dynamic behavior, alternating flooding, and drying of the northern and southern faces during the period of study ( Figure 5A-B). Where flooding is more intense, the inner parts of La Banya spit were subjected to loss of volume (Rodríguez-Santalla et al., 2021). The drying out of some of these parts may be explained by the accretion due to sand blown from the beach and dunes and trapped by vegetation growing on the foredune. This phenomenon is observed in Figure 5C, with most of the dried area concentrated behind barchan dunes. This highlights the importance of the barchan dunes for maintaining the shoreline during a storm (Rodríguez-Santalla et al., 2021). Regarding water quality dynamics associated with the Gloria storm, water turbulence increased, particularly in the inner north-eastern Alfacs Bay, related to the Trabucador breaching (Fig. 4B) and the shallow shelf formed behind it. This led to increased chla concentrations, although the model could be overestimating them in this scenario, as high TSS concentration is expected in this area behind the Trabucador barrier breaching stretch. In the last date shown (10.01.20 - Figure 4C), chl-a concentrations were again closer to pre-storm concentrations ( Figure 5C). Different from the Trabucador barrier, the water quality showed signs of recovery in a short term. In the long-term, water quality dynamics should be revised again since the breaching of the Trabucador barrier has important implications in the morphodynamics and hydrodynamics of the bay (Gracia et al., 2013). These are major factors for the development of phytoplankton populations related to the modulation of chl-a inside Alfacs Bay, with a direct effect on the farming of shellfish (Soriano et al., 2019).

CONCLUSIONS
This study presents a methodology for both chl-a (in water) and flooding mapping (on land) based on parallel processing of Sentinel 2 multispectral imagery. The methodology has been applied to the assessment of impacts of an extreme coastal storm affecting the Ebro Delta in January 2020. The Ebro Delta is a suitable place to apply the methodology developed here since it shares common threats with other coastal regions and, more specifically, with other deltas.
Despite constraints related to the validation of flooding maps and the uncertainty of chl-a retrieval under optical water types out of the range of development of the model, the results obtained are promising. The model presented for chl-a mapping showed an accuracy greater than 70 % and, the methodology proposed for mapping of land water surface changes demonstrated great agreement with visual interpretation from 10 m RGB imagery. The exposed methodological development allowed assessing changes both in water and land, through a well-explained processing chain which can be easily implemented almost with SNAP open software. Both, flooding delineation and water quality mapping methodologies can be applied fast, allowing for rapid evaluation of storm-related perturbances on terrestrial and aquatic ecosystems, with particular importance in the management of coastal areas. Sentinel 2 demonstrated great capabilities for complementing current methods and supporting management in diverse environments, through the potential for conducting studies with integrated approaches, which is currently a line to follow for achieving the sustainable development goals.
The research presented in this paper is one of the first steps with Sentinel 2 of a long-term strategy to understand the ecosystem response in Ebro Delta both in space and time in front of short extreme events. Further research should include an improvement of the current water quality modelling aiming to de-coupling TSS from chl-a to improve chl-a estimation. The flooding mapping must be validated, assessing the accuracy of the method, and arising potential constraints or improvements not included in this preliminary study. The proposed methodology is planned to be extended to all Ebro Delta region and tested also in other short extreme events that occurred recently such as Filomena storm (January 2021).