Combining Optical and Radar Satellite Imagery to Investigate the Surface Properties and Evolution of the Lordsburg Playa, New Mexico, USA

Abstract

Driven by erodible soil, hydrological stresses, land use/land cover (LULC) changes, and meteorological parameters, windblown dust events initiated from Lordsburg Playa, New Mexico, United States, threaten public safety and health through low visibility and exposure to dust emissions. Combining optical and radar satellite imagery products can provide invaluable benefits in characterizing surface properties of desert playas—a potent landform for wind erosion. The optical images provide a long-term data record, while radar images can observe land surface irrespective of clouds, darkness, and precipitation. As a home for optical and radar imagery, powerful algorithms, cloud computing infrastructure, and application programming interface applications, Google Earth Engine (GEE) is an invaluable resource facilitating acquisition, processing, and analysis. In this study, the fractional abundance of soil, vegetation, and water endmembers were determined from pixel mixtures using the linear spectral unmixing model in GEE for Lordsburg Playa. For this approach, Landsat 5 and 8 images at 30 m spatial resolution and Sentinel-2 images at 10–20 m spatial resolution were used. Employing the Interferometric Synthetic Aperture Radar (InSAR) techniques, the playa’s land surface changes and possible sinks for sediment loading from the surrounding catchment area were identified. In this data recipe, a pair of Sentinel-1 images bracketing a monsoon day with high rainfall and a pair of images representing spring (dry, windy) and monsoon seasons were used. The combination of optical and radar images significantly improved the effort to identify long-term changes in the playa and locations within the playa susceptible to hydrological stresses and LULC changes. The linear spectral unmixing algorithm addressed the limitation of Landsat and Sentinel-2 images related to their moderate spatial resolutions. The application of GEE facilitated the study by minimizing the time required for acquisition, processing, and analysis of images, and storage required for the big satellite data.

1. Introduction

Playas are one of the major sources of atmospheric dust [1,2,3,4,5], attributed to the availability of thick deposits of erodible clastic and chemical sediments [1], often inherited from wetter conditions during the Pleistocene and Holocene epochs [6], and unvegetated land cover type [4,7,8]. They are generally dry lake beds in internal draining basins, common in arid and semi-arid environments where evaporation exceeds precipitation and inflow [9,10]. Globally, ephemeral water bodies account for around 31% of dust emissions [11]. For example, Bodele Depression, part of a dry lake bed at the southern edge of the Sahara Desert, is the world’s greatest source of mineral dust, contributing between 6 and 18% to the global atmospheric dust composition [5,12]. Likewise, playas in southern Africa were responsible for more than half of the all the dust plumes detected over the region [13]. At the northern edge of the Mojave Desert, the dried bed of Owens Lake, California, United States, was described as the most potent source of mineral dust in North America [14,15]. At the eastern Great Basin of Utah, United States, barren playa surfaces were identified as the major dust sources, recording the largest number of observed dust plumes ( 60% of the plumes) [3]. Ephemeral lakes in the Chihuahuan Desert of Southwest North America, despite a very small areal coverage (only 4%), were the sources of 48% of the observed dust plumes in the desert [16].

The land surface’s response to the wind erosion is determined by a combination of physical processes including soil moisture, surface roughness, vegetation cover, degree of surface crusting, and climate [17,18,19,20,21]. Soil moisture increases threshold friction velocity for dust emission by increasing cohesion between soil particles [22,23]. Density of vegetation cover, rocks, and microtopography affects threshold friction velocity by modulating surface roughness [20,24,25]. Diminishing of grass cover below 9% can increase dust flux dramatically [26]. Bolstering this observation, climate-induced decrease in perennial vegetation cover coupled with the disturbance of biological soil crusts triggers an exponential increase of aeolian sediment flux [27]. The setting of land cover in terms of spatial pattern, distribution, density, and vegetation height also play a significant role in the land surface’s response to erosion and deposition of soil particles [21,28,29,30].

