Unmixing water and mud: Characterizing diffuse boundaries of subtidal mud banks from individual satellite observations

Mapping of subtidal banks in mud-dominated coastal systems is crucial as they influence not only shoreline and ecosystem dynamics but also economic activities and livelihoods of local communities. Due to associated spatiotemporal variations in suspended particulate matter concentrations, subtidal mudbanks are often confined by diffuse and rapidly changing boundaries. To avoid inaccurate representations of these mudbanks in remote sensing images, it is necessary to unmix distinctive reflectance signals into representative landcover fractions. Yet, extracting mud fractions, in order to characterize such diffuse boundaries, is challenging because of the spectral similarity between subtidal-and intertidal features. Here we show that an unsupervised decision tree, used to derive spatially explicit and spectrally coherent image endmembers, facilitates robust linear spectral unmixing on an image-to-image basis, enabling the separation of these coastal features. We found that resulting abundance maps represent cross-shore gradients of vegetation, water and mud fractions present at the coast of Suriname. Furthermore, we confirmed that it is possible to separate land, water and an initial estimate of intertidal zones on individual images. Thus, spectral signatures of end-member candidates, determined from relevant index histograms within these initial estimates, are consistent. These results demonstrate that spectral information from well-defined spatial neighbourhoods facilitates the detection of diffuse boundaries of mud-banks with a spectral unmixing approach.


Introduction
Accurately mapping boundaries of geomorphological landforms in mud-dominated coasts is crucial for integrated coastal management, and for understanding landscape evolution in these dynamic areas (Koohafkan and Gibson, 2018).These landforms are often associated with the distribution of large sediment concentrations, originating from highdischarge rivers such as the Amazon, Mississippi and Yangtze rivers (Murray et al., 2019).The coalesce of mud is typically favored when coastlines are confined, wave energy is low and tidal range is large (Anthony et al., 2013), resulting in characteristic landforms and coastline fringing wetland vegetation such as mangroves and saltmarshes.In some settings the presence of fluid mud can also trigger the formation of subtidal mudbanks along open coasts (Anthony et al., 2010).
These banks are associated with extreme spatiotemporal variations in suspended particulate matter (SPM) concentrations in the water column (Gratiot and Anthony, 2016).This can be related to their non-linear response to for example, wave climate, tidal stage, currents, proximity to the coastline and liquefaction processes (Vantrepotte et al., 2011).As a result, the boundaries between mudbanks and adjacent heterogeneous coastal features are inherently diffuse and change rapidly.Sediment exchange and increased wave damping potential associated with mudbanks, provide a window of opportunity for mangrove species to colonize large intertidal surfaces (Balke et al., 2011).These mangroves provide regulating services, such as protection against sea level rise and carbon sequestration (Anthony et al., 2013).Hence it is vital to monitor mud-bank dynamics, to support adaptive coastal management-and conservation strategies that account for material exchange between these banks and the coastline.
Due to the inaccessibility of mudbanks and the spatial and temporal scales involved in their dynamics, monitoring via remote sensing emerged as a cost effective alternative to field surveys (Augustinus, 1980;Vantrepotte et al., 2011) and aerial photographs (Augustinus et al., 1989).Various methods, including SPM inversion algorithms (Froidefond et al., 2004;Zorrilla et al., 2018), as well as supervised-and unsupervised classifications (Anthony et al., 2008), have been developed to estimate mud-bank characteristics.Some of these methods have been applied over increasingly large areas and improved the temporal resolution, especially since the introduction of data-cubes such as Google Earth Engine (GEE) (Gorelick et al., 2017).
At the same time, challenges emerged in the identification of these mudbanks from time-series of satellite observations due to missing data and uneven sampling across tidal stages.This is related to the sunsynchronous orbits of satellite constellations, such as those from Landsat, that never capture the extremes of low-and high-tide with favourable cloud cover (Murray et al., 2019).Previous attempts successfully focussed on extracting spatially coherent information of intertidal features by selecting (Murray et al., 2012) and aggregating observations for specific tidal stages (Sagar et al., 2017).Yet, for subtidal features, resuspension of mud and migration processes at seasonal timescales are responsible for the spatiotemporal variability of their footprints (Zorrilla et al., 2018).This suggests that pixel based image compositing, where one tries to overcome missing data by reducing multiple observations to a single 'best' observation, potentially misrepresents diffuse boundaries and temporal evolution of these features (Koohafkan and Gibson, 2018).This is related to the fact that reflectance signals describe subtle differences in composition, resulting from radiation interacting with multiple active substances in a pixel (Odermatt et al., 2012).Consequently, the spectral similarity in muddy coastal systems, together with differences in grain sizes, turbidity and soil moisture content, complicates differentiation between subtidal-and intertidal features (Ryu et al., 2002).Accordingly, image analysis techniques that involve semiautomatic unmixing of landcover fractions from distinctive reflectance signals, offer an approach for the analysis of gradients in muddominated coastal system that are associated with diffuse mud bank boundaries (Alcântara et al., 2009).Where per-pixel classifiers assign each pixel to a class based on its similarities, unmixing methods model abundance as a linear-or nonlinear combination of each provided spectral signature.This means that the presence of more materials within one pixel is estimated, based on the provided end-member signatures (Shanmugam et al., 2006).Especially linear spectral unmixing (LSU) is a convenient unmixing tool to handle mixed pixels, as it does neither require extensive training data nor a computationally demanding analysis (Somers et al., 2011).
The aim of this study was to develop a data-driven approach to analyse diffuse boundaries of mudbanks from individual Landsat images using LSU and cloud computing in GEE.The method is evaluated on its ability to consistently select end-member candidates for the purpose of unmixing landcover fractions of water, mud and vegetation, using a section of the Suriname coast as case study.

