Remote Sensing for Wetland Mapping and Historical Change Detection at the Nisqually River Delta

: Coastal wetlands are important ecosystems for carbon storage and coastal resilience to climate change and sea-level rise. As such, changes in wetland habitat types can also impact ecosystem functions. Our goal was to quantify historical vegetation change within the Nisqually River watershed relevant to carbon storage, wildlife habitat, and wetland sustainability, and identify watershed-scale anthropogenic and hydrodynamic drivers of these changes. To achieve this, we produced time-series classiﬁcations of habitat, photosynthetic pathway functional types and species in the Nisqually River Delta for the years 1957, 1980, and 2015. Using an object-oriented approach, we performed a hierarchical classiﬁcation on historical and current imagery to identify change within the watershed and wetland ecosystems. We found a 188.4 ha (79%) increase in emergent marsh wetland within the Nisqually River Delta between 1957 and 2015 as a result of restoration efforts that occurred in several phases through 2009. Despite these wetland gains, a total of 83.1 ha (35%) of marsh was lost between 1957 and 2015, particularly in areas near the Nisqually River mouth due to erosion and shifting river channels, resulting in a net wetland gain of 105.4 ha (44%). We found the trajectory of wetland recovery coincided with previous studies, demonstrating the role of remote sensing for historical wetland change detection as well as future coastal wetland monitoring.


Introduction
Coastal wetlands are ecologically and economically beneficial ecosystems, however, they are at significant risk of degradation and loss. Coastal wetlands contribute to habitat biodiversity, coastal protection, fishery health, and carbon sequestration [1]. In the contiguous United States, wetlands cover 5% of total land area, of which only 5% are tidally influenced coastal wetlands [2]. These are some of the most readily declining ecosystems globally. Within the United States alone, the rate of coastal wetland decline between 2004 and 2009 was 324.3 km 2 annually, with the majority of salt marsh loss attributed to coastal storms and land development [3].
Recently coastal wetlands have gained attention as the value and potential of coastal ecosystems for carbon sequestration and storage is being realized. Coastal wetland blue carbon is the carbon found in coastal vegetated ecosystems, such as mangroves, tidal marshes, and sea grasses, where it is stored in soils, above-and belowground biomass, litter, and dead wood [4]. Carbon sequestration rates for tidal saline wetlands and mangroves are globally estimated to be at 44. 6 TgC/year, and the release of greenhouse gases from these wetlands are minimal [5]. As a result, there is a growing interest in including tidal salt marsh in climate mitigation strategies and policies [6].
technologically appropriate spectral bands given the time period collected, and provided sufficient spaces in time to identify long-term changes. Secondly, we performed a change analysis over the 60 year period to identify changes in specific habitat types that can be distinguished from stable isotope analysis, including C3 salt marsh, C4 salt marsh, freshwater marsh, transitional marsh, and woody riparian habitats [37], to support work on the flow of carbon in Nisqually wildlife food webs. Finally, we aimed to identify and quantify trajectories of the emergent marsh habitat and forests to help understand the effects of restoration and other watershed-scale processes.

