Conceição Lagoon is a subtropical coastal lagoon located on Santa Catarina Island, southern Brazil (Figure 1). The water body is a popular natural attraction in the region, widely used by the local and tourist populations for multiple purposes: fishing, transportation, recreation, and sports (Martini et al., 2006). The lagoon has a surface area of about 20 km² and a heterogeneous bathymetry, with extensive sandbanks of less than 1 m water depth contrasting with muddy sediment regions of up to 8.7 m water depth (Dias et al., 2014). According to the morphological features and water properties, the lagoon is traditionally divided into three subsystems: North Lagoon, Central Lagoon and South Lagoon (Martini et al., 2006; Lisboa et al., 2008; Lüchmann et al., 2008; Fontes et al., 2010). Its only connection to the Atlantic Ocean (a 2.8 km narrow channel in the Central Lagoon) dissipates approximately 84% of the tide effects, especially the short frequency astronomic component (99% attenuation) (Godoy et al., 2008). The lagoon’s exchange with the ocean is thus limited, with hydrodynamics studies suggesting water residence times from 150 to 270 days (Silva, 2021). As such, water motion is controlled primarily by wind and hydrological conditions in the contributing watershed (i.e., precipitation, evaporation, runoff) (Silva et al., 2017). The urbanized area in the watershed increased from 13% to 23% between 1990 and 2019 (Odreski et al., 2021). The sanitation infrastructure, nonetheless, was not developed at the same rate. It serves only 57% of the population and there are frequent issues with pumping stations’ overflows and irregular sewer connections (Odreski et al., 2021). Both the deficient sewage system and non-sewered areas contribute to the sustained organic matter and nutrient loading into the lagoon (Cabral et al., 2019). This has led to the chronic eutrophication of Conceição Lagoon, with recurrent algal blooms and hypoxia/anoxia conditions (de Barros et al., 2017; Silva et al., 2017; Cabral et al., 2019). The rupture of the disposal pond of a Wastewater Treatment Plant (WWTP) in late January 2021, discharging 79,000 m³ of wastewater into Conceição Lagoon (CASAN, 2021), caused severe consequences to the lagoon’s ecosystem and evidenced the urgency of monitoring wastewater spills in the area. Due to its magnitude and the presence of a visible plume in the following days, this event serves as a suitable training sample to investigate signals that could be used to identify wastewater from space.

The 2SeaColor ⁴ (Salama and Verhoef, 2015) is an analytical model that retrieves Inherent Optical Properties (IOPs) from remote sensing reflectance ( R r s ). The retrieved IOPs reflect the composition and state of the water (IOCCG, 2006) and therefore represent Water Quality Indicators (WQIs). The model simulates light interactions within the water by analytically solving the two-stream radiative transfer equations (forward model), including the three radiation components: downwelling direct and diffuse fluxes and upwelling diffuse flux. To reduce the dimensionality of unknowns, parametrizations are implemented as in Yu et al. (2016). An inversion scheme, based on spectral optimization with non-linear least-squares fitting, is used to derive water surface IOPs. Held as important factors in non-linear optimization, the initialization values are defined based on Lee et al. (1999). The equations of the forward model, parametrizations and initialization values are detailed in the aforementioned references and summarized for convenience in Supplementary Table 1.

In this study, we adapted 2SeaColor to be applicable to the Sentinel-2 MSI images for deriving three WQIs: i). a c h l a 440 —absorption coefficient of chlorophyll-a, ii). a d g 440 —combined absorption coefficient of detritus and coloured dissolved organic matter (CDOM), and iii). b b s p m 440 —backscattering coefficient of suspended particulate matter (SPM). These indicators were selected considering that wastewater plumes have been previously associated with high levels of CDOM (Marmorino et al., 2010; Ayad et al., 2020) and SPM (Ayad et al., 2020). High chlorophyll-a concentration has also been observed in the wake of wastewater outfalls (Trinh et al., 2017).

The 2SeaColor has been shown to perform well in both freshwater (Salama and Verhoef, 2015) and coastal turbid waters (Arabi et al., 2016; Yu et al., 2016; Arabi et al., 2018). Since the proposed method for formulating the WCI relies on the WQIs anomalies and their spatial and temporal variability (see section 2.4), validation of IOP retrievals specifically for the study area was not considered a requisite. Areas with water depths shallower than 1.5 m, based on the bathymetric map provided by Horn et al. (2022), were masked to avoid the bottom reflectance effect on the retrievals. The resulting time series of WQIs’ images was arranged as a raster cube, with the third dimension representing time.

The steps taken for the formulation of the WCI from the WQIs raster cube are summarized in Figure 3 and the processes are detailed in the sub-sections below.

To evaluate the WCI capabilities in detecting wastewater plumes in Conceição Lagoon, the signals produced on images associated with known outfalls that occurred in the area between 2019 and 2021 were investigated. Although at least six incidents have been reported in that period, most of them did not coincide with satellite overpass or there was cloud cover preventing analysis. For this study, two reported wastewater outfall events were analysed by inspection of the WCI maps at the locations of the expected plume. Environmental conditions (i.e., wind, precipitation, water quality measurements) and photographic records were also considered in the analysis.