Hydrology linked inundation of playas is a crucial factor in controlling the amount of sediment loading from watersheds and their deposition on playa surfaces, leading to regional and local changes of dust storm hotspots [31,32]. Changes in the frequency and intensity of runoff and spatial extent of inundated area of playas may lead to a significant fluctuation of sediment availability on seasonal and inter-annual bases [33,34]. The movement of the edge of a playa’s water body also dictates the variation in the spatial pattern of sediment deposition across the playa surface and the area susceptible to aeolian erosion [34,35]. Heavy precipitation may eventually trigger an increase of dust events through the development of soft surfaces with thin and newly formed crusts, production of loose aggregates of evaporite minerals, and rise of water table’s hydraulic head [19]. These processes create erodible surfaces by increasing rates of vapor discharge, formation of loose evaporite minerals, and temporary softening of playa surfaces. During sediment runoff onto playas from surrounding uplands, coarse particles deposit at the edge of the playa and fine particles accumulate at the inner part of the playa, primarily influenced by their settling velocity [36]. When a playa surface rich with loose sediments is coupled with strong surface winds, dust events develop. Dust particles are entrained and transported from the surfaces when the wind speed reaches a critical threshold friction velocity and often assisted by saltation mechanisms [37,38].

Being unvegetated, flat, rich with unlimited supply of sediments, susceptible to wind erosion, and sensitive to anthropogenic disturbances, Lordsburg Playa (Figure 1) in Hidalgo County, New Mexico, United States; has been hotspot of dust events. This playa is crossed by Interstate 10 (I-10) and Union Pacific Railroad, making it the only playa hotspot in the United States bisected by an interstate highway [39]. It is usually inundated by highly variable seasonal rainfall events, mostly during summer North American Monsoon season, loading sediment into the playa through water erosion of upstream breached tanks and berms [40]. Hard crust surfaces are developed in the inner part of the playa as a result of very fine sediments and evaporation. An increase in the availability of loose sediments on the playa through anthropogenic disturbances (livestock and vehicular traffic, alteration of surface hydrology, depleted vegetation, and presence of exposed saline-sodic soils) routinely intensify dust generation [41]. The dust events from this surface have been causing major risks to traffic safety [40,42], public health [43], and economy [44]. For example, the New Mexico Department of Transportation (NMDOT) reported that dust storm (poor visibility)-linked traffic accidents have killed 21 and caused 39 closures at the Lordsburg Playa’s I-10 corridor between 2012 and 2020 [40].

Remote sensing imagery products have been invaluable resources for analyzing and monitoring LULC changes, hydrological responses, water erosion, and wind erosion [45,46,47,48]. Given their high temporal coverage, optical satellite remote sensing products have been widely employed for investigating playas’ long-term variability of inundation [9,49,50,51], water surface extent and area [52,53,54], water level [55,56], evaporative water loss [57], and spectral reflectance parameters [58]. However, few studies have been conducted to link the impact of hydrological changes to dust generation from playa surfaces [8,18,59,60]. The performance of optical satellite remote sensing is usually hindered during cloudy days, rainfall events, and nighttime. Due to its capability to penetrate through these obstructions, radar satellite remote sensing in combination with optical satellite remote sensing can significantly improve the characterization of playa surface responses to hydrological and LULC changes. Moreover, the signal detected by optical sensors into a single pixel is frequently a mixture of numerous disparate signals from different land cover types. This spectral mixing can occur for two reasons: if the spatial resolution is low enough and if the distinct endmembers are combined into a homogeneous mixture [61,62]. This type of mixing is most common in arid environments that are typically a combination of vegetation and soil in different proportions [63], and it can be partly addressed using spectral unmixing techniques [64].

Given this paradigm, the present study combines optical and radar satellite remote imagery products to characterize long-term changes of surfaces properties of Lordsburg Playa as a case study. It examines the fractional abundance of soil, vegetation, and water endmembers of each pixel of images using a spectral unmixing algorithm. This was done using optical remote sensing imagery products of Landsat and Sentinel-2 at Google Earth Engine. The study used the Landsat 8 and 5 at 30 m spatial resolution, and the Sentinel-2 Multispectral instrument at 10–20 m spatial resolution. It also studied the playa’s water level changes and the possible sinks for sediment loading and hotspots for wind erosion by employing Interferometric Synthetic Aperture Radar (InSAR) techniques. In this data recipe, the study analyzed a pair of Sentinel-1 images that bracket a monsoon day with high rainfall events and a pair of images representing a dry day and the monsoon day.

2. Study Area