Study Site
The Nisqually River Delta is located 16 kilometers east of Olympia, WA, in southern Puget Sound ( Figure 1). The Billy Frank Jr. Nisqually National Wildlife Refuge extends into both Thurston and Pierce County where the Nisqually River reaches Puget Sound to form an estuary that provides a nutrient rich habitat for fish and wildfowl. From its headwaters at Mount Rainier, the Nisqually River carries 1,200,000 metric tons of sediment annually, most which remains contained behind the La Grande and Alder dams, built in 1912 and 1945, respectively [38]. As a result, only 100,000 metric tons of sediment a year is transported in Puget Sound. The Nisqually River and McAllister Creek bring fresh water to the delta with average flow on the Nisqually River ranging from 1000 to 2000 cubic feet per second (cfs) during low seasonal flows, and flood rates reaching 13,000 cfs [38]. High water years and increased discharge during flooding account for the majority of river channel erosion and increased sediment delivery along the Nisqually River to the Delta [34,38].  Habitats within the estuary are defined by the U.S. Fish & Wildlife Service as estuarine, riparian, and upland habitats, based on the degree of tidal and freshwater influences, as well as elevation. Unconsolidated shore, aquatic bed, and estuarine emergent marsh are directly affected by daily tidal inundation, while freshwater and riparian wetlands act as transitional zones where the influences of freshwater increases moving closer to upstream riparian forest. Tide levels in the Nisqually estuary range greatly, with a mean higher high water line (MHHW) of 13.5 ft. and mean lower low water line (MLLW) of −3.5 ft. [34]. The impacts of this extreme tidal range extend 5.3 km upstream along the transitional marsh into riverine forests [39]. Plant associations in the estuarine emergent marsh are spatially organized by tide inundation exposure and proximity to the fresh water supplied by the Nisqually River [36]. In the Nisqually River Delta, salinity correlates to specific species, while elevation is related to species richness. Within the salt marsh Carex lyngbyei (Lyngbye's sedge), Juncus baliticus (Baltic rush), and Potentilla anserine (silverweed) are commonly found in brackish areas, while Distichlis spicata (salt grass), Salicornia pacifica (pickleweed), and Jaumea carnosa (marsh jaumea) occupy areas of saline marsh [40]. Upland habitats occur at higher elevations outside of tidal and riparian influence, and are made up of mixed deciduous and evergreen forests, grassland, and agriculture [34].
The Nisqually Indian Tribe has derived cultural, spiritual, subsistence, and economic value from the Nisqually River Delta and watershed for millennia. Over the last century the Tribe has suffered dramatic losses to their salmon fisheries as salmon populations throughout the region have declined. In 1999 the Nisqually Fall Chinook (Oncorhynchus tshawytscha) stock became one of the 27 stocks in the Puget Sound evolutionarily significant unit listed as threatened under the federal Endangered Species Act [41]. Recovery planning identified the continued preservation of the Nisqually River Delta coupled with extensive restoration of altered estuarine habitat as the highest priority habitat action for the recovery of naturally spawning, self-sustaining Nisqually Fall Chinook [41]. Chinook salmon rear extensively in estuaries and are thought to be the most estuary-dependent of the Pacific salmonids [42]. The estuary is also important to other salmon, which are significant for tribal fisheries, like the Nisqually winter chum (Oncorhynchus keta), which has one of the largest wild runs in Washington State [43] and is known to utilize the estuary for feeding and growth [44,45].
Extensive historical anthropogenic changes to the Nisqually River Delta and watershed began in the late 1800s as pioneers altered the landscape for agriculture and development, including draining and diking a large portion of the salt marsh for growing crops, raising animals, and dairy farming. Logging began around the same time, continuing through the early 1990s until efforts shifted towards reforestation. In 1974, the Billy Frank Jr. Nisqually National Wildlife Refuge was established, allowing the diked agricultural areas to transition to rain-fed freshwater wetland [34]. The spread of invasive weeds, such as reed canary grass (Phalaris arundinaceae), degraded the freshwater wetland and the external dike was costly to maintain. Through a series of public review, the Refuge completed their Comprehensive Conservation Plan that called for the removal of the old Brown Farm Dike and enhancements for managed freshwater units within the Refuge. Initial restoration by the Tribe began in the Pilot management site in 1996 ( Figure 2), with the removal of dikes from tribal property, followed by the Phase I site in 2002 and the Phase II site in 2006 [40,46]. Restoration efforts culminated in 2009, with the removal of the Brown Farm Dike, with a total of 308.37 ha restored to estuarine conditions and 9800 m of channels to tidal flow [47,48].

Imagery and Data
For mapping vegetation and habitats within the study site, we acquired imagery from July 1957, July 1980, and July/August 2015. The 1957 imagery was comprised of nine overlapping highresolution panchromatic images scanned at 1000 dots per inch from the Aerial Photography Single Frame Records collection from the EarthExplorer website (US Geological Survey, 1957, 1:32,000

Imagery and Data
For mapping vegetation and habitats within the study site, we acquired imagery from July 1957, July 1980, and July/August 2015. The 1957 imagery was comprised of nine overlapping high-resolution panchromatic images scanned at 1000 dots per inch from the Aerial Photography Single Frame Records collection from the EarthExplorer website (US Geological Survey, 1957, 1:32,000 (JALBTCX), NOAA Office for Coastal Management, 2014), and 2015 1 m National Agriculture Imagery (NAIP). Additional high-resolution imagery and sampling datasets were used for field confirmation.
Given the large tidal range at the Nisqually River Delta, tide stage during image collection was an important consideration for wetland classification. Dupont Wharf in the Nisqually Reach has a subordinate tide station based on the Seattle, WA reference tide station (ID 9447130) using a 49-min time offset and height offset of 1.04 ft. for low tides. The 2015 WV-3 images were collected on 10 August 2015 at 11:22 Local Standard Time (LST), with an estimated tide during acquisition of 1.8 ft. (National Oceanic and Atmospheric Administration, 2017). The 2014 color-infrared image was acquired on 14 July 2014, specifically to capture the low tide extent of the delta flat, which we calculated to be near −2 ft. Tide levels were unavailable for earlier image dates.

Pre-Processing
Data pre-processing included image georeferencing, atmospheric correction, and image mosaicking. The 1957 panchromatic, 1980 CIR, and 1980 true color images were georeferenced to a 2013 NAIP image using ground control points and the polynomial geometric model in ERDAS Imagine 2014 [49]. A minimum of 20 ground control points were identified for historical and reference image and rectification error for all images was acceptable with RMSE below 0.5. For the 2015 WV-3 image, we performed an atmospheric correction using ATCOR software [49] and georeferenced the images in ERDAS Imagine using the rational polynomial coefficients (RPC) ortho-rectification tool and ground control points. All images were projected to UTM NAD 83 zone 10 and the 1957 and 1980 images were resampled in ArcGIS using the 'nearest' technique to a spatial resolution of 1.28 m for alignment during change detection. Images for each year were mosaicked and clipped to the proposed and approved Billy Frank Jr. Nisqually National Wildlife Refuge boundaries [34] with an extension to the southwest to include the Nisqually Clear Creek Hatchery.

Ground and Ancillary Data for Training and Validation
Reference data for the habitat and species classification varied in quality and availability depending on classification year. The 2015 classification benefited from the availability of repeated field observations and high spatial resolution data sources for ground-truthing. Habitat classes were confirmed with visual interpretation of 2014 1 m CIR imagery and 2015 1 m NAIP imagery. Field collected species data were compiled from annual USGS 0.25 m 2 quadrat, circle plot, and transect vegetation sampling data, 2015 biomass and species study data [50], field observations from March 2016, and vegetation samples collected in 2012 by Thorne et al. [51]. Ideally, species samples would be collected around the same date as the 2015 image acquisition; however, the addition of 2012 and 2016 data improved the spatial coverage of samples across the study area. Species samples were compiled into a single, standardized database and data points were removed if soil cover was greater than 25 percent, or if the dominant species had less than 70 percent coverage. All samples were selected using a stratified random sampling method.
Historical image classifications are particularly challenging for assessing accuracy due to limited or non-existent historical reference data, in turn necessitating the reduction and simplification of classes [52] and the inclusion of several data sources when available [53]. Due to limited availability of historical data, 1957 habitat samples were ground-truthed with visual interpretation of 1961 panchromatic air photos from the Washington State Geospatial Data Archive in conjunction with the 1957 aerial photographs used for classification, while 1980 habitat samples were validated with the Thurston Met true color air photos from 1980. For all reference samples based on high resolution imagery, those that could not be confirmed through visual assessment were discarded. Early wetland vegetation data was unavailable; therefore, we used a Nisqually salt marsh vegetation distribution map of 12 plant associations created by Burg et al. [36] as species reference data for 1957 and 1980. The Nisqually salt marsh vegetation map [36] was created in July 1975, based on data collected from 138 0.25 m 2 quadrats 30 m apart, along a total of 12 transects. The plant associations were identified using an indirect polar ordination method, while the vegetation map was prepared using field sampling and a color aerial photograph. Accuracy for the vegetation map was not available. We georeferenced and digitized the vegetation map, using levee boundaries and marsh contours as ground control points, which were limited due to the generalized nature of the map. We reclassified plant associations on the salt marsh map into dominant species groups and randomly selected points to be used as training and validation samples.

Classification Scheme
Our classification schemes focused on habitats, photosynthetic pathway functional types, and emergent marsh species. The habitat classifications identified the presence of key land cover types and marsh habitats that included agriculture/grassland, developed land, emergent marsh, forest, riverbank, mudflat, and water. For the 1980 and 2015 classifications, the availability of additional spatial data allowed for further defined classes for wetland type (emergent, transitional) and forest type (riparian, upland, transitional), while in some cases, additional classes indicated newly present classes within the images, including outer mudflat, freshwater/managed wetland, and dead trees. The pathway classification focused on two species classes, with salt grass designated as a C4 plant, and all other emergent marsh plants labelled as C3 plants. The species classification identified dominant plant species located within the emergent marsh.
In the existing historical vegetation map, many of the 12 plant associations contained overlapping species and were spectrally inseparable, therefore we condensed plant association classes into six species classes based on their dominant species (Table 1) and spectral similarity as defined by the eCognitition Developer 9 [54] sample editor tool. As a result, the Distichlis Spicata-Salicornia pacifica association was included in the salt grass class. Still, the high variability of species distribution within the classes and limited spectral information provided by the 1957 aerial image resulted in poor spectral distinction between the six species groups during testing and as a result, emergent marsh species for 1957 were not mapped. The six species group classes were spectrally separable in the 1980 images and these classes were further combined into C3 and C4 plants classes.
Classes for the 2015 species classification were based on dominant emergent marsh species rather than plant association, based on the sampling methods used to collect the 2015 field data. We selected a limited number of key classes based on species of interest and marsh zone indicator Sustainability 2017, 9, 1919 9 of 32 species, including pickleweed, reed canary grass, algae, Lyngbeyi sedge, Typha latifolia (broad-leaf cattail), salt grass, as well as a combined single class for all other emergent plants.

Image Classification
We performed a hierarchical classification on each image using an object-based approach with eCognition 9 software [54] to identify and group pixels with similar feature values into segments that could then be classified. The object-based segmentation approach works well with hierarchical analysis in wetland environments, as it is suited for multi-level processing where discrete land cover types can be identified using expert knowledge [29]. Rule-based hierarchical classifications require initial set-up time, but the benefit of using ancillary data and expert knowledge can improve results over other statistical models [27].
For all classifications, our rule-set approach utilized hierarchical processing that alternated between image segmentation and threshold-based classification based on methods developed by Lewinski & Bochenek [55]. Image segmentation was performed using a multi-resolution segmentation tool that uses weighted spectral, brightness, and textural features to create image objects. Scale parameters were determined through trial and error to find appropriate object sizes based on hierarchical level and class type. A subset of training samples was used to select the reflectance threshold parameter we used to discriminate class membership. Classification occurred on a class by class basis beginning with the identification of land and water, which was then followed by further segmentation and classification of sub-classes. For the species classification we focused on the emergent marsh habitat to identify wetland vegetation types using the robust machine learning algorithm linear support vector machine classifier (SVM).

Classification Methods 1957 and 1980
The classification of historical air photos can be challenging due to the lack of historical ancillary data available and the limited spectral information of a panchromatic image; therefore we created image derivatives from our 1957 and 1980 images to improve classifier performance. A Euclidean distance texture image was produced in ERDAS Imagine from the panchromatic 1957 image to aid in distinguishing features with similar spectral values in the panchromatic band. To remove image noise and better identify variance within the images we performed principal component analysis (PCA) on the 1980 image. For the final 1980 classification we include the CIR, true color bands, brightness, and PCA bands. In some cases image objects could not be defined within our sample thresholds so we utilized visual interpretation and manual classification. Due to the lack of historical elevation data, upland and riparian forests were not classified in the final 1957 and 1980 habitat maps.

Classification Methods 2015
Availability of spectral and spatial data increased for the 2015 classification, resulting in more features used in the segmentation and classification based on a range of spectral bands, index ratio, and derivative features, including PCA. In addition to the eight WorldView-3 bands, a variety of normalized difference indices were produced from band combinations ( Table 2) [56]. Also for the 2015 classification we produced a 1m spatial resolution canopy height model (CHM) to further differentiate salt marsh vegetation from forest. The CHM was created by interpolating 2014 LIDAR points using LasTools and ArcGIS [57,58]. While a LIDAR-based CHM can be helpful in identifying vegetation heights in certain terrains, coastal wetland vegetation can cause errors and positive bias caused by litter and certain vegetation structures that can conceal the ground elevation [59]. As a result we only used the CHM to identify tall wetland vegetation, dead trees within the restored marsh, and forests.
Wetland inventories should include tidal mud flats, given their ecological importance of providing invertebrate prey resources for shorebirds and fishes, as well as their potential ability to sequester carbon [60]; therefore in the 2015 maps we used the low tide 2014 CIR image to define the exposed delta mudflats. The dead tree class identified trees that had grown within the levee boundary in the early 2000s and then subsequently died due to exposure to tidal inundation after the 2009 dike removal. The transitional classes occur within the buffer zone between riparian forest and emergent marsh and are commonly identified by the presence of cattails and a mix of other riparian and emergent marsh plants. Given the roughly defined continuum of the transition marsh between habitats, we used a post-restoration habitat zone dataset provided by the Tribe to define a generalized boundary of the transitional zone. In the final 2015 habitat classification, the transitional marsh extent was further refined with the addition of cattails as a transitional indicator species. We further delineated transitional forest and transitional marsh classes using a height threshold based on the CHM to separate tall woody vegetation from low shrub. Upland and riparian forest classes were identified using an elevation threshold, based on a 2010 DEM. Data produced for this study are available through ScienceBase.gov at https://doi.org/10.5066/F78G8JN3.

Accuracy Assessment
We assessed the accuracy of each classification using confusion matrices to determine user's accuracy, producer's accuracy, and overall accuracy. To improve error estimates of our habitat change accuracy assessments, we collected 50 additional hand-selected reference samples per image in areas along forest and wetland boundaries where change was most likely to occur [61]. After removing a subset of samples for classifier training, we had a total of 439 samples for validation of the 1957 classifications, 568 for the 1980 classifications, and 443 samples for the 2015 classification. The accuracy assessment was performed in ArcGIS, by comparing reference sample points to their predicted classes within each classification.
Accuracy assessments for multi-class change detection can be a time-consuming and complex process; therefore, we used a simplified, two-step methodology proposed by Congalton and Green [61] to determine accuracy. This method attempts to both determine the accuracy of the area mapped as change, and measure the success of identifying change. The first step in the change assessment is to determine the accuracy of predicted change by creating a single-date error matrix for the area mapped as having changed. This first step of the assessment allowed us to identify the accuracy of the class identified as changed, while also weighing the accuracy of the original classification. The second step of the change assessment is a simple change/no change matrix performed on the change map, to determine how well change was identified.

Change Detection
Change analysis within the Nisqually River Delta focused on three classes including wetland change, forest change, and photosynthetic functional type change between 1957, 1980, and 2015.
To perform the change analysis we created a Boolean (present/not present) raster for both the emergent wetland and forest classes for each year, and a three-class raster for the pathway classification (C3 plants, C4 plants, and other). Change analysis was completed in three series (1957 to 1980, 1980 to 2015, and 1957 to 2015), and while the resulting maps displayed net land cover changes, both net and gross change was calculated. In addition to the forest change detection map, we performed a change analysis on forest type. Based on the 1957 and 1980 forest extents, we overlaid riparian and upland forest class rasters to the 2010 DEM as had been done in the 2015 classification. While it is possible that historical terrain elevations differed from the 2010 DEM, we proceeded with the assumption that upland elevations remained relatively stable during our study period.

Habitat Classification
The 1957 habitat classification showed extensive areas of agriculture/grassland and forest as the dominant habitat type within the study area ( Figure 3). Emergent marsh extended north from the well-defined dikes to the delta flats, and along up McAllister Creek and the main channel of the Nisqually River. Overall accuracy for the habitat classification was 94.9% (Appendix A), with producer's and user's accuracies above 95% for almost all classes. The mudflat class resulted in the lowest accuracies of 70% and 77.7%, with mudflat identified as water in areas along the edge of deeper main sloughs and emergent marsh along exposed mudflat in first-order channels. In both cases, this error highlights the gradient of spectral ranges between water-wetland interfaces, which become particularly challenging in panchromatic images.
The 1980 habitat classification resulted in overall accuracy of 96.7% with all classes but developed land having both user's and producer's accuracy above 89%, The addition of spectral bands from the CIR image allowed for improved separation of water, vegetation, and agricultural land. Developed land increased along the southern edge of the managed freshwater wetland, coinciding with the building of Highway I-5 in 1967 and other nearby developments (Figure 4a).
The 2015 habitat classification also resulted in high overall accuracy (92.3%), and although this was lower than the 1957 and 1980 classifications, it was likely due to the increased number of classes in the 2015 classification ( Figure 4b).
Generally, most of the 2015 classes had accuracies above 80%; however, algae, aquatic bed, transition, and mudflat resulted in low class accuracies caused by confusion between these classes. Emergent marsh species most often misclassified as mudflat included pickleweed, Cotula coronopifolia (brass buttons), Lysimachia maritima (sea milkwart), and Scirpus cernuus (Seacoast bulrush). The majority of misclassified emergent marsh samples were located in the restored marsh where emergent marsh plants may be less established and where increased tidal exposure can cause spectral interference of mud and water ( Figure 5). The dynamic and somewhat ambiguous boundaries defining transitional marsh and emergent marsh help explain the confusion between the classes. The most common emergent marsh species classified as transitional marsh included Bolboschoenus maritimus (saltmarsh bulrush), Eleocharis acicularis (needle spikerush), Lyngbyei's sedge, and Juncus effusus (common rush). west accuracies of 70% and 77.7%, with mudflat identified as water in areas along the edge eper main sloughs and emergent marsh along exposed mudflat in first-order channels. In bot ses, this error highlights the gradient of spectral ranges between water-wetland interfaces, whic come particularly challenging in panchromatic images.  bands from the CIR image allowed for improved separation of water, vegetation, and agricultural land. Developed land increased along the southern edge of the managed freshwater wetland, coinciding with the building of Highway I-5 in 1967 and other nearby developments (Figure 4a).
The 2015 habitat classification also resulted in high overall accuracy (92.3%), and although this was lower than the 1957 and 1980 classifications, it was likely due to the increased number of classes in the 2015 classification (Figure 4b).  Generally, most of the 2015 classes had accuracies above 80%; however, algae, aquatic bed, transition, and mudflat resulted in low class accuracies caused by confusion between these classes. Emergent marsh species most often misclassified as mudflat included pickleweed, Cotula coronopifolia (brass buttons), Lysimachia maritima (sea milkwart), and Scirpus cernuus (Seacoast bulrush). The majority of misclassified emergent marsh samples were located in the restored marsh where emergent marsh plants may be less established and where increased tidal exposure can cause spectral interference of mud and water ( Figure 5). The dynamic and somewhat ambiguous boundaries defining transitional marsh and emergent marsh help explain the confusion between the classes. The most common emergent marsh species classified as transitional marsh included Bolboschoenus maritimus (saltmarsh bulrush), Eleocharis acicularis (needle spikerush), Lyngbyei's sedge, and Juncus effusus (common rush).  The photosynthetic pathway classifications highlight the challenges of plant type identification through historical imagery. The 1980 and 2015 classifications identified C3 plants as the dominant marsh type with C4 plants occurring along the outer marsh banks and in low marsh areas on the eastern edge of the marsh. Overall accuracy for the 1980 classification was 79.1%, with user's and producer's accuracy was between 80% and 89% for C3 plants and C4 plant accuracies at 74.6% and 58.7%. The 2015 classification achieved slightly higher overall accuracy of 81.5%, with all class accuracies above 70%. In both the 1980 and 2015 classifications C3 plant accuracies remained higher than C4 plants, and C4 plants were omitted from the maps more often than C4 plants were included erroneously.

Species Classification
The 1980 species classification showed large extents of the 'jaumea/carex/salt grass' distributed throughout the marsh, with the 'carex/hairgrass/festuca' class occurring in the Animal site, while 'salt grass' was found in larger areas on the outer eastern edges of the delta, near the Nisqually River ( Figure 6). The 1980 species classification had an overall accuracy of 81.6% for all classes, however, when the six dominant species were assessed alone, accuracy dropped to 62.1% (Appendix A). The individual accuracies of each species class varied greatly; however, the 'salt grass' and 'jaumea/carex/salt grass' classes had user's and producer's accuracies above 60%.  The 2015 species map was classified by single species classes, showing a more detailed distribution of classes (Figure 7). 'Pickleweed' was dominant in the restored marsh, while 'salt grass' was distributed throughout the reference sites. Similar to the 1980 classification, patches of 'carex' are shown grouped throughout the Animal site. The overall accuracy for the 2015 species classification was 81% for all classes, as a result of very low producer's and user's accuracies for most species classes (Appendix A). 'Salt grass' was the most successfully identified species, with 86.8% user's accuracy and 71.9% producer's accuracy. Often the confusion between species occurred between nonsalt grass classes and the 'other emergent marsh' class. The 2015 species map was classified by single species classes, showing a more detailed distribution of classes (Figure 7). 'Pickleweed' was dominant in the restored marsh, while 'salt grass' was distributed throughout the reference sites. Similar to the 1980 classification, patches of 'carex' are shown grouped throughout the Animal site. The overall accuracy for the 2015 species classification was 81% for all classes, as a result of very low producer's and user's accuracies for most species classes (Appendix A). 'Salt grass' was the most successfully identified species, with 86.8% user's accuracy and 71.9% producer's accuracy. Often the confusion between species occurred between non-salt grass classes and the 'other emergent marsh' class.
was distributed throughout the reference sites. Similar to the 1980 classification, patches of 'carex' are shown grouped throughout the Animal site. The overall accuracy for the 2015 species classification was 81% for all classes, as a result of very low producer's and user's accuracies for most species classes (Appendix A). 'Salt grass' was the most successfully identified species, with 86.8% user's accuracy and 71.9% producer's accuracy. Often the confusion between species occurred between nonsalt grass classes and the 'other emergent marsh' class.

Forest Change
The 1957 classification identified 1430.8 ha of forest within our study site, which decreased by 8.6% to 1307.4 ha in 1980 and 1229.6 ha in 2015 (Table 3). Although total net forest loss decreased over time, gross forest loss increased, resulting in 537 ha of forest disturbance for the study period (Table 4). This change occurred in both upland and riparian areas, with forest loss most commonly resulting from land development, agriculture, and river channel movement. The shifting river channel was distinctly visible in the forest change maps, and although each map shows forest loss associated with the new river channel, forest regrowth also filled in old channels. Forest increased from 1957 to 1980 in areas along the Nisqually River and near McAllister Creek near Highway I-5 ( Figure 8). Overall accuracy for the forest change maps were all above 87% (Appendix B).

Forest Change
The 1957 classification identified 1430.8 ha of forest within our study site, which decreased by 8.6% to 1307.4 ha in 1980 and 1229.6 ha in 2015 (Table 3). Although total net forest loss decreased over time, gross forest loss increased, resulting in 537 ha of forest disturbance for the study period (Table 4). This change occurred in both upland and riparian areas, with forest loss most commonly resulting from land development, agriculture, and river channel movement. The shifting river channel was distinctly visible in the forest change maps, and although each map shows forest loss associated with the new river channel, forest regrowth also filled in old channels. Forest increased from 1957 to 1980 in areas along the Nisqually River and near McAllister Creek near Highway I-5 ( Figure 8). Overall accuracy for the forest change maps were all above 87% (Appendix B).

Riparian and Upland Forest Change
The forest type change map identified a 7.7% net loss (Table 5) of riparian forest between 1957 and 1980 due to conversion to agriculture in the lowlands between McAllister Creek and the Nisqually River, while the minimal loss from river channel migration was negated by regrowth of riparian vegetation in the old river channel (Figure 9). From 1980 to 2015, riparian forest experienced a 2.3% increase located within the Refuge boundaries. Upland forest was minimally impacted between 1957 and 1980 with only a 0.9% loss caused by forest removal for quarry expansion, yet between 1980 and 2015 a more significant decline of 8.2% of upland forest occurred due to land development. riparian vegetation in the old river channel (Figure 9). From 1980 to 2015, riparian forest experienced a 2.3% increase located within the Refuge boundaries. Upland forest was minimally impacted between 1957 and 1980 with only a 0.9% loss caused by forest removal for quarry expansion, yet between 1980 and 2015 a more significant decline of 8.2% of upland forest occurred due to land development.

Wetland Change
The wetland change analysis showed an initial decrease in tidal marsh from 1957 to 1980, followed by wetland gains from 1980 to 2015 occurring within the former levee boundaries, indicating an expansion of tidal wetlands and a transition from agricultural land to emergent marsh within the restoration sites. This 54.2% net increase in emergent marsh in 2015 resulted in 342.7 ha of tidal wetland vegetation (Table 6). Total wetland conversion occurred on 271.5 ha, and despite trends of initial net decrease, then increase, wetland loss and gain occurred throughout the study period (Table 7).

Photosynthetic Pathway Change
The photosynthetic pathway change map showed transitions of C3 and C4 plants in the emergent marsh ( Figure 11). Conversion from salt grass dominated C4 vegetation to C3 vegetation is visible largely along the outer edges of the low marsh and along channels where erosion may occur, while the transition from C3 to C4 plants can be observed in areas of higher marsh.

Photosynthetic Pathway Change
The photosynthetic pathway change map showed transitions of C3 and C4 plants in the emergent marsh ( Figure 11). Conversion from salt grass dominated C4 vegetation to C3 vegetation is visible largely along the outer edges of the low marsh and along channels where erosion may occur, while the transition from C3 to C4 plants can be observed in areas of higher marsh.

Photosynthetic Pathway Change
The photosynthetic pathway change map showed transitions of C3 and C4 plants in the emergent marsh ( Figure 11). Conversion from salt grass dominated C4 vegetation to C3 vegetation is visible largely along the outer edges of the low marsh and along channels where erosion may occur, while the transition from C3 to C4 plants can be observed in areas of higher marsh.

Discussion
The classification of historical habitats resulted in high overall accuracies, demonstrating the potential for successful mapping of coastal wetlands at the habitat scale using both multispectral and panchromatic historical imagery. The addition of supplemental spectral and elevation data to the 2015 image enabled a more detailed classification of habitat types, yet even with minimal spectral band and texture information, the historical classifications successfully distinguished basic habitats.
Along the water/marsh boundary in the 1957 classification and in the restored marsh in the 2015 classification, some confusion of classes occurred between emergent marsh vegetation and mudflat, likely due to tidal influence and standing water. This can be expected for low marsh species with limited patch extent, as misclassifications commonly occur with the presence of water and mud interfering with the normal spectral response of a plant [23]. These types of misclassification in single image classification can present a challenge with change analysis. This was apparent in the photosynthetic functional type change maps series where C3 and C4 marsh plants showed inconsistent and reversing trends, which we attributed to the low accuracies in the 1957 classification caused by limited spectral information rather than actual functional type changes. For the well delineated forest and emergent marsh classes, change maps allowed us to measure habitat area and percent change over time to successfully quantify vegetation trajectories. Furthermore, the change maps show the morphological changes of channel movement and erosion, changes in vegetation cover, and land use/land cover change, which can be indicators of the effects of watershed processes, restoration, and local land management.

Limitations and Error
Data availability and image types pose noteworthy challenges for historical wetland mapping and change analysis. In addition to limited availability and low frequency of image collection, historical aerial photographs tend to be panchromatic or natural color resulting in insufficient data to inform highly-detailed classifications. To resolve some of the challenges of limited data availability, image derivatives such as texture may help improve classification results. More recent advancements in imagery and data include LIDAR and hyperspectral sensors that can improve wetland vegetation mapping; however, the cost of acquiring and processing these data can be a barrier to use [24]. As demonstrated by the accuracies of our 2015 species map, even high-resolution multispectral imagery may not be sufficient for species vegetation mapping in highly heterogeneous tidal wetlands. The complexity of change analysis in wetland promotes further challenges caused by seasonal and tidal variability in imagery. High tides captured in the 1980 image obscured the majority of mudflat, resulting in insufficient area to assess accuracy of the class. While historical tidal data may help inform classification and change analysis, in some cases this data may not be available. A visual comparison of the 1980 images to other images at known tide stages helped inform us that tidal inundation of the emergent marsh had not occurred during image acquisition. With insufficient tidal data, however, the impact of tidal inundation in the imagery can be a concern and the suitability of the historical image for tidal wetland analysis should be considered. Additionally, reference data can be challenging to come by, and when available, the missing accuracy and processing information can introduce uncertainty, as we found with the 1975 vegetation map. With such challenges and limitations of both historical data availability and potential issues with tides, special consideration should be made for image selection and potential introduction of errors.
Assessing accuracy for change detection between two classifications can be a time-consuming and costly task, due to the increase in complexity and potential for error propagation. Error sources from a single date classification can occur during image acquisition and processing, data analysis, error assessment, and even in the production of final maps. Because change detection requires two or more dates, the addition of the second classification compounds the error. Finally, change detection analysis can result in attribute error and positional error, which results in confusion of whether actual or apparent change has occurred. As a result of this complexity, many studies fail to assess change accuracy [61]. While a full error assessment for change detection in the refuge was not feasible, our two-step assessment was able to determine the accuracy of the areas identified as change and how accurate those areas of change were classified.

Wetland Change
Wetland change maps provide insight into long term trajectories and trends in wetland health, but these trajectories can also be indicators of the influence of tidal marsh restoration. In 2005, the Nisqually Comprehensive Conservation Plan estimated intertidal estuarine habitat in the Nisqually River Delta to be around 2511.8 ha in 1878, with a 19% loss by 2005. Salt marsh within the estuarine habitat declined from 590.0 to 272.8 ha, totaling a substantial 54% loss due to diking, fill, and channel movement [34]. Our 1980 estimates of pre-restoration salt marsh area were 50.6 ha less than the US Fish and Wildlife Service 2005 [34] estimates; however, this difference can be explained by the small restoration projects that occurred in the early 2000s after the 1980 classification. Despite the shift away from agriculture in the 1960s to freshwater wetland under Refuge management in the 1970s, we only identified a small 6.4% net loss in tidal wetland from 1957 to 1980 related to the loss of tidal wetland in the Phase II site. The final shift to tidally-influenced wetland occurred with the 2009 restoration, highlighted in the 1980-2015 change map with a 54% net increase in emergent marsh. During this time wetland areas increased from 222.2 to 342.4 ha, with most change located in the restored marsh sites.
Long-term success of the tidal wetland restoration requires marsh recovery, which is directly tied to elevation gain, and in turn, the establishment of vegetation. An estimated 1 m of subsidence has occurred at the Nisqually salt marsh since diking and drainage of the wetland [33], therefore sufficient vertical accretion is necessary to maintain an elevation where vegetation can develop without overexposure to tidal inundation. The impacts of the 2009 restoration on the emergent marsh are reflected in the 1980-2015 change map by both the regrowth of salt marsh vegetation and presence of mudflat in the restored marsh. As expected, the reintroduction of tidal exposure in the restored site resulted in a dramatic conversion from freshwater wetland to mudflat. The dieback of freshwater plants and litter buildup initially helped increase marsh elevation through sediment delivery and capture. After the dike removal, tidal exposure to newly accessed channels enabled sediment transport and deposition in the restored marsh, while at the same time propagating and extending erosion into inland channels [62]. To better understand sediment supply connectivity and estuary turbidity, water modeling and monitoring occurred between 2010 and 2011. Through modeling simulations it was determined that suspended sediments transported by the river were often deposited into channels before reaching the marsh, and as a result, sediment delivery into the restoration marsh occurred only during high tides [63]. Stable emergent vegetation may take some time to develop at the restored marsh, yet once established, sediment accretion should be accelerated due to stabilization and increased friction from vegetation [64]. Vegetation recovery after marsh disturbances often begins with growth from edge patches into bare areas, with salt grass and pickleweed recovering quickest [65]. Based on vegetation elevation modeling, Monden [62] predicted that the establishment of plants would begin in patches filling in on the eastern side of the restoration, with minimal vegetation developing on the western side due to low elevations. In 2010, half of the restored marsh had sufficient height for the establishment of vegetation, with most of it occupied by low marsh plants [40]. Two years after the restoration, colonization of tidal marsh vegetation was mostly limited to pickleweed (Salicornia pacifica) and sand-spurry (Spergularia canadensis), with only 20% of the monitoring quadrats containing emergent plants [66]. In line with these predictions and observations the restored site on our 2015 habitat classification identified limited and non-continuous growth of emergent marsh vegetation on the eastern edge of the restored wetland, with minimal development on the western mudflats. This may suggest that the restoration is on a slow trajectory towards reestablishment of vegetation and wetland recovery.

Edge Erosion
Tidal marsh change mapping not only identifies shifts in estuarine marsh habitats, but also can be used to indicate morphological changes, such as edge erosion. Despite gains in wetland regrowth in the restored marsh, marsh loss was apparent most commonly occurring along tidal meanders and slough edges. In the Nisqually River Delta, the transition between tidal flat and vegetated marsh occurs in low gradient slopes caused by tidal sediment accretion and in shear eroded banks along Nisqually River, McAllister Creek, and where strong tidal currents exist [36]. Wetland loss along the main Nisqually River channel and sloughs are distinct in both change maps, but to a larger extent between 1957 and 1980.
Horizontal instability in coastal wetlands is less understood than vertical stability, as lateral edge erosion and expansion of marshes can be in constant movement [67,68]. Lateral marsh scarp boundary movement is a complex function of waves, sediment transport and deposition, vegetation, wind, and tidal oscillation [69]. Although established vegetation is helpful for edge stabilization during gradual sea-level and tidal changes, soil type is a greater contributor to stability in more extreme events, such as breaking waves [70]. Both vertical and horizontal movement can contribute to loss of carbon storage capacity, as carbon stored in peats, soils, and biomass is removed due to erosion.

Photosynthetic Functional Type Shifts
Photosynthetic functional type shifts showed changes in dominance of C3 and C4 marsh plants with some stratification by high marsh and low marsh, reflecting the influence of marsh elevation on the placement of C3 and C4 marsh vegetation within the emergent wetlands. In addition, photosynthetic functional type change mapping can aid in research to evaluate how carbon sources from different habits provide trophic support to juvenile Chinook food webs using stable isotope analysis (SIA). Similarly, these maps can be combined with SIA of soil cores to identify the isotopic signatures of primary producers contributing to peat formation in both the restored and natural marshes.

Watershed-Scale Changes
Watershed-scale changes in forest and riparian habitat can be indicators of development and geomorphic changes in river channels as a response to both anthropogenic and natural drivers. Both upland and riparian forest in the Pacific Northwest play an important role for carbon sequestration [11] and as transitional zones for food webs and carbon movement [12]. The assessment of forest type change showed greater loss of riparian forest and limited upland loss between 1957 and 1980, and a slight increase in riparian forest and large decline in upland forest from 1980 to 2015. Much of the early riparian forest loss can be attributed to agricultural and grassland expansion, as well as shifting of the river channel. The period from 1980 to 2015 showed riparian forest conversion decreasing despite development near highway I-5 and the quarry. At the same time, riparian forest near the transitional zone of the Refuge accounted for a large portion of riparian regrowth. Upland forest loss after 1980 reflects an increase in housing development and conversion to open space on both the eastern and western bluffs of the Nisqually River Delta, which accounted for the majority of the 8.3% loss of upland forest. This regional development reflects a trend of population growth in Thurston County with an increase of 55,049 residents in 1960 to 269,536 in 2015 [71] resulting in a 389% increase in population in Thurston County and 162% increase in Pierce County. With continued projections of growth for these counties, direct watershed-scale drivers related to development and land cover disturbances can be expected.
Tracking both historical change within coastal wetland systems and the dynamic interface with watersheds is necessary for understanding the net ecosystem carbon balance (NECB) and tidal wetland sustainability with expected future sea level rise. The quantification of net changes on a landscape, while informative for habitat trajectories and understanding spatial extents of land alterations, results in an underestimation of total land change. When accounting for carbon storage or habitat stability, gross change is a necessary metric to consider dynamic changes on the landscape by accounting for long-term disturbance over time impacting wetland carbon storage potential, forest age, and sediment exposure [72]. The NECB of a coastal marsh is determined by examining historical extents and alterations within coastal marshes, assessing atmospheric carbon exchanges, and measuring tidal transportation of dissolved organic carbon, dissolved inorganic carbon, and sediments [73].
Future consideration of coastal marsh resilience depends on rates of accelerated sea-level rise and availability of sediments. Using WARMER, a wetland accretion model, Thorne et al. [74] examined the impacts of sea-level rise over the next century in the Pacific Northwest, projecting an eventual transition to un-vegetated tidal flat at the Nisqually River Delta under mid to high sea-level rise due to the limited sediment supply and accretion rates. In addition to rising sea levels, earthquake-related subsidence [75] and increased groundwater extraction could further decrease the wetland elevation in the future. Anthropogenic disturbances and climatic changes are expected to reduce carbon storage in wetlands; yet these same climate changes through extreme weather events may increase riverine transport and deposition of carbon and sediments [10]. Anthropogenic changes to the environment at the watershed scale and direct interventions at the local wetland scale may ultimately determine the future sustainability and resilience of the Nisqually marsh.
Habitat classification and change mapping at multiple scales can provide insight into landscape-scale changes. At the wetland-scale the tidal marsh responded to the reintroduction of tidal flow during the 2009 restoration with an initial vegetation dieback and conversion to mudflat, but also with eventual emergent marsh regrowth and reestablishment of plants in higher marsh elevations. Extending beyond the emergent marsh, up Nisqually River Delta we see consistent edge bank erosion in the tidal marsh and channel migration along the river occurring over the last 60 years of our study period. Along the river corridor, floodplain, and upland bluffs we see a steady decline in forest occurred during our study period beginning with riparian forest conversion for agricultural purposes and later shifting to upland forest loss for land development. Despite these trends in forest loss, more recent growth of riparian forest habitats have occurred in the transitional forested riverine tidal zones of the Refuge. By interpreting the wetland-to watershed-scale transformations over time, we can shed light on the integrated drivers of these changes to infer future trajectories of the Nisqually River Delta.
The long-term outlook for the Nisqually tidal marsh depends on continued sediment supply from the watershed and the wetland response to sea level rise. With expected changes in climate, along with population and industry growth, alterations to sediment delivery to the Nisqually estuary can be anticipated. Along river corridors, discharge regulates sediment transport from terrestrial systems, while climate ultimately dictates the rates of transport with both normal flows and extreme flood events [10]. Such flood events have potential to increase erosion, bringing greater amounts of sediment and carbon to the delta. While sediments captured at dams are still unavailable for wetland deposition, land disturbances from agriculture and development may continue to contribute to sedimentation of the coastal wetland. Bank erosion within the delta results in marsh loss, but in turn provides a redistribution of sediments into the tidal system. The reestablishment of vegetation in the restored marsh, although slow, may suggest movement towards future emergent marsh stability. With post-restoration vertical accretion rates of 2.5 mm/year in the restoration site between 2010 and 2014, it is possible that long-term recovery of the restoration site could match conditions of the nearby natural wetland sites under current conditions [47].

Conclusions
Historical mapping of wetlands can help inform long-term tidal wetland and land-cover trajectories. Tracking changes in coastal habitats, wetland emergent marsh, and forests is just a small step towards understanding long-term trends and short-term responses of carbon stocks in coastal wetlands. Furthermore, mapping marsh vegetation extents can improve aboveground biomass accounting and help inform understanding of channel erosion and sediment transport in soil carbon lateral flux, while species vegetation mapping can provide insight into high and low marsh extents for sea-level rise resiliency. Remotely-sensed marsh losses and gains in the restoration site were generally consistent with previous studies, demonstrating the feasibility of historical change mapping as a method for understanding past wetland conditions and potential marsh trajectories in the future. With the extent of marsh channel erosion occurring along the Nisqually River, further study is necessary to understand the movement and deposition of eroded bank sediment as a sediment source for the restored marsh. Additionally, soil core analysis of vertical accretion within the restored and reference sites may help quantify tidal marsh sediment deposition from normal flow and flood events preand post-restoration.  Table A2. Two-step accuracy assessment for forest change.

Appendix A
Step 1: Single Date Confusion Matrix Step  Table A3. Two-step accuracy assessment for wetland change.
Step  Table A4. Two-step accuracy assessment for photosynthetic functional type change. Step