Monitoring Tamarix Changes Using WorldView-2 Satellite Imagery in Grand Canyon National Park, Arizona

: Remote sensing methods are commonly used to monitor the invasive riparian shrub tamarisk ( Tamarix spp. ) and its response to the northern tamarisk beetle ( D. carinulata ), a specialized herbivore introduced as a biocontrol agent to control tamarisk in the Southwest USA in 2001. We use a Spectral Angle Mapper (SAM) supervised classiﬁcation method with WorldView-2 (2 m spatial resolution) multispectral images from May and August of 2019 to map healthy tamarisk, canopy dieback, and defoliated tamarisk over a 48 km segment of the Colorado River in the topographically complex Grand Canyon National Park, where coarse-resolution satellite images are of limited use. The classiﬁcations in May and August produced overall accuracies of 80.0% and 83.1%, respectively. Seasonal change detection between May and August 2019 indicated that 47.5% of the healthy tamarisk detected in May 2019 had been defoliated by August 2019 within the WorldView-2 image extent. When compared to a previously published tamarisk map from 2009, derived from multispectral aerial imagery, we found that 29.5% of healthy tamarisk canopy declined between 2009 and 2019. This implies that tamarisk beetle impacts are continuing to accumulate even though land managers have noted the presence of the beetles in this reach of the river for 7 years since 2012.


Introduction
Tamarisk (Tamarix spp.) is an invasive, Eurasian shrub that has substantially altered riparian ecosystems throughout the Southwestern United States. In the early 1820s, it was introduced to the United States for ornamental purposes [1,2] and four species (T. ramossisima, T. chinensis, T. parviflora, and T. aphylla) have since become invasive throughout arid Southwestern riparian riverbanks and floodplains [2]. Invasion of tamarisk has been aided by the initialization of regulated flow regimes and lowered abundance of native vegetation along floodplains due to dams built along riparian corridors throughout the western US [1].
Tamarisk can decrease groundwater supply [3] and alter hydrologic processes. Estimates produced nearly two decades ago indicated a potential cost of this in the Western US to be 1.4-3.0 billion m 3 of water per year [4], although tamarisk has since been shown to have low to moderate rates of water use relative to native species depending on stand age, leaf area, density, and access to groundwater [3,[5][6][7][8]. Tamarisk can procure groundwater and lower the water table via a deep tap root [3], transformation of soil chemistry via soil salinization and disruption of mycorrhiza mutualisms [9], and high rates of seed production [10]. Tamarisk also has physiological traits that allow it to thrive in disturbed environments [11]. Additionally, rapid hybridization of T. ramosissima and T. chinensis in the US has increased the genetic variability of this hybrid, potentially increasing fitness of it in US riparian corridors [12]. Tamarisk Our study goal was to separately classify three varying stages of tamarisk: (1) healthy tamarisk, (2) canopy dieback, and (3) defoliated tamarisk using the WorldView-2 multispectral satellite sensor with 2 m spatial resolution between 18 May 2019 and 2 August 2019. We hypothesized that the Worldview-2 classification would accurately detect tamarisk distribution, canopy dieback, and defoliation and provide a valuable data source for continued tamarisk monitoring. Canopy dieback in our study is defined as repeated observation of no green foliage on tamarisk, where tamarisk canopy does not appear to recover to a healthy state between observation dates.

Field Observations
This study was conducted along a 48.5 km reach of the Colorado River within Grand Canyon National Park, Arizona, United States ( Figure 1). While the tamarisk beetle first appeared in Grand Canyon in 2009, they were not detected in our study reach until 2012. Field data were collected repeatedly for this study during one trip taken down the river each month from 16 to 30 May, from 13 to 22 June, and from 8 to 17 July 2019. Observations were made about every 0.4 river km at tamarisk stands with total canopy area greater than approximately 16 m 2 .
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 18 Index (NDVI) [14,16,38,39] or Enhanced Vegetation Index (EVI) [6,40,41], or by detecting changes in supervised classifications [41]. However, we are aware of no studies that have taken advantage of supervised classification methods to discern these different stages of tamarisk within a single growing season.
Our study goal was to separately classify three varying stages of tamarisk: 1) healthy tamarisk, 2) canopy dieback, and 3) defoliated tamarisk using the WorldView-2 multispectral satellite sensor with 2 m spatial resolution between 18 May 2019 and 2 August 2019. We hypothesized that the Worldview-2 classification would accurately detect tamarisk distribution, canopy dieback, and defoliation and provide a valuable data source for continued tamarisk monitoring. Canopy dieback in our study is defined as repeated observation of no green foliage on tamarisk, where tamarisk canopy does not appear to recover to a healthy state between observation dates.