The Lordsburg Playa is located in the far northwestern portion of the Chihuahuan Desert, a few miles west of the City of Lordsburg, New Mexico, USA (Figure 1). The closed basin containing the Lordsburg Playa is bounded to the east by the Pyramid Mountains and to the west by the Peloncillo Mountains [65,66]. This pluvial ephemeral lake has a lake area of 374 km², catchment area of 5670 km², lake area/catchment area ratio of 0.066, annual precipitation of 25 cm, and annual evaporation of 184 cm [67]. Thus, it usually experiences high hydrological stresses, because the annual evaporation rate exceeds the overall inflow rate to the basin floors from precipitation, runoff, and groundwater [67]. As a result of hydrological and geological processes which occurred during the pluvial Pleistocene epoch, sediments were deposited and formed into fine-textured soil associations in the basin between the mountains [68]. As a result, the soils of the piedmont slopes (including the edges of the playa) and the inner part of the playa are dominated by highly fragile and erodible soil classes [69,70]. The playa is usually inundated by highly variable seasonal rainfall, mostly occurring during the summer North American monsoon season, loading sediments through water erosion and breached tanks and berms [40,42]. It develops crusts on its surfaces once the water evaporates because the outflow rate is greater than the inflow rate. Anthropogenic disturbances in the playa increases the availability of loose sediments leading to the generation of dust storms [42].

3. Data

This study employed a combination of optical and radar satellite remote sensing imagery products and spectral signatures of endmembers to characterize the surface properties of the Lordsburg Playa. Satellite remote sensors acquire information about the Earth system through a variety of basic physical principles, detecting and measuring the electromagnetic spectrum by the radiation reflected (optical sensors), emitted (thermal infrared or passive microwave sensors), or scattered (active radar sensors including radar) [72]. Given these properties and their different atmospheric penetration capabilities, optical and radar satellite remote sensing data were synthesized to cover wider portions of the electromagnetic spectrum and, in turn, to improve the monitoring performance of LULC changes and hydrological processes [73,74].

3.1. Optical Satellite Remote Sensing Imagery, Inventory, and Preprocessing

Optical satellite remote sensing utilizes the visible, near infrared, and short-wave infrared parts of the electromagnetic spectrum to rasterize the Earth’s surface by detecting the radiation (mainly solar) reflected from the target on the ground. From this category, Landsat mission (Landsat 5 and 8) and Sentinel-2 images were used. Landsat products were selected for their long-term data coverage and Sentinel-2 images for their finer spatial resolution. Landsat 7 images were excluded from this study due to the data loss of 22–25% in each image resulting from the mechanical failure of the Scan Line Corrector (SLC) on 31 May 2003 [75]. This study utilized the invaluable resources of Google Earth Engine (GEE) [76] to avoid the space required for big satellite data storage; minimize the time required for data acquisition, downloading, processing, and analyzing; and exploit its powerful algorithms. GEE is a cloud-based platform that provides multi-petabyte analysis-ready satellite remote sensing data, internet-accessible application programming interface, and an interactive development environment that enable rapid prototyping and visualization, powerful algorithms, and capabilities for disseminating results to other stakeholders [76]. It is in use across a wide range of disciplines, including flood mapping, global forest change, global surface water change, crop yield estimation, rice paddy mapping, urban mapping, fire recovery, and malaria risk mapping [76].

The NASA and USGS joint program of Landsat series of earth observation satellites has continuously acquired images of the Earth’s land surface since 1972 at a 30 m spatial resolution about once every two weeks, providing invaluable resources for agriculture, engineering, geology, forestry, regional planning, mapping, and global change research [77,78]. This study used the combination of Landsat 5 (ETM sensor) and Landsat 8 (OLI/TIRS Sensors) Surface Reflectance Tier 1 images for Worldwide Reference System (WRS) Row 38 and Path 34 covering the Lordsburg Playa. Landsat images include seven bands with four visible and near-infrared (VNIR) bands, two short-wave infrared (SWIR) bands, and one thermal infrared (TIR) band. On the other hand, Landsat 8 scenes contain nine bands with five visible and VNIR bands, two SWIR bands, and two TIR bands. Both datasets are corrected for atmospheric effects.

