Remotely Sensed Changes in Vegetation Cover Distribution and Groundwater along the Lower Gila River

Introduced as a soil erosion deterrent, salt cedar has become a menace along riverbeds in the desert southwest. Salt cedar replaces native species, permanently altering the structure, composition, function, and natural processes of the landscape. Remote sensing technologies have the potential to monitor the level of invasion and its impacts on ecosystem services. In this research, we developed a species map by segmenting and classifying various species along a stretch of the Lower Gila River. We calculated metrics from high-resolution multispectral imagery and light detection and ranging (LiDAR) data to identify salt cedar, mesquite, and creosote. Analysts derived training and validation information from drone-acquired orthophotos to achieve an overall accuracy of 94%. It is clear from the results that salt cedar completely dominates the study area with small numbers of mesquite and creosote present. We also show that vegetation has declined in the study area over the last 25 years. We discuss how water usage may be influencing the plant health and biodiversity in the region. An examination of ground well, stream gauge, and Gravity Recovery and Climate Experiment (GRACE) groundwater storage data indicates a decline in water levels near the study area over the last 25 years.


Introduction
Salt cedar (Tamarix chinensis) is invading riparian corridors in the semiarid southwest at an alarming rate, affecting ecosystem processes and services such as water availability, carbon sequestration, and nutrient cycling. Drought tolerant, flood tolerant, fire adapted, and with a tendency to grow in dense stands, salt cedars replace other arid and semiarid species [1][2][3]. The aggressive spread of salt cedar is decreasing resources for native cottonwood and willow species, ideal habitats for native wildlife [4,5]. However, some species, such as the southwestern willow flycatcher, have adapted to the salt cedar influx and the changes aiding its expansion [6,7].
Introduced as ornamentals and windbreaks in the early 1800s, salt cedar is an ideal plant for erosion control, both inland and coastally, due to its high salt tolerance [8]. Introduced, salt cedar spreads aggressively along the flood plains of rivers in the southwestern United States [8][9][10][11][12]. Coverage of salt cedar increased from about 4000 ha in the 1920s to more than 600,000 ha in 1998 [1]. It is unclear how much area salt cedar currently occupies in the southwestern United States, but it is clear that the species is present within most areas suitable to its growth [13]. Salt cedar continues to expand along rivers in the western United States at a rate of around 20-25 km per year [13,14]. Arizona is ideal for agriculture due to 300 plus days of sunshine per year. However, due to inadequate rainfall for agricultural production, Arizona's farmers rely on irrigation via river water or groundwater pumping [15]. Between 2001 and 2005, farmers used 775,858,920 m 3 of water per year to irrigate crops in the Lower Gila Basin [16]. Farmers grow various vegetables, alfalfa, melons, and wheat along the Gila River with 60 percent of the water diverted from the Colorado River. The future of irrigation via river water in the region is uncertain as water levels continue to decline due to lengthy drought conditions [16]. Diversion of water for agriculture creates a prime environment for salt cedar. However, uncertainty with regards to water leads to agricultural abandonment in many arid environments [17,18], allowing salt cedar to encroach further.
To improve remediation efforts of ecosystem services in tamarisk-invaded (salt cedar) riparian ecosystems, land managers need an inventory of current species distribution as well as historical information on the trend of salt cedar expansion. Such information will inform the resources needed for restoration and future control efforts along with the impact of such practices [19,20]. Historical remote sensing imagery and modern analysis methods are ideal tools for inventory and monitoring of salt cedar.
Arid environments typically have sparse vegetation cover and often consist of plants that evolved small leaves in order to decrease evapotranspiration. These physical properties challenge traditional satellite-based vegetation indices [21]. A better tool to map salt cedar (and other arid land plant species) is high-spatial-resolution aerial imagery [22]. With this imagery, the separation between plants actually aides our ability to identify individual plants (i.e., segmentation) [23]. We can further enhance species identification by fusing complementary remotely sensed data [24]. For example, light detection and ranging (LiDAR) can capture individual plant height and structure information [25][26][27][28][29].
Focusing on a 42 km stretch of the Gila River in southwest Arizona, our research had two main questions. First, which plant species are currently present in the study area? We developed a method combining high-resolution multispectral aerial imagery, aerial LiDAR, and small drone imagery to produce a contemporary species map of the area. We examine the accuracy of identifying different species and discuss how species size and spectral characteristics influence the results. Second, what is the trend of ecosystem health in the study area using vegetation as a proxy? We sought to characterize the changing extent of vegetation cover by examining a time series of aerial imagery from 1996 to present day. We examine how anthropogenic water use and a changing climate are impacting available water for vegetation. We compare the vegetation maps to ground well, precipitation, stream flow, and water equivalent thickness data to investigate what may be exacerbating salt cedar expansion and overcrowding. We also discuss trends in agricultural patterns in relation to groundwater pumping.