Study area
The study area is located near Paramaribo, the capital of Suriname (Fig. 1).This area was selected because it is part of the Guyana coastline, a prime example of a mud dominated coast (Augustinus, 1980).About 15-20% of the sediment migrates (0.5-5 km a year) alongshore as mud banks attached to the coastline, which are continuously reworked by waves and currents (Anthony et al., 2010).Due to the dynamics of these mud banks, the topography and bathymetry of the Suriname coastline experiences quasi-cyclic variations in erosion and progradation, related to inter-bank and bank phases that last up to 30 years (Allison and Lee, 2004).Trade winds and precipitation vary on a seasonal scale, with more south-easterly trade winds (3-9 m/s) during the major dry season (August-November) and more easterly to north-easterly winds during the major wet season (April -August).This change in wind direction results in higher swell waves in the wet period.The tide is semidiurnal with a range up to 2.5 m during spring tide.As a result of these environmental conditions, the migrating mudbanks continuously change in shape and orientation (Augustinus, 2004;Augustinus et al., 1989), adding to the diffuse character of their boundaries (Vantrepotte et al., 2011).

Selecting end-member candidates
The coastal system of Suriname is characterized by spatially heterogeneous mixtures of vegetation, mud and water.In order to define a collection of spectral signatures that represents these distinct substances, endmembers are identified.Because endmembers from spectral libraries and field surveys do not handle variability between image acquisition very well, a more suitable approach was chosen by deriving these spectra directly from end-member pixels in each image (Somers et al., 2011).For this we used Top of atmosphere (TOA) reflectance data available in GEE (Gorelick et al., 2017).This includes data from geometrically corrected Landsat-4 and -5 Thematic Mapper and Multispectral Scanner system, Landsat-7 Enhanced Thematic Mapper and Landsat-8 Operational Land Images sensors.Any pre-processing steps are described in this section, followed the step-by-step decisions (steps 1-4 in Fig. 2) made for each pixel, to determine whether it was a candidate endmember for either of the representative landcover types.Such an Unsupervised Decision Tree (UDT) was used to consistently separate land, water and an initial estimate of the intertidal area.By combining this spatial context with spectral criteria, we selected image endmembers that represent characteristics found in the image at the time of acquisition (Shi and Wang, 2014).The abundance maps that result from the subsequent LSU can therefore indicate the mixed pixels, in our case related to diffuse boundaries.

Pre-processing
For each image the flags from the pixel assessment band were used to mask out clouds and shadows, as detected by the automatic cloud mask algorithm (Zhu and Woodcock, 2014).The here used blue, green, red, near infrared (NIR) and shortwave infrared (SWIR) bands have 30 m resolution.The thermal infrared (TIR) band with a resolution of 120 m was resampled to 30 m to match the other bands.From these bands the Normalized Difference Vegetation Index (NDVI) and Normalized Difference Water Index (NDWI) were calculated.The NDWI uses the green wavelengths to maximize the reflectance of water features and simultaneously takes advantage of their low reflectance in the NIR range (McFeeters, 1996): As a result of higher NIR reflectance, both vegetation and mud usually have low or negative values.