Combining Landsat 5 and 8, the images cover almost 37 years of data records, Landsat 5 extending from March 1984 to May 2012 and Landsat 8 from April 2013 to present. The satellite overpass time (sensing time) over the Lordsburg Playa was around 17:12 UTC (between 16:45:23 and 17:40:50 UTC) for Landsat 5 and 17:47 UTC for Landsat 8 (between 17:44:43 and 17:50:31 UTC). A total of 637 images were collected, 487 images from Landsat 5 between March 1984 and May 2012, and 150 images from Landsat 8 extending from April 2013 to December 2019. The majority of the images showed very low cloud cover score (Figure 2a). After the screening, the study was left with 487 images representing 76% of the total images. The long-term monthly sum of cloud cover score distribution exhibited the mean value of 40.58% after the cloud screening, decreased by around 12.5% from the before-cloud screening scenario, which was 53.08% (Figure 2c). August, affected by the monsoon season, encountered a large drop in the number of images after cloud screening. The months in winter and late autumn experienced slightly fewer retained images compared with the remaining months before and after cloud screening. The collection, processing (including cloud screening), scaling, merging, and visualization of the images were performed in Google Earth Engine.

The European Space Agency (ESA)’s Multispectral Instrument (MSI) aboard two polar-orbiting Sentinel-2 satellites (Sentinel-2A and Sentinel-2B) provides global coverage (between latitudes 56°S and 84°N) at every 5 days with two satellites and 10 days with one satellite (European Space Agency, 2020). The MSI instrument detects the Earth’s surface using 13 spectral bands: four visible and NIR bands at 10 m spatial resolution, six red edge and SWIR bands at 20 m spatial resolution, and three atmospheric bands at 60 m spatial resolution. The current study used Level-2A images, which are atmospherically corrected Bottom of Atmosphere (BOA) surface reflectance derived from their corresponding Level-1C products [79]. This product temporally covers from 28 March 2017 to present, and each image consists of 12 bands (unlike in Level-1C, excluding cirrus band 10, as it does not contain surface information).

The study utilized Sentinel-2 images collected from March 2017 to December 2019 with the Orbit number of 141 and Tile Number Fields of T12SXA covering the Lordsburg Playa. A total of 77 images (39 images from Sentinel-2A and images from Sentinel-2B) images were acquired, 4 images during 2018 and 73 images during 2019 (Figure 2f). The satellite sensing time over the Lordsburg Playa was between 17:49:11 and 17:57:41 UTC for Sentinel-2A and between 17:49:09 and 17:57:39 UTC for Sentinel-2B. Majority of the images were not corrupted by clouds (Figure 2e). Before the cloud screening (<25%) the image counts were evenly distributed throughout the year, averaging at 6.42, except during December, in which all the 2018 images were acquired in December (Figure 2d). After the cloud screening, the number of images were reduced to 45, with the early and late months of the year showing a considerable drop. The collection of images from only 13 months (1 month from 2018 and 12 months from 2019) may have impacted the assessment of the cloud climatology. The collection and processing of all the images were performed in Google Earth Engine.

3.2. Spectral Signatures of Endmembers and Spectral Resampling

The spectral signatures of vegetation, soil, and water endmembers were acquired from the USGS spectral library with electromagnetic spectrum covering from the ultraviolet to the far infrared [80]. The measurements were acquired using laboratory, field, and airborne spectrometers for identifying and mapping materials on the earth and throughout the solar system through spectroscopic remote sensing. These measurements were already resampled to terrestrial multispectral sensors including Landsat 8 OLI and Sentinel-2 MSI. Thus, the resampled spectral signatures of the endmembers were directly applied to the linear spectral unmixing model. However, the original measurements were not resampled to the Landsat 5 ETM sensor wavelength spectrum. To fill this gap, we resampled the endmember spectra from Analytical Spectral Devices (which were convolved from the original measurements) to the wavelengths of the Landsat 5 ETM using ENVI version 5.3 (Exelis Visual Information Solutions, Boulder, Colorado) (Figure 3). In addition, the USGS spectral library does not include measurements from the Lordsburg Playa. To deal with this problem, we used spectral information of the endmembers from other playas across the southwestern United States (including Stonewall Playa, Cuprite, Nevada) with similar properties as endmembers from the Lordsburg Playa [80]. We used 6 bands from Landsat 5, 8 bands from Landsat 8, and 12 bands from Sentinel-2.

3.3. Radar Satellite Remote Sensing Imagery and Preprocessing