In Figure 5 we have the mean monthly WQI anomalies for the three subsystems of the lagoon (North, Central, and South), in addition to the local temperature and precipitation climatology. It is not possible to observe a clear seasonality pattern, although there are indications that high positive anomalies of all WQIs coincide with an increase in precipitation at the beginning of (austral) Spring in September. Still, since in this study we are analysing events with scenes from February and May, potential effects of seasonality would not cause interference on the interpretation.

Figure 6 shows the WQIs anomalies of the training sample, which was used as input to the spatial variation mode analysis. It is possible to observe that a c h l a has homogeneous positive anomalies in the entire lagoon, with little spatial variability. a d g , on the other hand, has extremely high positive anomalies throughout the Central and North subsystems, while values indicate average conditions in the South Lagoon. For b b s p m we see high positive anomalies originating from the area of the disposal pond burst, with hotspots in the Central Lagoon and on the margins of the North Lagoon. Reflecting the spatial variability of the WQIs anomalies, the first Principal Component ( P C 1 ) of the training sample, which explained 56.6% of its variability, was given as: P C 1 = 0.02   . a n o m _   a c h l a + 0.92   .   a n o m _ a d g + 0.92   .   a n o m _ b b s p m (7)

The comparison between the first eigenvectors of the reference (training sample) and the other 137 images showed angle differences ranging from 17˚ to 178˚, with a mean of 54˚. This indicates that the reference eigenvector provides a distinct rotation when compared to the other eigenvectors, and confirms the atypical conditions seen in the training sample.

The reference P C 1 rotation implies that the anomalies of SPM backscattering and absorption by combined detritus and CDOM are equally important in the differentiation of the wastewater plume (49.5% each), whereas absorption by chlorophyll-a provides a very low contribution (1%). This suggests that the wastewater plume in this case is characterized by simultaneous positive anomalies of a d g and b b s p m . For that reason, the optimal Linear Combination of WQIs anomalies (LC) for wastewater detection was defined as: L C = 0.50   .   a n o m _ a d g + 0.50   .   a n o m _ b b s p m (8)

Figure 7 shows the empirical cumulative distribution function (ECDF) of the LC raster cube, from which the minimum (1% probability) and maximum (99% probability) values for rescaling were defined as −1.35 and 2.28 respectively.

In Figure 8 we have the comparison between the empirical CDFs of the E. coli measurements and the WCI values. The E. coli CDF shows that 58% of the measurements fall within the threshold of low risk of contamination (<200 MPN/100 mL) and 80% within the medium risk upper threshold. Transferring these percentiles to the WCI CDF, the thresholds for the risk of contamination are defined as follows:

• Low: WCI < 0.41

As an indirect verification, we checked the distribution of the risk classes based on the WCI and E. coli measurements for each of the sampling sites during the studied period (Figure 9). For most of the sites, the WCI was able to capture very similar distribution features as the indicated by E. coli measurements. Exceptions are seen for points 38 and 62, where especially the high-risk class frequency was underestimated by the WCI. For point 43 the WCI overestimates the medium and high risk classes, and underestimates the low risk. It should be noted that these three points are located at the transition between the South and Central subsystems of the lagoon (see Figure 2), which per se has complex dynamics. In addition, the morphology of this region makes it highly susceptible to adjacency effects.

The spatial variation mode analysis of the training sample provided a linear combination that best captured the reference wastewater plume with equally high weights for a n o m _ b b s p m and a n o m _ a d g (50%), and no contribution from a n o m _ a c h l a . Note that during PCA the variables were standardized to avoid bias towards the one with higher overall values and variance.

The high weight obtained for a n o m _ b b s p m is physically sound considering that, given the nature of the reference plume (i.e., wastewater-slurry mixture from disposal pond flash-flood), elevated concentrations of suspended solids are expected. This is in agreement with Ayad et al. (2020), who showed that mixed stormwater and wastewater had medium to high turbidity levels associated with suspended sediments. As for a d g , wastewater is known to be closely associated with dissolved organic matter (Hudson et al., 2007), which makes it logical for its anomaly to have an important role in distinguishing the plume. This is also in line with the results of Marmorino et al. (2010), who found sewage discharge to be related to high levels of CDOM.

The negligible effect of a n o m _ a c h l a to distinguish this plume is in accordance with other studies that were not able to identify a detectable satellite-derived feature of chlorophyll-a in wastewater contaminated areas (Gierach et al., 2017; Seegers et al., 2017). Possible reasons given by the authors include the short duration of the outfall, chlorination in the case of treated wastewater, or strong mixing and dilution. We also argue that phytoplankton growth due to a wastewater spill would likely be a post-effect following nutrient enrichment and dependent on other environmental conditions, i.e., temperature, stratification and light availability.

