*2.2. Data Description and Correction*

The imagery used in this study is Sentinel-2 A and B data produced by the European Space Agency. The mission is composed of two multispectral satellites that collect images that cover a large spatial extent, have a fine spatial resolution (up to 10 m for half the electromagnetic bands they detect), and a five day revisit time using both sensors. Both satellites carry optical sensors that sample in 13 spectral bands at varying spatial resolutions (Table 1) [25,26]. The high temporal resolution of the two satellites allowed us to assemble a time series for quantifying marsh change despite the variable cloud cover and inundation of CSMR which would prevent accurate image analysis. Imagery dates were selected to represent similar times of year, tide, and cloud cover. Four dates were selected to establish pre-flow, immediate, and post-flow conditions (approximately one and three years after the initial Montecito Debris Flows). Imagery from 13 November 2017 was used for pre-flow conditions as it was the date closest to the debris flow in which the marsh was not flooded or covered by clouds. Imagery from 12 January 2018 represented the immediate conditions as the data were collected three days after the flow occurred and before mechanical clean up and king tides occurred. Lastly, 3 November 2018 and 12 November 2020 imagery represented two recovery steps and were chosen to be consistent with the pre-flow November image. No November 2019 imagery was selected, as all available images were collected when there was either dense cloud cover or the marsh was inundated by high tide. All Sentinel-2 imagery was downloaded from the USGS Earth Explorer portal [27].


**Table 1.** Sentinel-2 band descriptions.

Imagery was preprocessed in the Sentinel Application Platform (SNAP) prior to being implemented in ENVI Classic 5.5.3 [28,29]. First, the Sen2Cor SNAP add-on was used to perform atmospheric correction to obtain bottom of atmosphere L2A imagery from the top of atmosphere L1C imagery downloaded from the USGS [30]. This process produced 12 atmospherically-corrected L2A bands. Once corrected, bands 1, 9, and 10 were removed as they are primarily used for atmospheric properties and are too coarse (60 m resolution) to be used in assessment of the fine scale change in the marsh. The remaining 20 m resolution bands (bands 5, 6, 7, 8A, 11, and 12) were then resampled using pixel replication to match the 10 m resolution of bands 2, 3, 4, and 8. Resampled and native 10 m resolution bands were layer stacked for further processing in ENVI.

High density LiDAR was also used in addition to Sentinel-2 imagery to assess conditions immediately after the debris flows occurred (January 2018). The LiDAR data were collected over the areas affected by the debris flows soon after the event by the Federal Emergency Management Agency (FEMA) at a density of at least 4 points per square meter [31]. LiDAR data were corrected and processed using the BCAL add-on for ENVI [29,32]. Data were height filtered at a threshold of 30 m with a 10 m search window. Height filtered

data were then processed using last returns into a digital terrain model (DTM) with 10 m resolution to match Sentinel-2 data.

#### *2.3. Spectral Analysis*

Before classification, corrected Sentinel-2 images were processed to obtain fractional cover and to calculate the normalized difference vegetation index (NDVI) and modified anthocyanin reflectance index (mARI).

Fractional cover was obtained via MESMA using the following steps. First, two spectral libraries were generated using the November 2017 and January 2018 processed images. Endmembers were selected based on site knowledge and similarity of spectra to those that would be expected for each landcover class (November endmember spectra are shown in Figure 2). Both libraries had endmembers selected to represent four broad landcovers that are expected in the marsh: non-photosynthetic vegetation (NPV), green vegetation, bare soil, and subtidal (water). Libraries were optimized using the endmember average root mean square error (RMSE) (EAR) [33], minimum average spectral angle (MASA) [34] and count-based endmember selection (CoB) [35], or the EMC option in VIPER Tools, and included a minimum of four sample endmembers per landcover class [36]. The November library was used for the pre-debris and two recovery images. The January 2018 image had a separate spectral library for the unique conditions that were expected around the debris flow.

**Figure 2.** Spectral library for November 2017. Axes are (*x*-axis) wavelength in nm and (*y*-axis) reflectance values. Spectral signatures for (**a**) bare soil, (**b**) green vegetation, (**c**) non-photosynthetic vegetation (senesced), and (**d**) subtidal cover classes. Colors used for visual differentiation of the different endmembers.