During cloudy days, nighttime, and heavy precipitation, the performance of optical satellite remote sensors is impaired, creating noisy pixels in the sensed images and gaps in the continuity of information flow. Given this limitation, it is crucial to integrate other data products with the capability to penetrate through those obstructions. Thus, in this study we also employed radar images from Synthetic Aperture Radar (SAR) remote sensing instrument aboard Sentinel-1 satellite mission, which has a capability of operating in all weather and all time scenarios. European Space Agency’s Sentinel-1 mission comprises a constellation of two polar-orbiting satellites (Sentinel-1A and Sentinel-1B), operating at C-band SAR (at 5.4 GHZ frequency or 5.6 cm wavelength) imaging providing a revisit time of 6 days for two satellites and 12 days for a single satellite [81]. Sentinel-1 C-band SAR collects images at a variety of spatial resolutions and polarizations.

To investigate the deformation and displacement of the Lordsburg Playa surface due to the inflow of water and associated water erosion from a single rainfall event and monsoon season, four images were acquired from the Alaska Satellite Facility’s Vertex data portal [82]. A pair of Sentinel-1 images (pre-event from 24 July 2017 and post-event from 5 August 2017) bracketing a wet monsoon day (1 August 2017) characterized by a high rainfall event on the region. We used a pair of images representing a dry day (pre-event from 14 March 2017) from the spring (dry) season and a wet day (post-event from 5 August 2017) from the summer (wet) season. We used Level-1 Single-Look Complex (SLC) product type from Sentinel-1B acquired through the Interferometric Wide (IW) Swath mode implying the spatial resolution of 5 × 20 m and dual polarization of VV and VH.

4. Methods

The methodologies utilized in this study are divided into two parts, the optical and radar satellite remote sensing imagery products (Figure 4). In the first part, we used the linear spectral unmixing model to identify the fractional abundance of endmembers in each pixel, water pixel flagging approach to estimate the areal extent of the water body, and Canny edge detection algorithm to delineate the edge of the water body in the playa. In the second part, we applied the Interferometric SAR (InSAR) technique to estimate the deformation and displacement of the playa due to changes in hydrology.

4.1. Linear Spectral Unmixing, Water Pixel Flagging, and Edge Detection

To estimate the fractional abundance of the playa’s land surface endmembers, we used the linear spectral unmixing that decomposes mixed spectra into distinct endmembers. Equation (1) represents the linear unmixing model:

$$y_i=\sum _{j=1}^3{x}_{ij}{\beta }_j+{\epsilon }_i$$

(1)

where $y_i$ denotes the spectrum mixture of each pixel at each ith spectral channel from the optical remote sensing imagery; $x_{ij}$ denotes the spectral signature of endmember j at ith spectral channel from either laboratory, field, or airborne spectrometers or combinations of these instruments; ${\beta }_j$ denotes the fractional abundance of endmember at ith spectral channel; ${\epsilon }_i$ denotes the modeling residual error at ith spectral channel; j denotes the number of endmembers (in this case, 3 endmembers namely, vegetation, soil, and water); and i denotes the number of the spectral bands. In addition, because the fractional abundance parameter $\beta$ represents the fractional areal abundances shared by the endmember j, we constrained the fractional abundance values by the nonnegativity (Equation (2)) and sum-to-one (Equation (3)) conditions. The unknown parameters (${\beta }_j$) can be estimated by using different linear algebra approaches, which minimizes the residual error term (${\epsilon }_i$). However, in this study, we used the spectral unmix algorithm provided by the Google Earth Engine. This approach decomposes each mixed pixel by computing the pseudo-inverse and multiplying it through each pixel. The fractional abundance of water was further used for the estimation of playa’s flooded water extent.

$${\beta }_j\ge 0,forendmembers\left(j\right)vegetation,soil,andwater$$

(2)

$$\sum _{j=1}^3{\beta }_j=1$$

(3)

In order to estimate the areal extent of the flooded water in the playa, we applied the water flagging approach. The pixels were flagged as water when the fractional abundance threshold of water endmember was ≥ 0.35 [83]. For example, Du et al. [84] showed that the maximum water mapping threshold was between 0.35 and 0.4. Furthermore, we also used the Canny edge detection algorithm to locate the edge of playa’s water body due to direct rainfall and runoff from its watershed. This analysis was also performed in Google Earth Engine, in which the output is an image with non-zero value indicating edges of the water body, and the magnitude of the value is the gradient magnitude. This digital image processing algorithm has been widely employed in different edge detection studies and compared with other algorithms Canny approach showed better performance [85].