A major assumption in this study is that the signals associated with the reference wastewater plume (derived from one specific outfall event) would be representative of other wastewater plumes in the area. This was a large-scale event, with an extensive visible plume even after weeks of the rupture. Therefore, it constituted a useful training sample to investigate the satellite-derived features associated with wastewater. Nevertheless, the weights given to the WQIs anomalies in the linear combination to formulate the WCI carry embedded characteristics of this one event, such as the large influence of suspended sediments due to the flash-flood nature, in combination with organic material. We would therefore expect that plumes that resemble the reference would be more easily detected and classified with a higher risk of contamination. In contrast, when only a d g is positively anomalous, which could in many cases also indicate wastewater contamination, the WCI may not capture the risk level well enough (false negative alarm). It is also possible that b b s p m is anomalously high due to strong resuspension related to winds or extreme events, and the WCI captures a high risk of contamination that does not reflect the truth (false positive alarm). It should be noted that, the index being used as an operational tool, it is always possible to consult the WQIs anomalies maps individually (as these would have to be generated in the process), which could help assessing situation.

Examples of potential underestimation of contamination risk were seen in the WCI maps of 05-02-2021 (Figure 10B), 11-05-2020 and 16-05-2020 (Figure 11). All of these maps were associated with wastewater spill events, but the coefficient of backscattering by SPM was low in comparison with average values. Therefore, even with markedly positive anomalies of a d g dominating the light attenuation, the risk was still interpreted as medium-low (0.3–0.5) as a cancellation effect of negative b b s p m anomalies. It should be noted that CDOM is more conservative than SPM, with SPM decreasing rapidly due to settling and CDOM decaying over a time frame of weeks to months with photodegradation (Nezlin et al., 2008). The WCI does not take this into consideration at this point.

In addition, overestimations in SPM backscattering could also result in misinterpretation of the contamination risk. Overestimations of b b s p m in the 2SeaColor model are caused mainly by unmasked glint, clouds, bottom reflectance, foam, or boats, which all cause high NIR reflectance. These overestimations could be minimized by having a more efficient cloud masking algorithm for coastal and inland waters, such as WiPE (Ngoc et al., 2019), the inclusion of bottom reflectance in the model so that shallow waters can also be evaluated, and a routine to detect boats or other artificial objects.

Ideally, a larger number of images portraying known wastewater plumes in varied conditions should be used to have more representativeness and statistical strength in the weights of the linear combination and formulation of the WCI. It would also be preferred to have more water quality in situ measurements coinciding with satellite overpass on days of known contamination. Still, challenges arise from the fact that information about wastewater spills is not commonly disclosed by sanitation companies, especially in developing countries, as this might have negative implications for them. In addition, overflows from sewage pumping stations are typically associated with heavy rainfall and cloud cover. This constitutes a constraint of the method since optical remote sensing relies on clear skies. Yet, depending on the scale, duration and circulation conditions during the outfall, plumes might still be detectable after a few days, as was seen in the cases analysed in this study.

Ayad et al. (2020) highlight the importance of considering circulation patterns in wastewater studies, as there is a knowledge gap regarding residence time and dilution of the plumes. According to the study of Silva (2013), Conceição Lagoon’s circulation patterns are governed primarily by wind, flow in the contributing streams and the meteorological tide. The wind was the only parameter available for the study period and considered in the analysis, which highlighted its influence on the transport of wastewater plumes and the overall dynamicity of Conceição Lagoon system. Nevertheless, the inclusion of flow and meteorological tide data could enhance the robustness of the analysis by allowing the isolation of the processes that contribute to the detection of contamination.

The methodology applied here for the construction of the WCI is easily transferable to other locations and it is desired for further verification. For that, it is only necessary to have a time series of WQIs, with at least two dates capturing a wastewater plume (one for training and another for testing). It should be noted that the span of the time series influences the calculation of the WQIs anomalies (and consequently the WCI), as the mean and standard deviation will depend on that. It is also important to investigate the seasonality of the WQIs in the area of study and, if present, the anomalies should be calculated based on seasonal values. This is because if the WQIs have distinct seasonality, calculating the mean and standard deviation over the whole period will not remove it, and seasonal characteristics will still be present in the standardized anomalies, biasing the interpretation.

The choice of the training sample(s) is also important. It should consist of scenes that capture the spatial dynamics of wastewater plumes that are considered representative of the study area. A higher number of training samples should enhance the representativeness and is therefore preferred. The coverage of the training sample should comprise the whole area of interest in which the WCI will be applied, as to consider the site-specific characteristics in the differentiation of the background waters and the plume. The optical complexity of the area and the characteristics of the wastewater spill will influence the weights given to the WQIs anomalies in the spatial variation mode analysis, and provide an optimal linear combination to indicate the risk of contamination that is unique to the location. Although the weights are expected to vary from location to location, the overall methodology is transferable and could constitute a useful management tool for water management.