Land or water
After the pre-processing, land and water pixels were separated based on grey level histograms of the NDWI.Therefore a Canny Edge detection algorithm (step 1) was used on these NDWI values to distinguish locally between homogeneous land and water pixels (Liu and Jezek, 2004).This algorithm detects edges by looking for maximum gradients, or sharpest changes of values, in the NDWI image.A Gaussian pre-filter (GPF) filtered noise by smoothing the original image and a minimum gradient (MG) excluded edges with weaker gradients.Subsequently, edges were filtered with a minimal length of connected pixels to avoid selecting discontinuities between other land covers (e.g.urbanagriculture or cloudshadows).A buffer was applied around these edges, resulting in a set of spatial neighbourhoods.The pixels from these neighbourhoods form a bimodal histogram (Donchyts et al., 2016), emphasizing the difference in spectral properties in the Green and NIR reflection bands for water (NDWI values → 1) and land (NDWI values → -1) pixels.The adaptive Otsu thresholding algorithm (step 2) was then employed on the smoothed bimodal histogram, to separate the two dominant lobes with distinct mean values, corresponding to land and water values (Lu et al., 2011).This threshold was used to simultaneously mask terrestrial vegetation and bare ground.
Finally, the mask was updated with land cover classes urban, agriculture and forests from the yearly Modis landcover classification (Friedl and Sulla-Menashe, 2015).Derived from Terra and Aqua reflectance data, this layer (international Geosphere-Biosphere program classification) was used to ensure urban, agriculture and forest classes were included in the land mask.The relevant pixels were eroded by 120 m to avoid interference in the coastal waters and intertidal surfaces resulting from differences in pixel size.

Initial estimate intertidal zone
The spectral similarity between intertidal flats and turbid waters prevents the use of generic thresholds or aggregation of images when trying to separate them (Ryu et al., 2002).Thermal Infrared (TIR) bands are sufficiently sensitive to separate these two classes, but they lack the spatial resolution to extract the waterline (Sørensen et al., 2006).Therefore, an initial estimate of the intertidal zone was made using a thermal threshold (step 3) derived from the modal pixel value in the detected Canny Edge neighbourhood zone (from step 1).This temperature threshold was applied on all water pixels that remained after applying the land mask (step 2).In this way all exposed water pixels with a different radiant temperature were isolated, serving as the initial estimate of the intertidal zone.

Linear spectral unmixing
Unique image histograms, made up of from the canny edge neighbourhood zones and water mask, were used select end-member candidates (step 4).Thus, spatial information to extract the endmembers was added (Shi and Wang, 2014).For vegetation, NDVI values of all land pixels inside the detected canny edge neighbourhood zones were selected (see Section 2.2.1).This resulted in a histogram with a clearly distinguishable peak that corresponded to dense vegetation.Afterwards, pixels were selected based on a 5% buffer on both sides of this NDVI peak.For the water end-member candidates a similar approach was adopted, deriving the peak NDWI value from the water pixels (see Section 2.2.2) with a buffer of 1% around it.The driest mud pixels were selected from the initial estimate of intertidal zone (see Section 2.2.3), excluding vegetation by adopting a NDVI threshold < 0.3.From these sets of pixels, the mean values per band were derived to define the spectral signatures per end-member group.
Fractions of the reflectance signal were untangled by applying the spectral signatures for each endmember to a fully constrained standard linear mixture model (Alcântara et al., 2009): where n is the number of bands, p ij the surface reflectance of land cover type j in band i, f being the fraction of the pixel covered by type j and e i the error.This model assumes that the total reflectance is the sum of the components that make up the pixel, without interaction between them (Somers et al., 2011).This results in a set of fractions that linearly represent the proportion of each active component to the signal per pixel.Our LSU model was constrained, so that the fractional cover f j is always between 0 and 1, and such that the sum cannot be >1 for each pixel (R i ).