4.2. Deformation and Displacement Estimation

Using the Synthetic Aperture Radar (SAR) images from Sentinel-1B satellite, we investigated the deformation and displacement of the Lordsburg Playa surface. The analysis involves two steps as depicted in the flowchart (Figure 4). The first step employs InSAR technique to measure the changes in land surface deformation using a pair of SAR images from Sentinel-1B satellite. This imaging technique measures the topography of te land surface, its changes over time, and changes in the detailed characteristics of the surface [86]. It exploits such changes by measuring the difference between the phase signals of the SAR images. All the analysis stages in this step were accomplished in Sentinel-1 Toolbox (S1TBX).

As a pre-processing stage, duplicate SAR images were first co-registered into a stack in order to properly align the pixels of the images. Then, interferogram phases were created by cross-multiplying the pre-event image with complex conjugate of the post-event image. In order to perform interferometric processing, we used two pairs of images covering the Lordsburg Playa acquired at different satellite over pass times. One pair of images consisted of one pre-event image from 24 July 2017 and one post-event image from 5 August 2017, bracketing a wet monsoon day of 1 August 2017. Another pair of images represented a dry spring day of 14 March 2017 and a wet summer monsoon day of 5 August 2017. Following that, the noise introduced by the previous stage was minimized by applying multi-looking and phase filtering techniques. Finally, the processed images were geocoded (projected into a geographic coordinate system) and exported using a preferable file format. In the second step, using the SNAPHU algorithm [87], the generated interferogram phase images were unwrapped to extract displacement.

5. Results

5.1. Fractional Abundance of the Endmembers

The maps for the fractional abundance of soil, vegetation, and water endmembers are presented in Figure 5, Figure 6 and Figure 7, respectively. These maps provided a means to discern surface cover proportions of each pixel. The maps captured the natural distribution of the land cover across the playa. These fractions were retrieved from Landsat 5 and 8 remotely sensed satellite images using linear spectral unmixing approach for the spring seasons of 1984 to 2019. Similar outputs for the summer/fall seasons from the Landsat images (Figure A1, Figure A2, Figure A3), and monthly maps of 2019 from Sentinel-2 images (Figure A4, Figure A5, Figure A6) are provided in the Appendix A.

For the soil endmember map (Figure 5), the light brown areas represent higher fractions, and deep blue areas represent lower fractions of soil cover. White color indicates masked pixels due to cloud obstruction. The playa areas were dominated by highest fractions of soil endmembers during spring in the majority of years. This suggests that during these years, large parts of the playa surfaces, including the edge of the playa and its surroundings, could be exposed to wind erosion. For example, in spring 2011, a year with high number of traffic crashes due to dust emissions [69], the playa pixels were dominated by the highest fractions of soil endmember. Low fractions of soil were mainly exhibited in the inner and low-lying parts of the playa for a few years, primarily influenced by hydrological processes (monsoon inundations) during a current season and/or accumulated from the preceding seasons. In 2015 to 2019, soil cover in the majority of pixels exhibited a decrease in soil fractions. This decrease in exposed soil surface could contribute to the decrease in the number of traffic accidents in those years associated with a reduction in dust events on the playa [69]. On the other hand, summer/fall seasons experienced a decrease in the fractions of soil endmember across all years as compared to corresponding spring seasons. This is due to rainfall input to the hydrological system of the playa during the North American monsoon.

Vegetation endmember retrieved from the Landsat 5 and 8 images, in contrast to soil endmember, contributed low fractions to pixel mixtures across the playa and throughout the majority of spring seasons (Figure 6). Only the spring seasons from 1995, 1998, and 2005 experienced the highest fractions of vegetation cover. Those pixels with highest vegetation fractions were from the uplands bounding the playa on the west and east. For 1998 and 2005, the pixels were only from the Peloncillo Mountains to the west. Thus, during the majority of the years, the playa surfaces were not covered by active vegetation, potentially exposing the surface to wind erosion. In summer/fall seasons, the proportion of vegetation endmembers in the pixel mixtures were highest, giving the playa surface and surrounding some protection against erosion (Figure A2). The 2019 monthly vegetation fraction developed from the Sentinel-2 images was lowest from May to August and highest in September and October (Figure A5). Although in this region July and August are the wettest monsoon months, the phase lag in the response of vegetation to rainfall input influenced the vegetation cover availability across the summer and fall months.