Field Observations
This study was conducted along a 48.5 km reach of the Colorado River within Grand Canyon National Park, Arizona, United States ( Figure 1). While the tamarisk beetle first appeared in Grand Canyon in 2009, they were not detected in our study reach until 2012. Field data were collected repeatedly for this study during one trip taken down the river each month from 16 to 30 May, from 13 to 22 June, and from 8 to 17 July 2019. Observations were made about every 0.4 river km at tamarisk stands with total canopy area greater than approximately 16 m 2 . Ground-based observations were made using methods of Jamison et al. [24]. At tamarisk stands or continuous clusters of tamarisk with total canopy area greater than 16 m 2 , canopy dieback was visually estimated in 10 percent increments (e.g., 0-10%, 10-20%, etc.) by comparing the proportion of dead branches to all branches. Once percent of dead branches was recorded, the proportion of live-but defoliated-branches to all living branches was estimated and recorded, in 10 percent increments as well. Defoliated tamarisk can be identified by brown leaves, or red stems, which would normally be covered in green foliage during the growing season from April through October. Dead branches can Ground-based observations were made using methods of Jamison et al. [24]. At tamarisk stands or continuous clusters of tamarisk with total canopy area greater than 16 m 2 , canopy dieback was visually estimated in 10 percent increments (e.g., 0-10%, 10-20%, etc.) by comparing the proportion of dead branches to all branches. Once percent of dead branches was recorded, the proportion of live-but defoliated-branches to all living branches was estimated and recorded, in 10 percent increments as well. Defoliated tamarisk can be identified by brown leaves, or red stems, which would normally be covered in green foliage during the growing season from April through October. Dead branches can be further distinguished from defoliated branches based on their silver color (Figure 2), whereas defoliated stems after leaf drop are a reddish color [27].
Field-based spectral signatures were also collected using an ASD Fieldspec 3 (https:// www.malvernpanalytical.com/en accessed on 10 May 2019) for tamarisk as well as riparian shrub species that commonly grow adjacent to tamarisk in Grand Canyon: seepwillow (Baccharis spp.), arrowweed (Pluchea sericea), acacia (Acacia greggii), and mesquite (Prosopis glandulosa). Field spectra of tamarisk were taken of pure green and healthy foliage, as well as varying levels of canopy dieback, including pure dead branches, and mixed healthy foliage and dead branches. Each spectrum was an average of 10 scans acquired at the canopy-level close to solar noon. be further distinguished from defoliated branches based on their silver color (Figure 2), whereas defoliated stems after leaf drop are a reddish color [27]. Field-based spectral signatures were also collected using an ASD Fieldspec 3 (https://www.malvernpanalytical.com/en accessed on 10 May 2019) for tamarisk as well as riparian shrub species that commonly grow adjacent to tamarisk in Grand Canyon: seepwillow (Baccharis spp.), arrowweed (Pluchea sericea), acacia (Acacia greggii), and mesquite (Prosopis glandulosa). Field spectra of tamarisk were taken of pure green and healthy foliage, as well as varying levels of canopy dieback, including pure dead branches, and mixed healthy foliage and dead branches. Each spectrum was an average of 10 scans acquired at the canopy-level close to solar noon.

WorldView-2 Satellite Images and Analysis
WorldView-2 images were acquired for the reach of the river that extends from Parashant Creek at river km 315 to the Diamond Creek at river km 363 ( Figure 1). River kilometers correspond to the river distance downstream the Colorado River from Lees Ferry, AZ, a standardized distance system used by researchers focused on the Grand Canyon (https://www.usgs.gov/centers/sbsc/gcmrc accessed on 10 May 2019). Worldview-2 is a commercial satellite that collects data in eight spectral bands. Two WorldView-2 images were used in this study: one image acquired 18 May 2019, and a second acquired 2 August 2019. The May 2019 image was acquired during the time ground-based observations were taken in that month. The August 2019 image was acquired 10 days after the last groundbased observations were taken in July, and corresponds with seasonal patterns of highest defoliation, which are expected to be similar to what was observed in July.

WorldView-2 Satellite Images and Analysis
WorldView-2 images were acquired for the reach of the river that extends from Parashant Creek at river km 315 to the Diamond Creek at river km 363 ( Figure 1). River kilometers correspond to the river distance downstream the Colorado River from Lees Ferry, AZ, a standardized distance system used by researchers focused on the Grand Canyon (https://www.usgs.gov/centers/sbsc/gcmrc accessed on 10 May 2019). Worldview-2 is a commercial satellite that collects data in eight spectral bands. Two WorldView-2 images were used in this study: one image acquired 18 May 2019, and a second acquired 2 August 2019. The May 2019 image was acquired during the time ground-based observations were taken in that month. The August 2019 image was acquired 10 days after the last ground-based observations were taken in July, and corresponds with seasonal patterns of highest defoliation, which are expected to be similar to what was observed in July.
Image pixels were offset between the two 2019 WorldView-2 images. The May 2019 image was, therefore, coregistered to the August 2019 image and then both images were registered to the May 2013 USGS 20 cm resolution aerial imagery [42,43]. This coregistration allowed for direct comparison between the WorldView-2 images and USGS 20 cm imagery, which was used by Bedford et al. [13,14] to generate a 2013 tamarisk map of the Grand Canyon. The same area of small cloud was masked out from both images.
Both WorldView-2 images were classified using the Spectral Angle Mapper (SAM) supervised classification algorithm [44] to map three target classes: (1) healthy tamarisk, (2) canopy dieback, and (3) defoliated tamarisk. The classifications were trained using four canopy-level spectra from the ASD field spectra collected in May 2019: (1) healthy tamarisk leaf spectra, (2) canopy dieback spectra, (3) mixed healthy leaf and canopy dieback spectra, and (4) defoliated tamarisk spectra. A fifth image-derived training spectra for bare ground (n = 100 pixels) was used to accurately separate canopy dieback from bare ground, since it was found that these two endmembers have similar spectra. The original ASD spectra were collected in 1250 spectral bands from 350 to 2500 nm, which were subset to match the center of the wider bandwidths of the eight spectral bands in the WorldView-2 image (Figure 3). and 4) defoliated tamarisk spectra. A fifth image-derived training spectra for bare ground (N=100 pixels) was used to accurately separate canopy dieback from bare ground, since it was found that these two endmembers have similar spectra. The original ASD spectra were collected in 1250 spectral bands from 350 to 2500 nm, which were subset to match the center of the wider bandwidths of the eight spectral bands in the WorldView-2 image ( Figure 3). Ground-based hyperspectral training spectra for the four tamarisk classes, spectrally subset from the ASD Fieldspec 3 spectroradiometer to match the WorldView-2 spectral resolution.
Five training classes were classified using the five spectra described above (Table 1). Subsequently, the bare ground class was deleted, and the mixed healthy and dead branch and canopy dieback classes were combined into a single, canopy dieback class. Those Five training classes were classified using the five spectra described above (Table 1). Subsequently, the bare ground class was deleted, and the mixed healthy and dead branch and canopy dieback classes were combined into a single, canopy dieback class. Those postprocessing steps produced a 3-class classification map, which we edited manually to omit extra canopy dieback that was misclassified as bare ground. Finally, a sieve operation was used to improve the final 3-class classification map by reducing spurious individual pixels within groupings of each class. The post-processing steps were performed because of the applied nature of the science, where our goal is to provide resource managers such as the National Park Service with maps that are as correct as possible without obvious misclassifications that can be fixed through these post-processing methods. When using the ground-based observations from May to train the image classification for May 2019, we had multiple locations that could be used for the three target classes: (1) healthy tamarisk, (2) canopy dieback, and (3) defoliated tamarisk. However, when using the ground-based observations from August to train the August 2019 image classification, it was not possible to accurately discern defoliated tamarisk from tamarisk canopy dieback owing to pixels in which these endmembers were mixed. Therefore, the defoliated tamarisk and canopy dieback classes were combined into a single class for August to create a final classification map representing only two classes: (1) healthy tamarisk, and (2) "beetleimpacted" tamarisk.

Accuracy Assessment
To assess the accuracy of these maps, the ground-based observations were converted to a binary format based on thresholds described below for each of the SAM target classes. The SAM target classes for the May classification were healthy tamarisk, canopy dieback, and defoliated tamarisk. The SAM target classes for the August classification were healthy tamarisk and beetle-impacted tamarisk. For May, canopy dieback was defined as >50% dead branches, defoliation was defined as >50% defoliated branches, and healthy tamarisk was defined as >90% healthy tamarisk branches. For August, beetle-impacted tamarisk was defined as >50% defoliation, and healthy tamarisk was defined as >90% healthy tamarisk branches. Applying these parameters to the ground observation data, we identified 57 healthy tamarisk samples, 45 canopy dieback samples, and three defoliated tamarisk samples in May (Table 2), and 17 healthy tamarisk samples and 48 beetle-impacted tamarisk samples in August (Table 3).

Change Detection
We then detected changes between the two classification maps to determine where healthy tamarisk changed to the beetle-impacted tamarisk class between May and August of 2019. We also compared our data to that of Bedford et al. [13,14]. Using airborne data, Bedford et al. [13,14] created a 2009 healthy tamarisk map and a 2013 beetle-impacted map of our study area. The 2009 healthy tamarisk map was generated before tamarisk beetle activity in Grand Canyon, while the 2013 beetle-impacted map was the result of a change detection between May 2009 and May 2013 [13,14]. Both aerial images were from late May, so the 2013 aerial image could have contained canopy dieback and some seasonal defoliation. We compared the 2009 healthy tamarisk map and the 2013 beetle-impacted map with our May, 2019 WorldView-2 map to determine: (1) how much tamarisk mapped in 2009 has died back, and (2) how spatial patterns and distribution of canopy dieback have changed since 2013. The overall accuracy in the 2009 map [13,14] ranged from 70% to 79% in the three map quads that overlap our WorldView-2 images. The May, 2019 WorldView-2 classification map was clipped to the Bedford et al. [13,14] extent of tamarisk to ensure that the spatial extent of comparison was the same between the datasets. The Bedford et al. [13,14] maps were resampled to match the spatial resolution of the WorldView-2 classification map, and changes were detected between the compared images. Tamarisk pixels that changed from healthy in 2009 or 2013 to dead branches in 2019 were identified as tamarisk canopy dieback between these years. Pixels that changed from beetle-impacted tamarisk to healthy tamarisk between 2013 and 2019 were identified as having recovered between these years, either from a state of canopy dieback or seasonal defoliation in 2013. The area without change was also reported for each of these two change detections.  (Figure 4). In May 2019, defoliation was only observed at a few locations. In June, the number of defoliation observations increased, as did the average proportion of defoliated branches in each observation. In July, defoliation of tamarisk by the beetle was observed almost everywhere at a higher proportion per tamarisk stand.

Tamarisk Detection
Canopy dieback was observed to some degree almost everywhere along the study reach, with only small patches of tamarisk with 0% dieback ( Figure 5). Tamarisk canopy dieback from previous events and prior years was most easily observed in May 2019 both on the ground and in the imagery, when the lack of tamarisk beetle defoliation reduced the need to discern defoliated branches from dead branches.  Canopy dieback was observed to some degree almost everywhere along the study reach, with only small patches of tamarisk with 0% dieback ( Figure 5). Tamarisk canopy dieback from previous events and prior years was most easily observed in May 2019 both on the ground and in the imagery, when the lack of tamarisk beetle defoliation reduced the need to discern defoliated branches from dead branches. Canopy dieback was observed to some degree almost everywhere along the study reach, with only small patches of tamarisk with 0% dieback ( Figure 5). Tamarisk canopy dieback from previous events and prior years was most easily observed in May 2019 both on the ground and in the imagery, when the lack of tamarisk beetle defoliation reduced the need to discern defoliated branches from dead branches.  The ground-based ASD hyperspectral data clearly indicate distinct spectra for tamarisk compared to the other common vegetation species. Healthy, green tamarisk would be most likely to be misclassified as arrowweed, since the near-infrared region of each endmember's spectra peaks with reflectance at about 0.4. However, each of the other species has spectra that is highly discernible from tamarisk. Furthermore, the field-based hyperspectral spectra indicate unique spectral signatures for varying tamarisk conditions. Compared to the healthy, green tamarisk, the canopy dieback spectra has less variability, and lower reflectance across all spectral bands.
The May, 2019 WorldView-2 image classification had an overall accuracy of 80.0% ( Table 2). The healthy tamarisk class had producer's and user's accuracies of 87.7% and 78.1%, owing to errors of omission and commission, respectively ( Table 2). The canopy dieback class had producer's and user's accuracies of 71.1% and 81.1%, respectively ( Table 2). The defoliated tamarisk class had producer's and user's accuracies of 66.7% and 100.0%, respectively (Table 2). We estimate that the post-processing steps improved the accuracies by 10-15% and were especially important for this class due to the mixing with bare ground. The August 2019 WorldView-2 image classification had an overall accuracy of 83.1% (Table 3). The defoliated (beetle-impacted) tamarisk class had a producer's and user's accuracies of 85.4% and 91.1%, respectively ( Table 3). The healthy tamarisk class had producer's and user's accuracies of 76.5% and 65.0%, respectively (Table 3), with the lower user's accuracy possibly due to spectral mixing of green tamarisk with defoliated tamarisk later in the growing season.

Change Detection
In May, 2019, the SAM classification detected 43.6 ha of healthy tamarisk, 18.2 ha of canopy dieback, and 2.9 ha of defoliated tamarisk within the WorldView-2 study reach. The detected tamarisk amounted to 64.1% of the known total vegetated area within the study reach from Durning et al. [45]. In August 2019, 20.4 ha of healthy tamarisk and 39.7 ha of beetle-impacted tamarisk were detected, which amounted to 59.5% of the total vegetated area [45]. After applying a matching cloud mask from August to both images, May detections increased to 36.7 ha of healthy tamarisk, 14.4 ha of canopy dieback, and 2.7 ha of defoliated tamarisk, resulting in a total of 53.8 ha of tamarisk detected. August detections increased to 19.1 ha of healthy tamarisk and 35.6 ha of beetle-impacted tamarisk after applying the mask, resulting in a total of 54.7 ha of tamarisk detected. We found that 17.4 ha (47.5%) of the healthy tamarisk detected in May 2019 was classified as beetle-impacted by August, 2019 within our WorldView-2 study reach. As dieback is a slower process than defoliation, we interpret this measure as an estimate of defoliation that occurred during this period. Of the 2.7 ha of defoliated tamarisk mapped in May, 7.4% (0.2 ha) had recovered by August, and the other 92.6% (2.5 ha) remained defoliated. Of the 14.4 ha of canopy dieback mapped in May, 20.1% (2.9 ha) changed to healthy tamarisk in August, which indicates potential error within this classification scheme, reflected in the 71.1% producer's accuracy for the canopy dieback class.

WorldView-2 Image Classification of Tamarisk Distribution, Canopy Dieback, and Defoliation
The WorldView-2 image classifications had high accuracies in detecting tamarisk distribution and beetle impacts for 2019 within the 48 km river stretch of Grand Canyon National Park. Our results indicate that the high spatial resolution (2 m) of the WorldView-2 image is appropriate for mapping tamarisk distribution in the Grand Canyon, despite the inherent challenges associated with remote sensing, particularly from satellite sensors, in deep canyon topography. Having a reliable satellite sensor for such challenging study areas allows for rapid analysis of landscape changes over tens of km of riparian corridor, where biodiversity, phenologic patterns, and land cover types can be dynamic.
The overall accuracies of the May and August WorldView-2 image classifications are comparable to the overall accuracy of 74% achieved by Bedford et al. [13,14]. The reach of the river covered by the Worldview imagery is within a segment of Grand Canyon that is

WorldView-2 Image Classification of Tamarisk Distribution, Canopy Dieback, and Defoliation
The WorldView-2 image classifications had high accuracies in detecting tamarisk distribution and beetle impacts for 2019 within the 48 km river stretch of Grand Canyon National Park. Our results indicate that the high spatial resolution (2 m) of the WorldView-2 image is appropriate for mapping tamarisk distribution in the Grand Canyon, despite the inherent challenges associated with remote sensing, particularly from satellite sensors, in deep canyon topography. Having a reliable satellite sensor for such challenging study areas allows for rapid analysis of landscape changes over tens of km of riparian corridor, where biodiversity, phenologic patterns, and land cover types can be dynamic.
The overall accuracies of the May and August WorldView-2 image classifications are comparable to the overall accuracy of 74% achieved by Bedford et al. [13,14]. The reach of the river covered by the Worldview imagery is within a segment of Grand Canyon that is oriented North-South and is relatively wide. Remote sensing with Worldview might be more challenging in more sinuous and narrow segments of Grand Canyon. Our classification accuracies are also comparable to the accuracies reported by Narumalani et al. [36], who achieved 83% accuracy using SAM in Lake Meredith, Texas. However, other studies report higher accuracies than observed in our study likely due to the wider riparian corridors and larger tamarisk stands found at their study sites using airborne images and Quickbird and Landsat satellite data [33][34][35].
The lowest accuracies in the WorldView-2 maps were in the producer's accuracy of defoliated tamarisk in May and the user's accuracy for healthy tamarisk in August.
The lower user's accuracy in August could be the result of mixing of healthy tamarisk foliage with defoliated tamarisk that occurred within the same pixels, whereas the converse might be an explanation for the low producer's accuracy in May. The low producer's accuracy in May might also be the result of mixing of green understory vegetation beneath defoliated tamarisk canopies. The sample size for defoliated tamarisk in May was very small because there is typically not as much tamarisk beetle activity at that time of year compared to later in the summer; owing to the small number of observations, the accuracy measure for this class might not be reliable.
Some nuances associated with the ground-based observation and classification methods might have impacted the accuracies of the WorldView-2 classifications. The defoliation observations were made 8-9 days apart between ground-based and image-based detections allowing for defoliation levels to potentially change. Additionally, the mismatch in classification levels between the ground-based and WorldView-2 image-derived methods introduced more room for error than using discrete training data representing a single class, collected in the field via GPS and translated to training and validation pixels. Furthermore, the 52.7% agreement in tamarisk detection between the WorldView-2 map and the 2009 supervised classification map can be attributed to the fact that the 2009 data comprised USGS aerial imagery rather than a consistent satellite image source. Such differences in spatial resolution can result in varying rates of detection and different area estimates of invasive species in remote sensing efforts [46,47]. Furthermore, the 2009 aerial image used for Mahalanobis supervised classification by Bedford et al. [13,14] had 100 times finer resolution than the WorldView-2 imagery, allowing for detection of shrubs throughout the study area that were smaller and covered less than a single WorldView-2 pixel. An advantage of using imagery with resolution as high as 20 cm [13,14,48] is that different plant species can sometimes be identified by visual analysis of the imagery, allowing for more control over post-classification processing and manual cleaning of non-target species.
This study also highlights the importance of using multiple images to tease out seasonal patterns of tamarisk defoliation. It was difficult to separate the canopy dieback class from the defoliated tamarisk class in a single image due to mixing of these classes within a 2 m pixel in August, when there was greater defoliation to potentially mix with other classes, compared to the May image. Combining the canopy dieback and defoliated tamarisk classes in August eliminated the possibility of discerning canopy dieback from defoliation in a single image, leading to a classification of healthy versus "beetle-impacted" tamarisk, as was done in Sankey et al. [16] and Bedford et al. [13,14]. The changes detected between May and August 2019 were useful for measuring defoliation within the growing season, because they take advantage of the healthy-tamarisk and dieback classes in May from a time when little seasonal defoliation was observed. This approach can be useful for land managers in monitoring these phenomena remotely and across large landscapes.

Seasonal and Multi-Year Change Detection
The 2019 WorldView-2 classification maps were used to detect seasonal changes over a single growing season as well as to estimate multi-year changes in tamarisk since 2009. First, the May to August 2019 change detections indicated that 47.5% of the healthy tamarisk in May were defoliated by August, which suggests significant defoliation by the beetle had occurred within this stretch over the growing season ( Figure 7). Second, the 2009-2019 change detection indicated that 29.5% of the healthy tamarisk classified in 2009 had experienced canopy dieback by May 2019, while 66.6% remained healthy (the remainder were defoliated). Similarly, the 2013 to 2019 change detection found 64.9% of the healthy tamarisk in 2013 were still healthy by May 2019 (Figure 7), and 28.5% of the beetle-impacted tamarisk (i.e., defoliated or with dieback) in 2013 was still impacted in May 2019. The tamarisk that changed from beetle-impacted back to healthy tamarisk from 2013 to 2019 (Figure 7) were possibly either the result of tamarisk recovering new canopies after being defoliated, or growing new canopies within the dieback during periods where defoliation pressure was low in the years between maps, or possibly from new growth from plant species with similar spectral signatures like arrowweed. Previous studies have shown that the impacts of the tamarisk beetle are the greatest within the first years following the introduction of the beetle to a new environment [43]. Between 2009-2013, Bedford et al. [13,14] detected 15% of the canopy to have died back or been defoliated within the 48 km river stretch of the WorldView-2 imagery, whereas in 2019 we found 29.5% of the canopy was classified as dieback. This implies that tamarisk beetle impacts are continuing to accumulate despite beetles being present in the study reach for 7 years.
In other river systems of the SW USA, Jarchow et al. [28] and Snyder & Scott [7] have shown declines in growing season tamarisk canopy greenness measured with NDVI from moderate resolution satellites during the first 2 to 10 years after introduction of the tamarisk beetle. Jarchow et al. [28] and Snyder & Scott [7] found that the declines in greenness correlated with decreases in ET and decreases in overall riparian ecosystem health, but also determined that in many but not all of the river reaches examined, tamarisk greenness, ET, and LAI recovered to pre-beetle levels within a decade. It remains an open question as to how tamarisk will recover in the long-term from beetle impacts in our study reach of the Grand Canyon. Kasprak et al. [49] recently used habitat suitability models developed by Butterfield et al. [50] to show that much of the unvegetated riparian habitat in Grand Canyon is currently suitable for colonization and expansion of tamarisk, but the modelling did not explicitly incorporate effects of, or feedbacks between, the tamarisk beetle and tamarisk.

Ecological Implications
Our multi-year comparison maps can be used by land managers to determine where tamarisk defoliation and canopy dieback are observed due to the tamarisk beetle, and where restoration efforts thus could be focused across a large landscape. Both tamarisk defoliation and canopy dieback have important implications for wildlife. For example, leaf loss of tamarisk reduces habitat for songbirds including the endangered Southwestern Willow Flycatcher [31,32], and other riparian nesting birds [51]. Herpetofauna that use the understory as habitat can also be negatively impacted by defoliation when alternative habitat is not available [30]. Groups like the USDA Animal and Plant Health Inspection Service (APHIS) Technical Advisory Group for Biological Control Agents of Weeds, or RiversEdge West, who communicate science to the public, provide guidance to researchers, or review petitions for biocontrol of weeds in the US, could potentially use these data to help inform future biocontrol efforts.
Changes in tamarisk and landscape-scale dynamics as illustrated in Figure 7 also directly influence plant cover and species richness [52][53][54]. Tamarisk abundance and cover decreases following tamarisk beetle defoliation and canopy dieback and can result in regeneration of other plant species [55]. Some increase in native plant cover following tamarisk removal or dieback is attributed to replacement by arrowweed in other river systems [53,54], which was also documented in our field observations along the Colorado River in Grand Canyon. González et al. [54] found that tamarisk in the southwestern U.S. is often replaced by secondary invasions of noxious weeds, including Russian knapweed (Acroptilon repens), cheatgrass (Bromus tectorum), summer cypress (Bassia scoparia), and prickly Russian thistle (Salsola tragus), though the success of each weed depends on if control measures are taken. Many of those same weeds occur in Grand Canyon, and thus as tamarisk cover is lost through canopy dieback in Grand Canyon, the replacement of tamarisk with native shrub species that perform similar ecological functions might be important for maintaining ecosystem health.

Implications for Future Remote Sensing Research
Bedford et al. [13,14] demonstrated the utility of 20-cm pixel resolution four-band imagery acquired from a multi-spectral sensor onboard a fixed-wing aircraft to accurately classify healthy tamarisk and beetle-impacted tamarisk, and detect changes over time in those classes along the Colorado River in Grand Canyon. Other invasive species studies across North America have used hyperspectral data such as AVIRIS imagery with 20 m resolution and 224 bands and HyMap imagery with 3.5 m resolution and 126 bands [56][57][58][59][60][61]. However, hyperspectral imagery can be expensive to acquire and tends to cover small spatial extents. Even multispectral, airborne images are not acquired frequently in Grand Canyon (e.g., annually or intra-annually) because they are also expensive and logistically challenging. In this study, we show that 2-m resolution multispectral satellite imagery, which is available at much greater temporal frequency, can be used to detect intra-annual seasonal changes in tamarisk owing to the tamarisk beetle by detecting the inherent spectral differences between healthy tamarisk, defoliated tamarisk, and tamarisk canopy dieback. In other arid and semi-arid environments around the world without the very steep topography, coarser resolution, multi-temporal satellite data such as Landsat have been used to similarly leverage the seasonal variability in invasive species detection [62].
A potential improvement on our methods would be to incorporate plant height estimates from lidar or photogrammetry along with the spectra in the classification to improve accuracies. Other studies have demonstrated that lidar-derived topographic data and unmanned aerial vehicle (UAV)-based photogrammetry-derived vegetation height data can be integrated with high resolution multispectral images to better detect invasive species distribution and estimate expansion rates [46,63,64]. We did not have lidar data for plant height or topography for this reach of the river, though Sankey et al. [16] were able to use that approach with lidar and multispectral imagery to detect tamarisk beetle impacts along the Colorado River in Glen Canyon, Arizona.

Conclusions
Using high spatial resolution WorldView-2 imagery from 18 May 2019 and 2 August 2019, we classified tamarisk distribution in the topographically complex study area of Grand Canyon National Park and demonstrate tamarisk beetle impacts in the park since 2009. A strength of our approach was that the 2019 WorldView-2 image classifications could be compared with previously published maps derived from aerial multispectral imagery acquired in 2009 and 2013. Taken together, our long-term results show that tamarisk beetle impacts are continuing to accumulate even though land managers have noted the presence of the beetles in the park since 2009 and specifically in the study reach of the river for 7 years since 2012. Our methods and maps [65] provide a valuable tool for land managers to track defoliation and canopy dieback of tamarisk from the tamarisk beetle over the region. These data can help determine what restoration efforts are needed in different areas. In areas where tamarisk is not being defoliated or killed by the beetle, thinning and other treatments could be preferred, while areas with high levels of mortality could require the removal of dead shrubs. This study is the first known effort of mapping decadal-scale tamarisk canopy dieback and seasonal-scale defoliation in the Grand Canyon with satellite remote sensing, setting the stage for future studies to conduct spatial analysis of defoliation, dieback, and eventual multi-year tamarisk mortality in Grand Canyon National Park. Further research and monitoring will provide useful insights regarding how the impacts of the beetle on tamarisk will affect this iconic riparian corridor.