Decision tree performance
Because of limited availability of field observations that can relate end-members and LSU fractions to sediment concentrations and the subtidal footprint of mudbanks in Suriname, validation is restricted here to a robustness check of the described decision tree.Especially in muddy coastal environments the combination of canny edge detection and Otsu thresholding methods can be sensitive to the input parameters (GPF and MG), the buffer width and minimum edge length (Section 2.2.2.) (Bishop-Taylor et al., 2019).This sensitivity expresses itself in spatial variability of the spatial neighbourhood zones from which to sample end-member candidates, resulting from variability in the Otsu threshold, vegetation peak, water peak and temperature threshold.The temporal consistency of end-member signatures was therefore quantified by comparing changes in Otsu thresholds, vegetation peak, water peak and temperature threshold for parameter combinations for a subset of images, acquired between 2008 and 2010 over Paramaribo, Suriname (Fig. 1).Values of 0.3, 0.5, 0.7, 0.9 and 1 were used for the GPF; for the MG values of 0.7, 0.9, 1.1 and 1.3; minimum length values of 10, 25, and 75 pixels; for the buffer values of 5, 10 and 15 pixels were tested.

Image histograms
To show the benefits of using an UDT and the application of the resulting fractions for the purpose of analysing cross-shore gradients of mud, water and vegetation fractions, two relatively cloud-free Landsat-5 images were selected from the subset (see Section 2.4).The image acquired on 12 September 2009 was captured near high-tide while the second image, acquired on 15 November 2009 was captured near lowtide.A GPF of 0.7, MG of 0.9, buffer width of 10 pixels and minimum length of 25 pixels were used to derive the canny edge neighbourhood zones.Based on their different NDWI histograms, made up of pixels in these zones, image specific thresholds of − 0.186 and − 0.246 were derived to separate land and water pixels (Fig. 3, panel A).In the high- tide and low-tide images the NDWI peaks, between 0.422 and 0.438 and 0.524-0.534,correspond with water end-member candidates (Fig. 3, panel B).For vegetation end-member candidates the peak values, derived from histograms with pixels from the canny edge neighbourhood zones, were between 0.684 and 0.696 and 0.678-0.700(Fig. 3, panel C).Mud pixels were separated from water pixels that remained after applying the land mask, based on a radiant temperature threshold of 296 and 294.2 Kelvin for low-and high-tide (Fig. 3, panel D).These observations of unique threshold values illustrate the rationale behind the decision to apply an UDT procedure that automatically separates land, water and an initial estimate of the intertidal zone.Threshold values for the low-and high-tide image, derived with alternative input parameters are shown in Table A1 (Appendix A).

Fractions
In Figs. 4 and 5 fractions of vegetation, water and mud, from the selected high-and low-tide image respectively, are compared along a 30 km cross-shore transect (see Fig. 1).These profiles reveal that the outlined approach was able to generate abundance maps of the endmember classes that match expected cross-shore patterns.Namely, the mud fractions show a discontinuous decrease for the visually estimated subtidal part, from 0.95 to 0.50 during low tide (Fig. 4).The water fractions show a contrary pattern, implying variable SPM concentrations.At around 6000-6041 m the ratio of mud and water fraction changes more quickly; more specifically this rearmost rapid decline of mud abundance reveals a diffuse, and thus seaward mud-bank boundary.
During high tide (Fig. 5) absolute mud-fractions were significantly lower over the visually estimated subtidal mudbank, ranging from 0.30 to 0.10.This difference can potentially be attributed to the amount of SPM, variable tidal elevation, wave climate (Zorrilla et al., 2018) or difference in spectral signatures.For this reason, it is more difficult to use an absolute fraction value as indication of diffuse mud-bank boundaries.Still, like the low-tide transect, the rearmost decline in mud abundance around 5982-6000 m suggests the same seaward mudbank boundary.
The edge of the land mask remained fixed at its position around m, between 12 September 2009 and 15 October 2009.This location coincided with a decreasing vegetation fraction, from 1 to ±0.25 for both images, indicating that the detected boundary between land and water follows the mangrove fringe.The lower vegetation fractions seaward of this fringe may correspond to a change in vegetation type, the presence of microphytobenthos or even sparse mangroves standing in turbid waters.
The intertidal extent visible during low tide, that is, the cross-shore distance between the land mask edge and sea-mudflat boundary, coincides with scattered mud fractions between 817 m and roughly m (Fig. 4).For the high-tide image (Fig. 5), the intertidal zone was also visible from the initial estimate based on the TIR threshold.These results highlight that the image's TIR threshold did not match the local thermal difference between water and land locally, due to weather conditions