The fractional abundance of the water endmember for the spring seasons from the Landsat 5 and 8 images captured the playa’s water dominated pixels (Figure 7). The pixels from the inner and low-lying part of the playa showed higher fractions. The map was also able to successfully capture water-dominated pixels from the very small water bodies along the sides of the major transportation systems (ditches) from the spectral mixtures. Complementing these results, using the fractional abundance of water endmember as input, the areal extent of the flooded water within the playa was identified using the fractional abundance threshold of 0.35 (Figure 8). This approach aided the detection of very small, shallow, and turbid water bodies due to the coarse spatial resolution of the optical images used in this study. It would have been very difficult to detect these types of water bodies using the widely used spectral indices such as the Normalized Difference Water Index (NDWI) [53,88]. For example, Herndon et al. [53] documented the failure of NDWI to detect many smaller, more turbid water bodies in the study area of Tahoua Region of Niger. Supporting this finding, Acharya et al. [88] also pointed out that NDWI showed better results for only pure water pixels. In contrast, the highly seasonal water bodies in the Lordsburg Playa are shallow and turbid (Figure 1). The water bodies are usually created as a result of strong summer monsoon seasons and can sustain standing water until the following spring. Similarly, the maps from Landsat images for summer/fall seasons captured the flooding water with larger areal extent compared with those from spring seasons (Figure A3). The 2019 monthly fractional abundance maps from the Sentinel-2 images provided similar results to the Landsat image (Figure A6), although the corruption of the maps by clouds was higher in Senintel-2 images. Thus, the linear spectral unmixing algorithm demonstrated its potential for extracting different surface properties from pixel spectral mixtures of optical images.

5.2. Synthetic Aperture Radar (SAR) Images

The results from the synthetic aperture radar (SAR) images from Sentinel-1 satellite products identified locations prone to phase and displacement changes via water level change. These locations may become hotspots of dust generation during dry seasons due to sediment availability. The output from the interferogram phase analysis revealed that the surface deformations extracted from pairs of images was primarily due to hydrological inputs. The deformation between the pair of images from spring (acquired on 14 March 2017) and summer (acquired on 5 August 2017) was due to the seasonal rainfall input and sediment loading transported from the watershed onto the playa through runoff (Figure 9a). The deformation between paired images of 24 July 2017 (pre-event) and 5 August 2017 (post-event), bracketing the rainfall event of 1 August 2017, was mainly due to the event’s rainfall and runoff causing sediment loading on the playa (Figure 9b). The movement of the edge of the playa’s water body between dry and wet seasons increases the availability of sediments for wind erosion along the edge of the playa [7,43]. The displacement maps (Figure 10) that were developed from the interferogram phase also highlighted the locations experiencing elevation change due to rainfall and sediment loadings. The highest positive elevation change reached approximately 5.0 cm between spring and monsoon seasons (Figure 10a) and up to around 5.4 cm due to a single monsoon rainfall event (Figure 10b) at the inner part of the playa. The lower elevation changes for the pair of images between the spring and summer than for the pair of images from summer indicates the intensity of evaporative loss in this region. Thus, the combination of optical and radar satellite remote sensing products aided us in assessing the horizontal and vertical changes within the playa and improved the mapping of hotspots susceptible to wind erosion as well as inundated areas, which would not be susceptible to erosion.

6. Discussion

In summary, the assimilation of optical and radar satellite remote sensing imagery products improved discernment of playa surface changes in space and time. The identification of regions vulnerable to these fluctuations assisted the process of delineating hotspots susceptible to wind erosion. Based on these delineations, land managers can plan their strategy in controlling wind erosion and establish sustainable field-level conservation plans. The utilization of Google Earth Engine facilitated acquisition, processing, and analysis of the optical images by providing cloud computing platform, long-term satellite imagery, and powerful algorithms. The application of linear spectral unmixing algorithms to the Landsat and Sentinel-2 optical images discriminated the contribution of each endmember (soil, vegetation, and water) into the pixel mixtures. The hydrological processes (monsoon rains, dry seasons, and playa inundation) occurring during the specific season and the preceding seasons were the main factors influencing the composition of each pixel.