With the libraries generated, MESMA was performed to obtain fractional cover for all four dates. Endmember models used were 2-, 3-, 4-, and 5-endmember models to ensure the inclusion of the model approaches recommended by Rosso et al. (2005) [12]. All models were constrained to fractional cover between 0.0 and 1.0, shade fraction between 0.0 and 0.8, and a maximum RMSE of 0.025. This process produced fractional cover for the four endmember classes for each date (Figure 3 for example).

**Figure 3.** Fractional cover for November 2020 shown as an example of MESMA results. Scale indicates proportion of pixel that is made up of respective endmembers or shade. Pixels with a 0% cover were removed for visual clarity.

To further build on the data that would guide the classification of CSMR, two vegetation indices were calculated from the Sentinel-2 imagery: NDVI (Equation (1)) [37] and mARI (Equation (2)) [38,39].

$$\text{NDVI} = \frac{\left(NIR - RED\right)}{\left(NIR + RED\right)} = \frac{\left(Band8 - Band4\right)}{\left(Band8 + Band4\right)}\tag{1}$$

NDVI is one of the vegetation indices recommended in the literature for wetland analysis and was found to be one of the more important factors in classifying landcover classes in CSMR in prior work [6,9,10,19].

$$\text{mARI} = 800 \text{ nm} \cdot \left(\frac{1}{550 \text{ nm}} - \frac{1}{700 \text{ nm}}\right) = Band8 \cdot \left(\frac{1}{Band3} - \frac{1}{Band5}\right) \tag{2}$$

mARI is used to detect the levels of anthocyanins, a family of red pigments that can be related to stress and senescence in plants [40]. Anthocyanin content in *Salicornia pacifica* has been found to increase in the fall and winter [41]. Therefore, mARI has the potential to further help the classification of both senesced vegetation and a dominant marsh plant in CSMR.

Once the Sentinel-2 and LiDAR products were produced, data were layer stacked prior to the creation of training data. Training data were produced for five landcover classes bare soil, high marsh, mid marsh, senesced, and subtidal—by selecting reference polygons that matched regions of corresponding landcover from an expert map and from a report of landcover prior to the debris flow developed by Myers et al. (2017) (see Table 2) [23]. The high marsh class represents a mixed plant community of *Salicornia pacifica*, *Arthrocnemum subterminale*, *Frankenia salina*, and *Distichlis spicata*. Mid marsh represents portions of the marsh dominated by *S. pacifica*. The senesced landcover is composed of upland regions dominated by non-native shrubs and grasses [23]. Training data were collected for each date by creation of rectangular polygons in ArcGIS (Table 2 and Supplementary Materials Table S1) [42]. Training data and layer stacked images were analyzed in R [43].


**Table 2.** Polygon counts and decision metrics for training data generation.

#### *2.4. Random Forest and Change Detection*

We used a random forest classifier to assign landcover class to each pixel. Random forest is a machine learning technique that automates the categorization of data by running a datapoint (e.g., a pixel) through a set number of decision trees and picking a finalized landcover class via majority vote. Pixel values were first extracted from the layer stacked images with the values and associated landcover recorded into a data frame, which was then filtered to remove variables with NA/NULL values. The data frame was then read into the random forest algorithm, with n = 500 decision trees. This process produced classified maps of the five landcover classes. Additional outputs include: (1) variable importance, a measure that identifies which layer stack inputs were important in the landcover classification; (2) mtry accuracy, an accuracy assessment metric of the final model; and (3) kappa values, a secondary accuracy assessment metric. Final model selection and accuracy assessment were done via k-fold cross validation. Post-classification change detection was performed in ENVI using the change detection statistics option. Dates were compared to each other in both chronological order (i.e., November 2017 to January 2018, November 2018 to November 2020) and net order (November 2017 to November 2020). Comparing the dates this way allowed us to track landcover and to obtain extent for all five classes as time progressed, thus obtaining net change for each landcover class in the system. ENVI reported change statistics in terms of pixel count, area in square meters, and percentage change. These statistics include class differences and image differences. Percentage change was recalculated using both pixel count and area and used in place of the ENVI reported percentages.

#### **3. Results**