Robustness of the approach
To test the robustness of the UDT approach, all available Landsat images (88) between 2008 and 2010 were selected.A total of 68 images resulted in explicit end-member candidates that could be used to extract spectral signatures for LSU.The excluded images did not contain enough detected canny edge neighbourhood zones to sample endmembers from, while matching the UDT informed decisions.Fig. 6 shows the reflectance variation of all endmembers of vegetation, water and mud.Although images originated from different sensors, they show similar values and recognizable signatures within the same order of magnitude.
The NDWI threshold to separate land and water proved to be an important parameter in the UDT workflow.With each parameter combination in the canny edge computation the median value, range and distribution of threshold values changed slightly (Fig. 7), especially when comparing them to the stable NDVI peaks and TIR thresholds for the same parameter combinations (Figs.A1 and A2, Appendix A).For all parameter combinations the NDWI threshold ranged from − 0.326 to − 0.050, with outliers up to 0.034 (Fig. 8).Especially with MG values of 1.1 or 1.3, a higher GPF reduced the Otsu threshold, while the opposite can be seen for an MG of 0.7 (Fig. 7).This implies that the combination of MG and GPF influenced the edge detection, and thus the NDWI histogram used to separate land and water pixels.The NDWI threshold reduced also with an increasing buffer size (Fig. 8), because of the larger amount of land and mixed pixels included in larger buffers.The opposite can be seen when increasing the minimum length of detected canny edges (Fig. A3, Appendix A).This effect disappears for MG > 1.1, suggesting that the effect of filtering on a minimal edge length is comparable with applying a higher MG threshold.Hence, the higher the MG, the more likely these larger magnitude changes in NDWI are detected; with the potential of not detecting any boundaries at all with a MG > 1.1 and GPF > 0.5.These results show that it is possible to consistently select boundaries from NDWI images that correspond with transitions from land to water.This robustness is supported by Table A1 (see Appendix A) where 4 parameter combinations were applied on the selected low-and high-tide images.The variation in NDWI threshold resulted in differences in land-water boundaries, and thus selection of end-member candidates.Yet, both the signal to noise ratio of the end members and mean LSU error (see Eq. ( 2)) remained consistent.

Gradients and diffuse boundaries
By using spatial selection criteria for the definition of robust image endmembers, we developed an UDT for analysing diffuse boundaries of subtidal mudbanks.The pixel fractions resulting from the LSU represent mud, vegetation and water gradients, as shown in the selected transect for a low-tide and high-tide image (Figs. 4 and 5).Thus, we account for mixed pixels and avoid aggregating multitemporal image observations with different tidal heights.
The mud fractions show a discontinuous decrease in the offshore direction for both selected images (Figs. 4 and 5).Yet, more abrupt changes at 6041 and 5982 m in both mud and water fractions indicate changes in the presence of mud in the upper water column.These are considered fuzzy transitions between higher and lower mud Fig. 6.Variation in spectral signatures for the baseline scenario with a buffer of 10 pixels, a minimal length of 25 pixels, a Gaussian pre filter sigma 0.7 and a minimal gradient of 0.9.End-member signatures are taken from selected images between 2008 and 2010.concentrations, corresponding with diffuse boundaries of mudbanks.This indicates that, despite differences in tidal elevation and the time between acquiring the images, this seaward mud bank boundary remains relatively stable at approximately 6000 m offshore.The remaining difference can be attributed to mud bank migration and the temporal variability of SPM concentration, resulting from waves, tides and currents (Zorrilla et al., 2018).Absolute values of mud fractions in the estimated subtidal area also varied between the two images, with lower mud fractions for the high-tide image than for the low-tide image.This reflects the difference in sediment concentrations, preventing the use of an absolute mud fraction threshold to indicate diffuse boundaries.

