Image Acquisition

Sentinel-2 imagery is freely available from the European Space Agency’s Copernicus Data Space Ecosystem (https://dataspace.copernicus.eu/). The most efficient strategy for searching the Sentinel-2 archive for the purposes of mapping eelgrass beds in the north Pacific is to determine which dates from June to September (ideally in July or August when eelgrass abundance peaks) had low-tides (less than 0.0, mean lower low water) during late-morning local time. We used the National Oceanic and Atmospheric Administration Tides and Currents website to acquire tide information for the station at Grant Point, Izembek Lagoon, Alaska (station ID 9463058, https://tidesandcurrents.noaa.gov/noaatidepredictions.html?id=9463058). After we created a list of dates with late-morning low tides, we used a satellite imagery viewing portal like National Aeronautics and Space Administration (NASA) Worldview (https://worldview.earthdata.nasa.gov/) to ascertain the cloud conditions on each date. We found most dates had cloudy conditions, but a few dates with complete or partial visibility of Izembek Lagoon were advanced to the next step of checking the Sentinel-2 archive to determine if (1) a Sentinel-2 image was collected on the respective date, and if so (2) does the image quality appear acceptable to meet the eelgrass-mapping objectives.

For mapping eelgrass, we specifically downloaded the Sentinel-2 Level-1C product (https://browser.dataspace.copernicus.eu/) that provided top-of-atmosphere reflectances in cartographic geometry. The Sentinel-2 imagery is distributed as Universal Transverse Mercator (UTM)/World Geodetic System 1984 (WGS84)-projected tiles of approximately 110 by 110-km dimension on a 100-km grid, so there is ample overlap with adjacent tiles if necessary. Sentinel-2 products consist of 13 bandwidths of spectral wavelengths: four bands at 10-m resolution, six bands at 20 m resolution, and three bands at 60 m resolution (fig. 1).

Our searches of the Sentinel-2 archive for the eelgrass growing seasons of 2008 to 2023 found two dates when clear sky, low tide, and image acquisition all coincided; the data from these dates were suitable for spectral classification and eelgrass mapping throughout Izembek Lagoon (table 1).

Image Pre-Processing

We downloaded the Sentinel Level-1C image tiles containing Izembek Lagoon for each of the two dates shown in table 1. The Level-1C product was delivered as terrain-corrected, top-of-atmosphere reflectances (the ratio between the incoming light and the light reflected after traveling through the atmosphere) in a georeferenced format (UTM Zone 3N), but otherwise, the Level-1C had the fewest alterations applied to the raw satellite image swath data compared to other Sentinel-2 image data formats. Each single tile from the two dates was cropped to the same 5,000-line by 6,000-sample image subset that fully encompassed Izembek Lagoon with 10-m pixel resolution and an upper left origin-pixel center-coordinate of ULX=6100005 and ULY=6155995. We resampled the 20-m bands to 10-m resolution using bilinear interpolation and stacked all bands into a 10-band composite image (table 2) for subsequent spectral analysis.

We originally attempted to identify areas of change by clustering a combined 2-date 20-band image into 80 classes using the unsupervised algorithm ISODATA in the open-source software package MultiSpec© (version 9.2011, https://engineering.purdue.edu/~biehl/MultiSpec/). The class statistics were used with a Gaussian maximum likelihood classifier to assign each image pixel membership into one of the cluster classes. Evaluation of that classification, however, showed that the classes most often associated with change were primarily affected by changes in the physical environment such as alterations in mudflats, high-tide areas, and deep-water channels as opposed to changes in the biological environment. Consequently, we only used the combined 2-date classification to isolate the lagoon area by grouping all classes within the lagoon (except for obvious deep-water classes) and masking all remaining classes (that is, upland vegetation and deepwater). An elevation mask (elevation greater than 0 m) was additionally applied to exclude onshore water. The resulting mask included the intertidal areas of Izembek Lagoon (not areas of deep water), and some neighboring embayments (that is, Cold Bay and Kinzarof Lagoon) that were within the footprint of the Sentinel-2 image subsection. This mask was applied to both 10-band sub-sectioned images (that is, 2016 and 2020), and thereafter the two masked images were analyzed independently.

Image Spectral Analysis

Image analysis involved two parts. First, image pixels were split into cluster classes, and second, cluster classes were grouped into spectral classes. The splitting phase used ISODATA and the maximum likelihood classifier to derive and map 40 cluster classes from each 10-band masked image; hence, the cluster classes were specific to Izembek Lagoon proper. Forty clusters were prescribed—roughly twice the desired number of final classes.

The second phase grouped the cluster classes into a set of classes that were homogenous both spectrally and spatially. A series of iterations were performed to group the cluster classes into spectral classes. The overarching strategy was to maximize the number of spectral classes while reducing duplication of very similar classes. We applied five approaches to aid in the grouping process:

These approaches were collectively used to guide a visual and numerical interpretation of the cluster classes with the primary goal of grouping similar cluster classes into spectral classes. First, cluster classes clearly associated with pure water or mud flats (but not eelgrass classes) were identified and merged. This merging resulted in 36 and 35 cluster classes for the 2016 and 2020 images, respectively. Next, classes with a small or very small number of pixels were grouped into the spectrally and spatially closest (that is, “nearest”) class. Then, cluster classes that were spectrally and spatially very similar were grouped. Finally, classes with no or very few ground data were grouped with their most similar class. These steps served to produce spectrally and spatially homogeneous spectral classes by aggregating similar entities from the clustering step and ensuring that each class had field data. Figure 2 illustrates the assignment of cluster classes to spectral classes, mostly by grouping but with a few instances of splitting and masking. The biplots in figure 2 show an example of spectral relations between just two of Sentinel’s 10 spectral bands, meaning figure 2 shows an incomplete representation of the many spectral relations among clusters. Figure 2 is presented here to show how biplots can be a useful tool for guiding cluster groupings.

Following critical inspection of the spectral classes, two types of targeted adjustments to “problem classes” were applied to attain better spectral continuity. One adjustment type involved splitting a spectral class and the other involved masking pixels for supervised class reassignment. Specifically, in the 2020 image, one cluster appeared to contain at least two different classes based on its spectral and spatial characteristics. That cluster was then split into 10 classes using the 10-band spectral data and the ISODATA clustering algorithm. The 10 classes were assessed spectrally and spatially and grouped into 2 classes resulting in one class primarily in the southernmost part of the lagoon and seemed to be associated with deep water, and the other class more in the northern lagoon associated with shallower water. The 2020 clustering also required a supervised reassignment to the classification. A few isolated cloud shadows occurred on the mud flats in the northeastern lagoon, so a mask boundary was digitized by hand to encompass each area of cloud and cloud shadow, and then the cluster classes within the cloud mask (in this case, two classes) were changed to match the class of the surrounding mud flats. A similar supervised fix was applied to the 2016 map that involved remapping affected classes to a surrounding mudflat class. The last adjustment to both maps targeted areas on mud flats where large drifting accumulations of dislodged eelgrass and green seaweed often get deposited by tidal actions, and not surprisingly, get included in spectral classes associated with eelgrass. Using expert knowledge, field experience, and prior eelgrass maps, we hand digitized a custom mask of areas where large entanglements of drifting eelgrass and green seaweed often wash up onto mud flats, and then recoded eelgrass-associated spectral classes within that mask to a mud flat class. After all the groupings and custom adjustments described above were applied, there were 29 spectral classes for the 2016 Sentinel-2 image and 24 classes for the 2020 image. The sequence of steps and outcomes described above are summarized in table 3.

Field Data Preparation

Field data (Ward, 2021) that were collected onsite at spatially coincident sample plots were used to attribute each respective spectral class with information about the ground cover characteristics (that is, surface types like eelgrass, mud, water, etc.). Many years of field data have been collected at Izembek Lagoon (Ward, 2021) to monitor changes in eelgrass abundance and distribution (2007–11, 2015–19, 2022–23). Recognizing that changes are intrinsic to eelgrass communities (Ward and others, 2003; Bartenfelder and others, 2022; Munsch and others 2023), using ground data collected within 1 or 2 years of a satellite image collection date reduces the chances of the field data misrepresenting the actual ground conditions as recorded by the satellite image. For the 2016 Sentinel-2 image, we used field data collected in 2015 and 2016, and for the 2020 Sentinel-2 image, we used field data collected in 2018, 2019, and 2022 (table 1). We also added ground data from plots sampled in 2007 and 2008 where no eelgrass was observed, but only if those plots were not resampled in subsequent years because they had since remained in high intertidal areas beyond the growth range of eelgrass in the lagoon (Ward and others, 1997). Including the 2007 and 2008 ground data collected at persistent mudflats bolstered sample sizes used to characterize associated spectral classes.

Data from each field sampling plot were recorded in a four 0.25-square meter quadrat formats, arbitrarily positioned diagonally 5–10 m from the plot location in each of four compass-oriented quadrants (northwest, southwest, northeast, southeast; Ward, 2021). We diagonally positioned the UTM location of each sampled quadrat by adding and (or) subtracting 10 m to the plot’s UTM easting and northing, respective to the quadrant’s diagonal (northwest, southwest, northeast, southeast) compass direction. Including each of the four sampled quadrants at each plot, as well as the 2007–08 unvegetated plots, bolstered the number of samples used to characterize the spectral classes with information about the associated ground conditions.

Four variables recorded at each plot quadrat were specifically relevant for quantifying eelgrass conditions associated with each spectral class (table 4; percent cover, Braun-Blanquet cover category [Braun-Blanquet, 1972], presence, and abundance index). If a plot was sampled in more than 1 year, averages were calculated for each quadrant subplot.

Spectral Class Attribution

The ground-plot quadrat data were spatially intersected with the spectral classes to calculate class-specific eelgrass cover statistics (mean, standard deviation, median, and interquartile range). Sometimes spectral classes were associated with notably different characterizations of ground cover. There are several explanations why a spectral class could be characterized by different ground data conditions, including (1) real changes in ground conditions between the dates of field sampling and satellite image collection; (2) accuracy of the field plot locations; (3) co-registration accuracies between the field plot locations and the georeferenced satellite image; (4) representativeness of the 0.25-square-meter quadrat relative to the 10 m×10 m Sentinel-2 grid cell; and (5) the choice of which Sentinel-2 grid cells to assign to each of the four sampled quadrants at each plot. For the last explanation, we assigned the sampled quadrants to the four grid cells on the diagonals of the grid cell that intersected the plot’s center location as opposed to the four cells in the cardinal directions.

Recognizing the high likelihood that some misrepresentative ground data gets included in calculations of class-specific eelgrass cover statistics, before those statistics were finalized, we first excluded the more egregious outliers of this ground data for each variable (table 4) within each spectral class. Outliers were defined as any data records with values outside (greater than or less than) the 1.5×interquartile range. After the outliers were excluded, we recalculated and graphed the class-specific statistics for each eelgrass variable as violin plots (app. 1, figs. 1.1–1.4) and tabulated the corresponding means, medians, and variances (app. 1, tables 1.1–1.8). The modeled eelgrass biomass estimates and 95-percent credible intervals (Patil, 2024) were calculated for each field plot location (not quad) so sample size intersecting each spectral class was four times smaller than the other field data metrics (app. 1, fig. 1.5). Consequently, we derived biomass maps from median values only (app. 1, tables 1.9, 1.10) because the mean values were sometimes unrealistically skewed.

Eelgrass Mapping

Our eelgrass-mapping strategy at Izembek Lagoon began with unsupervised derivation of statistical cluster classes that we subsequently grouped into spectral classes and then attributed with environmental data from field studies. In other words, the spectral classes established a gradient (or scale) of environmental conditions, and the ground data served to calibrate that scale. The underlying assumption of this strategy is that surface areas with similar cover characteristics will also possess similar spectral characteristics, which in turn will be grouped (clustered) into like spectral classes and thus provide useful information about the spatial distributions of surface-cover types throughout the study area.

Eelgrass maps were produced by attributing each spectral class with an eelgrass statistic of choice derived from field observations (after outliers had been excluded). Maps were generated for each of the five eelgrass variables in table 4. For example, a map of median eelgrass percentage of cover was produced by deriving the median for each spectral class from the ground plot data that spatially intersected the respective classes (fig. 3A, 3B), and similarly, a map of median eelgrass biomass was produced (fig. 3C, 3D).

All maps of spectral classes, eelgrass attributes, and derivative maps of eelgrass metrics in geoTIFF and CSV format are available in Douglas and others (2024).

Earlier Landsat-Derived Eelgrass Maps

Two prior eelgrass maps of Izembek Lagoon were produced from Landsat satellite imagery (Ward and others, 1997; Hogrefe and others, 2014). Both Landsat maps used ISODATA methods initially, like what was done for the Sentinel imagery described above, but thereafter, the methods of eelgrass mapping differed and are described here for comparison.

The first eelgrass map of Izembek Lagoon (Ward and others, 1997) was produced from a July 28, 1978, Landsat-3 multispectral scanner satellite image consisting of 4 spectral bands with 50-m pixel resolution collected during a late morning (11:14 a.m. local Alaska time) tide at −0.3 ft. ISODATA was applied to partition the image data encompassing the lagoon into 33 cluster classes, and then each pixel in the image was assigned to a cluster class using a maximum likelihood estimation (MLE) classifier. Each cluster class was subsequently labelled as water, eelgrass, unvegetated, or upland based on personal knowledge of the area, a review of black and white aerial photographs taken in June 1976, and ground data collected in 1986. Two final steps involved changing pixels classified as eelgrass in upland areas to the upland class, and similarly, changing small areas of eelgrass (≤2 pixels) surrounded by unvegetated areas to unvegetated.

The second eelgrass map of Izembek Lagoon (Hogrefe and others, 2014) was produced from a combination of two Landsat images: one collected on August 2, 2002, by Landsat-7-ETM+ at a tide level of +0.03 feet and the other on July 20, 2006, by Landsat-5-TM at a −0.27 tide. Two images were used (instead of one) to overcome some isolated cloud contaminations that did not overlap spatially. Each image was classified independently by first masking dry land and clouds and then using ISODATA to derive 35 cluster classes. Clusters were assigned to one of eight substrate classes to create eight training regions: four for eelgrass and two each for bare ground and water (or left unassigned if uncertain). Clusters spatially coinciding with areas of assumed eelgrass based on false-color enhancements of Landsat bands 4, 5, and 1 were manually assigned to one of four eelgrass classes based on the relative strengths of band-4 (near infrared) and band-3 (red) radiances: (1) high NIR, (2) moderate NIR, (3) low NIR, and (4) minimal NIR. Similarly, classes representing bare ground or water were determined by visual interpretation and partitioned into two respective subclasses, also based on spectral radiance profiles: dry or wet bare ground and shallow/turbid or deep/clear water. Then, all assigned clusters were merged to create discrete training clusters for a MLE classification that assigned every pixel to one of the eight substrate classes. To overcome small amounts of cloud contamination, the two maps were mosaiced with the 2002 map covering the southern one-third of the lagoon and the 2006 map covering the northern two-thirds. A final manual step changed small, isolated patches of eelgrass on bare ground to a ninth map class, seaweed. We refer to this map hereinafter as the 2006 Landsat map.

The initial unsupervised derivation of roughly 35–40 cluster classes by use of ISODATA was common to the Landsat and Sentinel satellite mapping approaches used in this study. Thereafter, however, both the Landsat approaches used a greater degree of subjective supervision than did the Sentinel approach. For the 1978 Landsat-3 map, the initial suite of 33 cluster classes were each manually assigned to one of four map classes based on expert interpretations of aerial photographs, spatial context, and field data. For the 2006 Landsat-7 map, the initial 35 cluster classes were manually assigned to one of eight training classes based on spatial context as well as spectral relations among the red and infrared imagery bands. Clusters assigned to each of the eight training classes were then merged, and those eight merged clusters were used with MLE to produce a map specifically of those eight defined substrate classes.

The Sentinel method began by developing a set of cluster classes with ISODATA, and then some were split or merged based on spectral relations, spatial juxtapositions, and sample sizes. Few presumptions or supervised attributions about surface cover were imposed on cluster assignments. The set of clusters was used instead to produce a map of spectral classes, and then each spectral class was attributed with surface-cover information by quantifying ground data collected within each class’s spatial footprint. In other words, our Sentinel mapping method transformed a static map of derived cluster classes into a map of surface cover by attributing each class with a summary of ground data about a metric of interest (for example, fig. 3). Having an adequate representation of ground conditions to robustly characterize the spectral classes, this strategy provides objectivity and utility owing to the breadth of map themes that can be derived for a variety of mapping goals.

Eelgrass Change Detection

Eelgrass changes at Izembek Lagoon can be estimated and quantified by subtracting two eelgrass maps depicting conditions in different years. The two maps must represent the same eelgrass metric, so the areas where the two maps disagree should presumably show areas where eelgrass had increased or decreased. However, map-detected changes are susceptible to false positives caused by factors other than real on-the-ground changes. Two influential factors stem from differences in (1) the precise tide level at the time when the images were collected; and (2) the seasonal phenology of eelgrass growth on the date when images were collected. Areas of eelgrass growth that are partially submerged at low tide and (or) possess only sparse eelgrass cover will be more challenging to robustly detect in imagery that is collected at a time when the tide is only slightly negative. Similarly, analyzing imagery of eelgrass that is many weeks from the time of peak growth will tend to diminish eelgrass detectability. These two factors (tide level and phenology) should be considered when interpreting map-derived changes, especially changes involving areas associated with marginal (sparse) eelgrass cover or deeper water.

The 2006 Landsat map of eelgrass at Izembek Lagoon delineated four spectral classes associated with areas of eelgrass with varying levels of near-infrared reflectance (Hogrefe and others, 2014). For the purposes of detecting eelgrass changes between the 2006 Landsat map of eelgrass and either the 2016 or 2020 Sentinel-2 eelgrass maps, we considered all four Landsat eelgrass classes to be representative of eelgrass presence. For the two Sentinel-2 eelgrass maps, we considered all spectral classes with a median Braun-Blanquet cover value greater than 1.0 (app. 1, fig. 1.6) to be comparable to the four combined Landsat classes of eelgrass. With those definitions, eelgrass presence/absence maps were generated for each of the three mapping years (2006, 2016, and 2020). Each pair of these derivative two-class (present/absent) eelgrass maps was then differenced to spatially depict where eelgrass was gained, lost, or unchanged from the earlier map year to the later map year. For illustration, we show results for the three difference maps (fig. 4), as well as the net areal extent (in square kilometers [km²]) of those mapped changes in eelgrass presence (table 5). Roughly 40 km² of area that supported eelgrass in 2006 was lost by 2020, whereas about 15 km² of eelgrass was gained, resulting in a net loss of about 25 km² (table 5). From 2006 to 2020, most areas where the presence of eelgrass was lost were in the central part of Izembek Lagoon (fig. 4).

Spatial changes in estimated modeled biomass from the first Sentinel-2 imagery (July 1, 2016) to the second imagery e (August 14, 2020) are shown in figure 5 and were calculated by subtracting the two biomass maps from figure 3. Although eelgrass biomass changed substantially over the course of a growing season (Ward and others, 2022), the biomass model adjusted for day of year and tide level (Ward and Amundson, 2019), so the two maps were better standardized with respect to phenology and thus better suited for revealing between-year changes in community density or vigor as might be affected by storm events, water temperature fluctuations, or winter ice scouring.

We multiplied the biomass estimate for each spectral class in both Sentinel-2 maps (as mapped in fig. 3) by the respective class’s aerial extent, and then summed across all classes to derive an estimate of the total eelgrass biomass throughout Izembek Lagoon on each respective image date (table 6). The same methods used to attribute spectral classes by intersecting the plot estimates of biomass and assigning the spectral classes with median values after outliers were removed were similarly applied using each ground plot’s upper (97.5-percent) and lower (2.5-percent) bounds of the 95-percent biomass credible interval (table 6). The modest net increase in lagoon-wide biomass from 2016 to 2020 (table 6) was inconclusive considering the broad 95-percent credible intervals of the annual point estimates.