During the spring seasons, the playa pixels were dominated by high fractions of soil endmember, except for inner low-lying parts of the playa. Due to this dominance, the majority of the playa surface was potentially exposed to wind erosion. Vegetation endmember showed low fractions across the playa and throughout the majority of the spring seasons. However, during summer/fall seasons, the proportion of vegetation endmember in the pixel mixtures was improved, giving the playa surface some protection against wind erosion. Water endmember exhibited high fractions in the inner part of the playa. The linear spectral unmixing algorithm demonstrated its capability by capturing the water-dominated pixels from the very small water bodies (ditches) along the sides of major transportation systems crossing the playa. Using the water fractions as an input, the pixels with greater than or equal to 0.35 threshold were flagged as water, and thus the areal extent of the flooded water within the playa was determined. By using SAR images from Sentinel-1, land deformation and displacement due to hydrological processes were estimated. The locations prone to phase and displacement changes due to water level change were identified and related to possible dust hotspots. Therefore, assimilation of optical and radar satellite remote sensing products significantly improved information about the long-term evolution of the playa and mapping of hotspots susceptible to wind erosion.

Satellite imagery products have been invaluable for large-scale investigation of historical LULC and hydrological changes. However, spectral indices retrieved from moderate spatial resolution satellite imagery products such as Landsat and Sentinel-2 have limited applications in detecting small-scale changes in areas dominated by highly seasonal and shallow water bodies like the Lordsburg Playa [53,88]. This is because many those water bodies possess areal coverage finer than the spatial resolution of the satellite imagery considered in this study. These water bodies and vegetation cover are critical factors that determine the availability of soil particles for wind erosion in global playas such as Lordsburg, and identifying these land cover types is important for land managers and policy makers in mitigating dust emission. In addition, for assessing the vegetation health in semi-arid regions, the spectral unmixing approach showed better performance compared to the Normalized Difference Vegetation Index (NDVI) due to being less sensitive to soil background [89].

This study identified the locations within the playa exposed to long-term hydrological and LULC changes, which are critical parameters that determine the sediment supply for wind erosion. The application of linear spectral unmixing algorithm and InSAR technique helped us to identify patterns and features of changes buried in the optical and radar imagery products, which are reflected in the final results of this study. Accordingly, the seasonal hydrological inputs and changes have paramount importance in determining the Lordsburg Playa’s capacity in generating and emitting dust, and in turn traffic safety.

7. Conclusions

The results of this study illustrate how introducing cloud computing platforms and accompanying algorithms into able to on the Earth system over the past decade has facilitated and improved the mining of LULC and hydrological changes buried within large volumes of satellite remote sensing imagery products [90]. It explored the potential of spectral unmixing in addressing limitations of moderate spatial resolution optical remote sensing imagery and discriminating the contribution of surface endmembers into the spectral mixture.

This study took advantage of long-term records of optical imagery; cloud-, nighttime-, and precipitation-insensitive radar imagery, and cloud-based cyber infrastructure to discern LULC and hydrological changes in the dust-emissive Lordsburg Playa. While the linear spectral unmixing algorithm demonstrated its capability to discriminate the surface endmembers from spectral mixture, non-linear approaches, including machine- and deep learning methods could also be used to accurately identify the fractional abundance of endmembers, because the LULC and hydrological variables are non-linear in nature. The spectral unmixing approach detected numerous small, shallow water bodies, including those alongside major transportation systems which are difficult to detect with traditional spectral indices due to their coarser spatial resolution. The study also showed that vertical changes within the playa can be estimated using radar images, compensating for limitations of optical remote sensing. Complementing SAR images with optical images aided investigation of both horizontal and vertical changes within the playa. The application of field and airborne spectrometers for the measurement of spectral signature of the endmembers over highly dynamic land surfaces such as playas should be further explored for potential improvement of the performance of both linear and non-linear spectral unmixing models. The detection of seasonally exposed ephemeral lake surfaces in greater detail can significantly assist the quantification of potential soil loss from playas through erosion and the resulting dust emissions’ impact on air quality and nearby infrastructure. Identification of locations sensitive to seasonal land use/land cover changes and hotspots for aeolian processes are valuable findings for supporting the efforts towards erosion control plans.