Intermediate results
Besides the pixel fractions, the intermediate results, so far mainly used in the UDT procedure, show promising signs of added value.For example, Fig. 5 indicates that the intertidal zone can be estimated with the TIR band, when image quality is sufficient.In our case this allowed us to isolate pure mud pixels as endmember for the entire image.The undetected intertidal area in Fig. 4 is an example where the image TIR threshold locally didn't result in an estimate of the intertidal area.However, when the threshold is determined locally from well-defined spatial neighbourhoods, the intertidal zone can be demarcated more precisely.This is in line with findings from earlier studies that relate exposure time and intertidal topography to de difference in thermal irradiance between land and water (Ryu et al., 2002;Sørensen et al., 2006).
The boundaries of the terrestrial land masks, for example between 12 and 09-2009 (Fig. 5) and 15-11-2009 (Fig. 4), were consistently estimated at 817 m.As a result it becomes possible to analyse transitions from intertidal-to terrestrial landcover classes without using a static global or regional land mask product (Laengner et al., 2019).Especially in Suriname, where analysis of coastal morphology is hampered by limited data availability and changes are rapid and omnifarious, this type of spatial information can be beneficial to coastal managers.It allows for example to detect changes in vegetation composition for terrestrial and intertidal zones and incorporate that in coastal conservation-and protection measures.

Robustness & sensitivity
The sensitivity analysis on the NDWI threshold provided a better understanding of the capabilities of combining canny edge detection algorithm and Otsu thresholding for separating land and water pixels.The range of unique NDWI thresholds supports earlier observations (Bishop-Taylor et al., 2019) that especially in mud-dominated systems, a dynamic threshold is required to separate sea from land.
We found that the canny edge algorithm complemented with buffer Fig. 7. Variation in the NDWI threshold, used to separate land and water, for different input parameters sets variation of the minimal gradient (0.7, 0.9, 1.1 and 1.3) and Gaussian pre filter sigma (0.3, 0.5, 0.7, 0.9 and 1).Minimum length of the edges and buffer size around them are fixed at the defaults of 25 and 10 pixels respectively.
and length filters results in a consistent, spatially explicit and robust neighbourhood zone when the MG and GPF are set appropriately.By applying a sufficiently large length filter, only true land-water boundaries are included (Fig. A3).For a larger buffer, when more pixels are included, the threshold values shift towards the land lobe of the histogram.This implies that more pixels are included in the terrestrial land mask, resulting in the boundary shifting seaward.Nonetheless, the resulting neighbourhood zone consistently creates histograms that can be used to estimate representative image endmembers from ranges in both NDWI and NDVI values.This is supported by the relatively small temporal variation in spectral endmembers, sampled from these neighbourhoods, which we observed for all images between 2008 and 2010 in the area (Fig. 6).This facilitates the comparison of the linear spectral unmixing outputs for different images.

Limitations & improvements
A disadvantage of the outlined approach is that resulting fractions are only approximations of SPM or vegetation cover, as non-linear responses in reflectance spectra are not accounted for (Somers et al., 2011).Moreover, the spectral signatures of the image endmembers were not selected based on their pixel purity but rather on spatial and spectral selection criteria.However, monitoring requires an objective and standardized end-member extraction technique that allows for spatiotemporal analysis of LSU outputs (Figs. 4 and 5).The here discussed UDT approach consistently selects these endmembers and thus results in comparable LSU outputs between different dates and environmental conditions during acquisition.
Also, automatically estimating the position of the seaward boundary of mud banks from the LSU outputs remains a future improvement.As we showed in the fraction profiles (Figs. 4 and 5), the cross-shore transects reveal multiple rapid declines of mud fractions, indicating diffuse boundaries.Especially automatically selecting the rearmost rapid decline requires additional advancements in processing the LSU outputs and sufficient field observations for validation.This diffuse boundary position estimate is required to facilitate multitemporal position analysis of mud-bank footprints.The thresholds and index ranges used to define end-member candidates can then be used to assess the suitability of the image in a timeseries analysis, compared to for example only using estimates of cloud cover.

Conclusions
We developed a data-driven method that facilitates the detection of diffuse boundaries of subtidal banks along mud-dominated coastlines, as exemplified for a section of the Suriname coast.The developed unsupervised decision tree includes (1) advances in separating land and water with Otsu thresholding, (2) a novel approach to automatically estimate the intertidal zone and (3) an assessment of cross-shore gradients in mud, water and vegetation fractions from individual Landsat observations.We reaffirm that efficiently separating land and water is possible in mud-dominated coastal systems by defining a spatial neighbourhood with an edge detection algorithm applied to the NDWI index.However, selecting the appropriate input parameters is a non-Fig.8. Variation in the NDWI threshold (n = 68) used to separate land and water for different input parameter set: variation of the buffer size (5, 10 and 15 pixels), different minimum edge length values (10, 25, 50 and 75 pixels) and MG values (0.7, 0.9, 1.1 and 1.3) and the GPF is set at its default of 0.7.trivial exercise.The resulting terrestrial boundary allows the separation of the remaining exposed intertidal zone from water, based on a temperature threshold.Consistently sampling potential end-member candidates from these initial estimates of water, intertidal and land surface can be done from their index histograms.The resulting spatially explicit and spectrally coherent image endmembers facilitate multitemporal LSU outputs with fraction maps of water, vegetation and mud.From these, spatiotemporal differences in their sub-pixel proportions can indicate much needed information about changing gradients and the presence of coastal features.