Study Area
Our study area encompassed a nearly 42 km stretch of the Lower Gila River about 100 km northeast of Yuma, Arizona between Sentinel and Dateland, Arizona ( Figure 1). The annual average temperature in the region is 23 degrees Celsius with annual average highs of 31 degrees Celsius and annual average lows of 14.5 degrees Celsius. Between 1994 and 2019, the average low has increased by about 2.25 degrees Celsius. The range between the minimum and maximum temperatures has decreased by about three degrees Celsius.
Average annual precipitation is just under 107 mm following a bimodal pattern with the majority of rainfall in the winter and summer months. Year-to-year precipitation is highly variable with a low of 15 mm in 1996 and a high of 166 mm in 2018. There has not been any clear decline in precipitation levels in the region over the last 25 years. Streamflow is minimal in the study area. Between 1994 and 2019, the U.S. Geological Survey (USGS) recorded less than 0.03 m 3 per second of annual discharge 21 times (https://waterdata.usgs.gov/nwis). encompassed about 50% of the study area [10].
According to the Bureau of Land Management (BLM), the current ecosystem along the Lower Gila River is composed of native species including mesquite (Prosopis velutina), creosote bush (Larrea tridentate), brittlebush (Encelia farinose), and arrowweed (Pluchea sericea). Cottonwood (Populus fremontii) and desert willow (Chilopsis linearis) are present near perennial waters, but those are rare or disappearing. The invasive salt cedar (Tamarix chinensis) encompasses most of the study area. Figure 1. Location of Study Area. Our study area, indicated in black on the inset map, is composed of a nearly 50 km stretch of the Gila River located between Sentinel, AZ (east) and Dateland, AZ (west). This 60 cm natural color National Agriculture Imagery Program (NAIP) image from 2019 shows private agricultural land use both north and south of the Gila River (blue line). The area outlined in yellow was converted from agriculture into a solar energy facility after 2010. Lands managed by the Bureau of Land Management (BLM) are outlined in red, which also delineates the extent of the light detection and ranging (LiDAR) collection used in this paper. Cities of interest are shown as black stars. Figure 1. Location of Study Area. Our study area, indicated in black on the inset map, is composed of a nearly 50 km stretch of the Gila River located between Sentinel, AZ (east) and Dateland, AZ (west). This 60 cm natural color National Agriculture Imagery Program (NAIP) image from 2019 shows private agricultural land use both north and south of the Gila River (blue line). The area outlined in yellow was converted from agriculture into a solar energy facility after 2010. Lands managed by the Bureau of Land Management (BLM) are outlined in red, which also delineates the extent of the light detection and ranging (LiDAR) collection used in this paper. Cities of interest are shown as black stars.
Elevation at the northeastern side of the study area is 184 m above mean sea level sloping down to 110 m at the southwest side of the study area. The majority of the region is composed of entisols because the area is essentially all flood plain or riverbed (https://websoilsurvey.sc.egov.usda.gov/). The Lower Gila River is a riparian system altered over the years through water diversion for agricultural purposes, electricity, and drought, creating a prime ecosystem for salt cedar. In the 1700s, fishing was possible on this stretch of the river and in the mid-1800s; cottonwoods grew to 9 m tall [10]. A 1972 survey identified six communities: cattails (Typha), salt cedar and arrowweed (Tamarix-Pluchea), mesquite (Prosopis), seepweed and pickleweed (Suaeda-Allenrolfea), saltbush (Atriplex), and creosote bush and mesquite (Larrea-Prosopis) [10]. The 1972 survey also stated that salt cedar encompassed about 50% of the study area [10].
According to the Bureau of Land Management (BLM), the current ecosystem along the Lower Gila River is composed of native species including mesquite (Prosopis velutina), creosote bush (Larrea tridentate), brittlebush (Encelia farinose), and arrowweed (Pluchea sericea). Cottonwood (Populus fremontii) and desert willow (Chilopsis linearis) are present near perennial waters, but those are rare or disappearing. The invasive salt cedar (Tamarix chinensis) encompasses most of the study area. The BLM acquisitioned the collection of digital orthoimagery over the study area, around 59 km 2 . On 22 May 2016, the Atlantic Group, LLC captured four-band multispectral aerial orthoimages. The nominal spatial resolution of the BLM data was 15 cm.
We created NDVI images from the red and near-infrared bands of all datasets.

Light Detection and Ranging (LiDAR)
On March 25, 2016, the Atlantic Group, LLC collected LiDAR data over the study area. The LiDAR data have a nominal pulse density of 16 pulses per m 2 with a nominal pulse spacing of 0.2 m.
The Atlantic Group, LLC provided a LiDAR point cloud containing all returns along with one representing only bare earth data. Using FUSION [30], we created a Digital Elevation Model (DEM) raster with a spatial resolution of 50 cm based on the mean value of ground returns in each 50 cm cell. This was the finest resolution we could achieve without creating a surface with missing data due to a lack of returns in a cell.
We subtracted the bare earth DEM from the first return point cloud to create an Object Height Model (OHM). The OHM provided the normalized height of objects, including trees, above the ground. This allows the user to compare trees at different elevations as the subtraction normalizes for elevation.
We used the DensityMetrics tool in FUSION to compute point cloud densities at various height ranges. The DensityMetrics tool outputs a raster for each height interval dictated by the user wherein each pixel contains the proportion of LiDAR returns within that interval [30]. We calculated the densities of the point cloud within 50 cm cells at 1 m intervals, 0-1 m, 1-2 m, 2-3 m, 3-4 m, 4-5 m, and 5-6 m, for a total of six density measurements. There were no trees in the study area with a height greater than 6 m.
We used the IntensityImage tool in Fusion to convert the point cloud intensities to grid form. LIDAR intensity raster records the return strength of each LiDAR pulse in relation to the strength of the LiDAR pulse when the sensor emits it [30]. We calculated the average intensity within 50 cm cells.
We used the GridMetrics tool in Fusion to compute kurtosis and skewness from the original LiDAR point cloud and output the metrics as rasters with 50 cm cells. Both kurtosis and skewness characterize the distribution of LiDAR returns for each canopy in the study area. Skewness summarizes the asymmetry of points within each canopy, while kurtosis indicates how "peaky" the point distribution is in each canopy [30]. We selected kurtosis and skewness because of the compactness of mesquites compared to the more ambiguously shaped salt cedar.

Drone and Field Data
We collected drone imagery and ground-based photographs on 30 September 2016, 10 April 2017, 9 February 2018, and 8 May 2018 to identify vegetation species to train and validate our aerial imagery classification products. We collected the data in different seasons hoping to capture phenological differences between species. Our ground-based training methods started out as technicians walking through the study area on foot and geolocating class training data (tamarix, mesquite, other) using Land 2020, 9, 326 5 of 18 a GPS-enabled, hand-held camera or an iPad loaded with Galileo Offline Maps. We subsequently discovered that we could collect far more training data over a larger area with UAV imagery.
We flew a Phantom 3 Quadcopter and a Phantom 4 Pro Quadcopter 60 m above ground level within portions of the study area. The quadcopters carried a 12-and 20-megapixel camera, respectively, that captured standard red, green, blue color photographs. We surveyed eight ground control targets per acquisition with a hand-held Trimble GPS unit in order to produce orthomosaics in the same coordinate system as the 15 cm manned imagery. In total, we surveyed just under 20 hectares of the study area using the UAVs.
We stitched together~1000 photos from each of the four drone flights using Agisoft Photoscan to create four high-resolution orthomosaics with ground sampling distances (pixel size) of 1.7 cm. The imagery was at a high enough resolution to visually identify and distinguish vegetation species including salt cedar, mesquite, and creosote.

Groundwater, Streamflow, Precipitation, and Water Equivalent Thickness Measurements
We acquired groundwater data from the Arizona Department of Water Resources (ADWR-https: //waterdata.usgs.gov/nwis). We identified 15 different wells with measurements from 1994 through 2019. We calculated the percent change in depth to water for the 25-year period.
We also acquired streamflow data for the Gila River from a set of USGS (https://waterdata.usgs. gov/nwis) gauges located near Dateland (09520280) and Painted Rock Dam (09519800). The records include data from 1994 through 2020.
We used Parameter-Elevation Regressions on Independent Slopes Model (PRISM-https://prism. oregonstate.edu/) data to assess rainfall in the study area. PRISM data has a spatial resolution of 4 km. We extracted total monthly precipitation for the study area from 1994 through 2019. We calculated quarterly and yearly summations. We examined trends in precipitation and correlations between streamflow and precipitation.
We also examined water equivalent thickness measured by Gravity Recovery and Climate Experiment (GRACE) sensors [31][32][33]. GRACE measures the distance between twin satellites circling the earth to make gravitational field measurements. Collection of GRACE data began in 2002 and continues through to present day. Water equivalent thickness acts as a proxy for underground water storage provided as a three-degree equal-area grid.

Classification of Vegetation in 2016 Using Acquired LiDAR and Multispectral Data
We used the C50 classification and regression tree (CART) R package to perform our classifications in this paper [34,35]. The C50 CART algorithm is a supervised machine learning technique that uses training data to create classification trees used to assign each pixel in a raster to a class. The tree outputs allow the user to assess confusion between classes and make adjustments to the training data to help distinguish classes.
The spatial resolution, 15 cm, of the 2016 BLM image data was fine enough to produce a map of vegetation species in the study area. Our workflow to create a classified species map began with a C50 pixel-based classification into three general classes (detritus, vegetation, bare ground) (Figure 2A). We loaded the drone and field data into ArcGIS to select pixels representing detritus, vegetation, and bare ground on the BLM orthoimagery. Using ArcGIS, we created separate point shapefiles for each class containing 100 points representing 100 pixels. We performed the classification using the four-band orthoimagery and derived NDVI image. The classifier only needed one trial and winnowed three of the five inputs. It separated the three classes using information from the red band and NDVI layer. The classifier achieved an overall accuracy of 100 percent. The bright soil in the region made it easy to distinguish living vegetation and detritus.

Vegetation Segment Metrics
We used zonal statistics in ArcMap to calculate several statistical measures for each vegetation segment ( Figure 2C). We calculated the mean, maximum, minimum, range, and standard deviation of each band of the orthoimagery and the NDVI layer. We also calculated the mean, maximum, minimum, range, and standard deviation of the LiDAR height, skewness, kurtosis, and intensity information for each segment. Using the LiDAR density slices, we calculated the percentage of LiDAR returns in each segment for each of the six height levels. Finally, we counted the total number of 15 cm pixels within each vegetation segment. We used the various metrics in hopes of capturing the variance within each vegetation segment that would allow us to separate species (Table 1). We then merged the detritus and bare-ground classes to produce a vegetation and nonvegetation class. We used the region grouping tool, on pixels classified as vegetation, to combine adjacent pixels into segments (sometimes called objects). Each segment represents essentially one plant (except when canopies of multiple plants touched) ( Figure 2B).

Vegetation Segment Metrics
We used zonal statistics in ArcMap to calculate several statistical measures for each vegetation segment ( Figure 2C). We calculated the mean, maximum, minimum, range, and standard deviation of each band of the orthoimagery and the NDVI layer. We also calculated the mean, maximum, minimum, range, and standard deviation of the LiDAR height, skewness, kurtosis, and intensity information for each segment. Using the LiDAR density slices, we calculated the percentage of LiDAR returns in each segment for each of the six height levels. Finally, we counted the total number of 15 cm pixels within each vegetation segment. We used the various metrics in hopes of capturing the variance within each vegetation segment that would allow us to separate species (Table 1). Using imagery collected via drone and GPS camera, we selected 258 segments representing individual salt cedars, 6 individual mesquites, and 24 individual creosotes. The amount of segments selected to represent each class varied because of the abundance of each species present in the study area. We could not accurately identify arrowweed and brittlebush on the 15 cm resolution orthophotos. However, it might be possible to pinpoint them on the higher resolution unmanned aerial vehicle (UAV) image data.
We partitioned the data in R, taking 50 percent of each class' samples for training and setting aside 50 percent of each class' samples for validation. We ran the C50 classifier on the various metrics in Table 1 with a maximum of 10 iterations. We also enabled boosting and winnowing of attributes.

Classification of Vegetation in 1996, 2007, and 2019 Using NAPP and NAIP Data
The spatial resolutions of the 1996, 2007, and 2019 data were not detailed enough to separate vegetation species from one another. However, we could identify areas of live vegetation, detritus, and bare ground.
Using the C50 package in R, we performed supervised classifications of pixels representing vegetation, detritus, and bare ground. We used the same protocol to train the C50 models for the 1996, 2007, and 2019 imagery sets. For training data, we selected 100 pixels, unique to each year, to represent the three classes. For validation purposes, we selected an additional 100 pixels per class, also unique to each year. We selected training samples using ESRI ArcMap 10.8 with Google Earth imagery as a reference. Class and image values (blue, green, red, NIR, NDVI) from each pixel were stored in point shapefiles.
C50 assesses the modeled output using both the training and validation pixels. Each classification achieved an overall accuracy of 100% when assessed with both sets of pixels. We manually delineated areas of active agricultural land in each of the three images.

Change Assessment
We created a tessellation, composed of 250 m squares (62,500 m 2 ), using ArcGIS. We then used the Zonal Statistics tool to calculate the total number of vegetation pixels in each 250 m square from the 1996, 2007, and 2019 data. Using the Raster Calculator tool in ArcGIS, we calculated the percentage of vegetation in each square for each year. We subtracted the three tessellations to identify areas of increased vegetation and areas of decreased vegetation between each year and over the entire Land 2020, 9, 326 8 of 18 study period. We used the 250 m tessellation in hopes of identifying hotspots of vegetation loss and expansion.

Buffer Analysis
We obtained a shapefile representing the Gila River from the USGS National Hydrography Dataset. We created four concentric rings around the shapefile at 0-100 m, 100 m-500 m, 500 m-1000 m, and 1000 m-2000 m to examine vegetation in relation to the central river channel.

Living Vegetation Segmentation Classification
The classification of the high-resolution aerial imagery and LiDAR from 2016 achieved an overall accuracy of 94 percent ( Table 2). The C50 algorithm needed 13 out of 52 variables to separate the four different vegetation classes. The mining procedure used minimum, maximum, mean, and range of the blue band, mean of the green band, maximum, range, and standard deviation of the near infrared band, maximum and mean of the NDVI layer, density of LiDAR returns between 0 and 1 m, maximum of the OHM, and minimum skewness of the LiDAR returns. The product exhibited confusion within all three classes. The classification process confused some of the larger creosote bushes for salt cedars. The algorithm also confused some of the small salt cedars for larger creosote bushes. The C50 algorithm correctly identified two of the three mesquite segments. The classification process misidentified only four of the segments representing salt cedars. It confused two with creosote and two with mesquite. The algorithm, in general, did an excellent job with the salt cedar class (    Figure A3). This equates to around a 6 percent net decrease in vegetation over the entire study area during the last 23 years. There are certain portions of the study area showing pockets of increased vegetation, but most of the area has seen a decrease in vegetation ( Figure 4). It seems most gains in vegetation are evident along stream channels. We used the buffer analysis to assess this phenomenon.    Figure A3). This equates to around a 6 percent net decrease in vegetation over the entire study area during the last 23 years. There are certain portions of the study area showing pockets of increased vegetation, but most of the area has seen a decrease in vegetation ( Figure 4). It seems most gains in vegetation are evident along stream channels. We used the buffer analysis to assess this phenomenon.

Vegetation Change and River Buffer Analysis
Based on the buffer analysis, the amount of vegetation was always greatest near the river and declined moving outward ( Figure 5). The 0-100 m buffer also registered the greatest amount of vegetation loss, declining from 37% in 1996 to 15% in 2019. All four buffer levels exhibited the greatest amount of vegetation coverage in 1996 and the least amount of vegetation coverage in 2019.
Based on the results of the classification processes, vegetation covered about 10 percent of the study area in 1996 ( Figure A1), about 7 percent in 2007 ( Figure A2), and around 4 percent in 2019 ( Figure A3). This equates to around a 6 percent net decrease in vegetation over the entire study area during the last 23 years. There are certain portions of the study area showing pockets of increased vegetation, but most of the area has seen a decrease in vegetation (Figure 4). It seems most gains in vegetation are evident along stream channels. We used the buffer analysis to assess this phenomenon.

Analysis of Water Data
Fifteen of the twenty-one wells around the study area exhibit declining water levels ( Figure 6). Based on the monthly water equivalent thickness measured by the GRACE, the amount of water in the region is declining (Figure 7).

Analysis of Water Data
Fifteen of the twenty-one wells around the study area exhibit declining water levels ( Figure 6). Based on the monthly water equivalent thickness measured by the GRACE, the amount of water in the region is declining (Figure 7).

Analysis of Water Data
Fifteen of the twenty-one wells around the study area exhibit declining water levels ( Figure 6). Based on the monthly water equivalent thickness measured by the GRACE, the amount of water in the region is declining (Figure 7). Figure 6. Location of wells and stream gauges in relation to the study area. There are 21 indexed wells, shown as circles, located in and around the study area. Larger circles represent wells with lower depth to water and smaller circles represent greater depth to water. Wells shown in red indicate wells where the depth to water has increased since 1992, while wells shown in blue indicate wells where the depth to water has decreased since 1992. The well displayed in black did not register a change in depth. USGS stream gauge locations are shown as green triangles. The Dateland gauge (09520280) is on the western side of the map and the Painted Dam gauge (09519800) is on the eastern side of the map. The study area shown in Figure 1 is represented by the black box.  Declining levels of water make it difficult for vegetation to remain established in the region. There is a clear decline in the amount of vegetation near the river channel and areas further from the river that already demonstrated low levels of vegetation show decline as well ( Figure 5).

Limitations and Improvements
For the 2016 species classification, we had a difficult time separating mesquite from salt cedar likely due to the limited point density of the LiDAR data, reducing our ability to characterize tree shape and height differences. The LiDAR used in this paper did not have a point density conducive to vegetation mapping in arid regions. We would need a LiDAR collection with a point density closer to 22 points per m 2 in order to create LiDAR raster products with a pixel size near the 30 cm resolution of the multispectral data [36]. Ideally, to best separate tree species, multiple LiDAR returns, above, within, and below the canopy, would improve characterization of canopy structure.
UAV-based photogrammetric point clouds may aid vegetation classification in ways that LiDAR cannot. Specifically, those point clouds are typically much denser than LiDAR and may be able to better capture structural differences between vegetation species. Future work may include this technology with the recognition that UAVs can only cover small areas relative to the study area of this project. It may be relevant to use a UAV approach for prime conservation treatment and restoration area candidates.
Other work demonstrates the ability to separate tree species using hyperspectral data [37][38][39]. However, hyperspectral data are not readily available, are expensive to acquire, and require large storage capacity and powerful computing facilities. The user must employ advanced analysis techniques to reduce the data cube and extract the relevant information for species separation.

Vegetation and Water
It is clear from the species classification that salt cedar dominates our study area along the Lower Gila River. Based on a 1972 survey [10], researchers estimated that more than 50 percent of the plants in the floodplain were salt cedar. From the maps produced in our study, the amount of salt cedar in the study area has increased. Based on our 2016 classification, salt cedar now encompasses close to 90 percent of the living vegetation in our study region.
The 1972 survey also identified communities of Typha and Suaeda-Allenrolfea. Typha, commonly known as cattail, is reliant on standing or permanently running water. The stream gauge data for the study area indicates that stream flow is intermittent and minimal. Based on examination of the aerial photographs and our drone data, no standing water is present in the region. We did not encounter either of these communities on our many trips to the study area or identify them on the aerial photography and drone data. They may no longer be present in the study area, however, it would be difficult to identify either of these communities using our high-spatial-resolution aerial photographs or UAV data.
The change analysis showed that a majority of the study area lost vegetation cover between 1996 and 2019. We see the greatest loss in the 0-100 m buffer zone ( Figure 5). We posed that agricultural activities were drawing water from the region. The amount of agriculture in the study area was low but doubled from 1 to 2 percent. Active agriculture made up one percent of the study area in 1996, increasing to two percent in 2007 and holding constant in 2019. However, the location of agricultural endeavors shifted. A solar facility replaced one of the larger sets of fields between 2007 and 2019 ( Figure 1). A series of central pivot fields appeared between 1996 and 2007 and were still present in 2019.
Pumping of groundwater in the region may explain the vegetation loss in areas inside and outside of the stream channel. Groundwater well records indicate that the depth to water has increased throughout the area since 1994. We examined 21 wells around the study area and 15 of them show Land 2020, 9,326 13 of 18 that the water table has declined between 1994 and 2019. Five of the wells surveyed are labeled for agricultural use and only one demonstrates an easier access to water ( Figure 6).
Streamflow occasionally occurs in the region with extreme events happening intermittently. A large amount of water flowed along the Gila River in 1993/1994, 1995, 2005, and 2010 according to a USGS stream gauge (09520280) located on the river just north of Dateland. These flows only occur when water is released from the Painted Rock Dam upstream as the measurements from the surface water gauge (09519800) just below the dam have a 97 percent correlation with the measurements taken at the Dateland gauge. The Painted Rock Dam was built for flood control purposes with releases occurring due to storm events upstream [40].
Invasive salt cedar chokes stream flow further up the Gila River near Buckeye, Arizona. Community members are concerned about habitat degradation and increased wildfire risks [41]. A study conducted along the Pecos River in Texas discovered that removal of salt cedar along the river did not influence the stream flow. The authors state this is likely related to low sapwood area and the limited size of the salt cedar stands [42]. The situation on the Pecos River seems slightly different from the one along the Gila River as the salt cedar grows along the riverbank compared to in the main channel. Salt cedar is also present in the northern U.S. along the Musselshell River and Fort Peck Reservoir in Montana where scientists examined dispersal mechanisms [43].
There are no naturally occurring predators to salt cedar in North America to help limit the spread of the plant [43]. Managers have pursued various control methods over the years including cutting, burning, and removal to no avail [8,44]. Aerial spraying is a treatment option that seems to work, but requires years of reapplication costing millions of dollars [44,45]. In 2001, the United States Department of Agriculture (USDA) released the tamarisk beetle in California, Colorado, Utah, and Texas in an attempt to reduce the spread of salt cedar [44,[46][47][48]. The beetles defoliate the salt cedar, disabling its ability to create seeds and weakening it to other elements. However, even successful management does not always mean that native species will return to dominance [11].
When the water does flow along the Gila River, it acts like a forest fire, clearing out dead detritus and allowing new vegetation to sprout. The extended length of the flows allows the water to infiltrate the soil and hydrate vegetation along the river.

Conclusions
We conducted this research to identify species along the Lower Gila River and determined salt cedar to be the dominant woody plant species. It is clear that the Lower Gila River is a prime environment for invasive salt cedar. The diversity of habitats documented in the 1970s has decreased. Species dependent on permanent water sources are no longer present in the region. Mesquite is present, but is having a hard time competing with salt cedar capable of drawing water from the surface and subterrain aquifers. Creosote, not dependent on water from the river, continues to be moderately present in the region.
We also examined changes in total vegetation cover over a 23-year period (1996 to 2019). The results show a decline of about 6 percent. The high accuracy of vegetation classifications indicates that vegetation loss is likely. While in the field, we observed numerous salt cedar skeletons, shown as detritus on the 2016 classification, to corroborate this conclusion. Locations of die-off may be an indication of water shortages in the region. If salt cedar can no longer survive along the Lower Gila River, few other plant species can.
The ability to measure vegetation conditions and species composition provides information about past and current ecosystem health. Once managers begin treatments along the Lower Gila River, we can use the techniques in this research to monitor and assess their impacts on the landscape. Creating a time series of maps will provide valuable information about the impact of restoration practices on the landscape.
A combination of environmental and human influences has led to a decline in biodiversity and a decline in healthy vegetation in general. The combination of damming upriver, groundwater withdrawal Land 2020, 9, 326 14 of 18 for agricultural purposes, and climate change has increased the chance of ecosystem collapse. Based on projections, the study area is expected to see an increase in overall temperatures, creating an even more arid environment [49]. If these trends continue, the amount of vegetation in the area will continue to decline and the amount of bare soil will continue to increase. This could lead to dust problems seen in other parts of the state and the world [50,51]. In the worst-case scenario, complete ecosystem transformation is possible.