Table A1
For the two example images five scenarios with different input parameters result in a different threshold, and thus a unique land mask, indicated by the absolute and relative amount of land pixels.Pure vegetation, pure water and pure mud indicate the number of pixels used to derive the mean end-member signatures.SNR is the sum of the standard deviation in reflectance for each band divided by the mean value, added up together for the 3 end member classes.The RMSE is the mean pixel error from the LSU (e i in Eq. ( 2)) analysis.

Fig. 1 .
Fig. 1.Median composite image of available Landsat images between 2008 and 2010 for the region of interest with in the bottom left panel, in white the intertidal area as estimated by Murray et al. (2019).The red line indicates the transect location used in this study, black pixels indicate masked clouds.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Workflow applied to individual image: (1) canny edge detection with the relevant parameters, (2) the Otsu thresholding on the NDWI histogram resulting from canny edge detection, (3) TIR thresholding on the TIR histogram to separate intertidal zone from turbid water, (4) defining the spectral ranges used to select end-member candidates.The resulting end-member graphs are used to apply linear spectral unmixing, resulting in maps representing fractions of water, vegetation and mud.

Fig. 3 .
Fig. 3.The histograms for the high-tide (left column) and low-tide (right column) image that are used to derive end-member candidates.Panel A shows the histogram (blue) used for separating land and water with the corresponding Otsu threshold for the selected images.In grey the histogram is shown corresponding to all pixels in the coastal zone.Panel B and C indicate the ranges of index values used for selecting the water and vegetation end-member (NDWI for water, NDVI for vegetation).The temperature threshold (in Kelvin), shown in panel D, indicates the value used to separate intertidal mud from water pixels and is based on the histogram derived from land pixels in the detected canny edge neighborhood zones.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Cross-shore patterns in pixel fractions for the defined end-member classes for the low-tide image (15-11-2009).The land mask, derived by applying a threshold value of − 0.246, is shown in panel A. Panel B shows the LSU fraction outputs in RGB (R = mud, G = vegetation and B = water).Panel C shows the crossshore development of fractions and the land mask boundary at 817 m.The subtidal extent of the mud-bank is visually estimated to align it with the rearmost rapid decline of mud abundance.Patches in panel A and B indicate masked clouds.

Fig. 5 .
Fig. 5.For the same location as in Fig. 5, now during high tide (12-9-2009).The in initial estimate of the intertidal extent and land mask, derived by applying a threshold value of − 0.186 are shown in panel A. Panel B shows the LSU fraction outputs in RGB (R = mud, G = vegetation and B = water).Panel C shows the crossshore development of fractions and the land mask boundary at 817 m.The subtidal extent of the mud-bank is visually estimated to align it with the rearmost rapid decline of mud abundance.Patches in panel A and B indicate masked clouds.

Fig. A2 .
Fig. A2.Temperature threshold variation (n = 68) for images acquired between 2008 and 2010 for different parameter sets.This threshold value was used to make an initial estimate of the intertidal area for the entire coastal area in the concerning Landsat image.Parameter sets include different buffer values (5, 10 and 15 pixels), minimum edge length thresholds (10, 25, 50 and 75 pixels), MG values (0.7, 0.9, 1.1 and 1.3) and the GPF at its default of 0.7.

Fig. A3 .
Fig. A3.Variation in the NDWI threshold (n = 68) used to separate land and water for different parameter sets: variation of the buffer size (5, 10 and 15 pixels), different minimum edge length values (10, 25, 50 and 75 pixels), MG values (0.7, 0.9, 1.1 and 1.3.)and the GPF at its default of 0.7.Difference with Fig. 8 is that the x-axis contains the buffer values to clearly show the effect of varying minimal length values on the thresholds.