Non-Vegetated Playa Morphodynamics Using Multi-Temporal Landsat Imagery in a Semi-Arid Endorheic Basin : Salar de Uyuni , Bolivia

Playas in endorheic basins are of environmental value and highly scientific because of their natural habitats of a wide variety of species and indicators for climatic changes and tectonic activities within continents. Remote sensing, due to its capability of acquiring repetitive data with synoptic coverage, provides a unique tool to monitor and collect spatial information about playas. Most studies have concentrated on evaporite mineral distribution using remote sensing techniques but research about grain size distribution and geomorphologic changes in playas has been rarely reported. We analysed playa morphodynamics using Landsat time series data in a semi-arid endorheic basin, Salar de Uyuni in Bolivia. The spectral libraries explaining the relationship between surface reflectance and surficial materials are extracted from the Landsat image on 11 November 2012, the collected samples in the area and the precipitation data. Such spectral libraries are then applied to the classification of the other Landsat images from 1985–2011 using maximum likelihood classifier. Four types of surficial materials on the playa are identified: salty surface, silt-rich surface, clay-rich surface and pure salt. The silt-rich surface is related to crevasse splays and river banks while the clay-rich surface is associated with floodplain and channel depressions. The classification results show that the silt-rich surface tends to have a positive relationship with annual precipitation, whereas the OPEN ACCESS Remote Sens. 2014, 6 10132 salty surface negatively correlates with annual precipitation and there is no correlation between clay-rich surface and annual precipitation. Salty surfaces seem to consist primarily of clay due to their similar characteristics in response to precipitation changes. The classification results also show the development of a crevasse splay and avulsions. The results demonstrate the potential of Landsat imagery to determine the grain size and sedimentary facies distribution on playas in endorheic basins.


Introduction
As a common feature of internal drainage basins, playas are vital for understanding the natural habitats of endemic microbes, plants and animals [1] and for interpreting the interrelation between climate and tectonic activity within continents [2].Studies also proved that playas sediments could be potential hydrocarbon reservoirs [3], as they may contain significant numbers of thin fluvial sheetflood deposits.Many researchers have suggested that a framework for surface sedimentary facies or depositional environments is needed in playas [4][5][6][7].Over the past decades, a number of studies have applied remote sensing techniques to study playas and have demonstrated the usefulness of remote sensing in such studies (e.g., [1,2,[8][9][10][11][12][13][14][15]).However, most of these studies have focused on mapping evaporite minerals and water bodies.There are few studies focusing on grain size and sedimentary facies distribution using remote sensing [10,[16][17][18][19][20], mostly due to difficult accessibility [2,18] and hazardous environments during peak discharge events [21].In this paper, we analyse playa morphodynamics by Landsat time series data in the semi-arid endorheic (internal drainage basin) basin of the Salar de Uyuni in Bolivia.
A basic theory about the relationship between sediment and reflectance is used in this paper.If the sediment is very coarse, the porosity is high and thus the scatter effect is strong and the reflectance is low.The high porosity also indicates a high capacity for storing water.As grain size decreases, the water content simultaneously decreases due to a decline in porosity.When the sediment is fine sand or silt, the reflectance is high due to little scattering effect and a large specific area.However, with an increasing content of clay, the capacity of the sediment for storing water increases because of the large internal surface, which is related to adhesive, cohesive and osmotic forces [22].Due to strong evapo-transpiration in the study area (1500 mm/year), the water in the sediments precipitates salt in the long dry period (normally from April to November in the study area).When the surface is covered by salt it has a strong reflectance, while for clay-rich areas, the reflectance is low due to soil moisture.In addition, crust cracking at clay-rich surface reinforces its roughness and therefore contributes to a low reflectance.
This study documents the use of Landsat imagery in combination with precipitation data and field investigation to analyse playa morphodynamics.First, we establish the relationship between the various types of sediments (e.g., salt-, sand-and clay-rich sediments) and spectral reflectance, and then build spectral libraries for each type.Such spectral libraries are then used to classify the Landsat time series data from 1985-2011 over the study area, with Maximum Likelihood Classification (MLC) method being used as a supervised classifier.Then, we analyse the relationships between different classes and sedimentary facies, and investigate the changes of these classes in combination with precipitation data.Finally, we analyse the geomorphological changes with a time series of classified maps.

Study Area
The study area is at the southeast edge of the world's largest salt lake, Salar de Uyuni, Bolivia, with an area of 10,582 km 2 and it is at an elevation of 3656 meters above mean sea level (Figure 1).It is located in the southern part of the Altiplano Basin.The basin formed as part of the Andean ocean-continent convergent margin with eastward subduction of the oceanic Nazca Plate beneath the continental South American Plate [23,24].The central Andes developed eastward during the Cenozoic and regional horizontal shortening led to increasing thickness of the continental crust [25,26].The Altiplano Basin, an endorheic basin, is filled with Tertiary to Quaternary fluvial and lacustrine sediments and volcaniclastic deposits [27,28].The Altiplano Basin has a semi-arid climate with an annual precipitation of more than 800 mm in the north and less than 200 mm in the south [30] and an evapo-transpiration potential of 1500 mm [31,32].Precipitation shows this clear-cut north-to-south decrease in Altiplano plateau due to the low pressure systems of the Altiplano, strong low-level north-westerly winds with warm, moist, unstable air flow along the eastern flank of the central Andes giving rise to convection precipitation.Moreover, pole-ward low-level air flow helps to maintain the intense convection [33].The rainy season in the study area is from December to March, and the average annual rainfall is about 185 mm, which is greatly exceeded by potential evapo-transpiration.With SE-to-NW flow direction, the Río Colorado is one of the main inflow rivers for the Salar de Uyuni.In this study, we focus on the terminus of Río Colorado river system, where human interference is limited.

Data Acquisition
Daily precipitation data were available over the study area for the period 1985-2012 from the Bolivian Servicio Nacional de Meteorología e Hidrología (Figure 2).Using the precipitation datasets from three meteorological stations in the study area over the period 1980-2010, Li et al. [21] has demonstrated that the rainfall distribution in Uyuni is spatially consistent with those of other meteorological stations using the double mass curve technique.This implies that such daily precipitation data in Uyuni can be considered representative of the whole study area.The field survey was carried out from 13-18 November 2012 before the rainy season.The field sampling sites were determined using Landsat ETM+ imagery and land cover samples were collected with geo-located GPS.Landsat Surface Reflectance Climate Data Record (CDR) data was collected from the USGS website.The surface reflectance CDR, including many Landsat TM and all ETM+ scenes, is a high-level Landsat data product and is derived from specialised software called Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS).The surface reflectance data were available in the form of 16-bit signed integer with the range of −2000 to 16,000.The studied river terminus is included entirely within Landsat path 233, row 74.Thanks to the available daily precipitation data, we accurately determined the rainy and dry periods of 1985-2012.This allowed us to select Landsat images in the dry periods (Table 1).

Landsat CDR Data Pre-Processing
As seen in Figure 3, first of all, image-to-image registration is done to make sure that a certain pixel covers the same area in all the images.The Landsat surface reflectance data are then normalized with Iteratively Reweighted Multivariate Alteration Detection (IR-MAD) transformation to eliminate the unwanted effects due to illumination variation for different dates, although, surface reflectance data as used in this study (Landsat CDR) are supposed to be corrected for such effects which implies that such normalizations may not be strictly required [34].The false-color composite of the base image in combination with field data analysis were used to select training samples.Prior to classification, it is necessary to check the separability of the target classes.This has been done using the transformed divergence technique, where a low value between two classes indicates similarities between them and hence difficulties in distinguishing two distinct classes by the classifier.The maximum likelihood classifier was then applied for the classification.This is a well-known approach that has been frequently used for classification in semi-arid environments [35][36][37].The accuracy assessment of the approach was evaluated by different metrics including overall accuracy, kappa coefficient, producer's accuracy and user's accuracy.More detailed explanation of the methodology will be given in the following subsections.
Spectral libraries of different classes (end-members) were built using the Landsat image on 11 November 2012 and ground truth data collected during the field campaign (Figure 3).This image was selected as the base image since this is the closest date corresponding to the ground measurements and there was no TM image available for this year.Bands with a low fraction of saturated pixels have been selected for later normalisation and classification.In this paper, we used three Landsat bands 2, 4 and 7, which have low fractions of saturated pixels for the period 1985-2011 according to digital number statistics of all bands.All these selected bands were subsequently normalized to a referenced image with IR-MAD transformation [38].The IR-MAD requires a full dataset since it uses the spatial data over the entire scene for the normalization.Therefore, the ETM image of 2012 due to the missing lines (Scan Line Corrector (SLC) failure) cannot be used and we have taken the TM image of 2011 as the reference image for normalization of the other images.This might introduce error into the analysis, however, since the data are surface reflectance (Landsat CDR data) that is supposed to be corrected for illumination, sensor-related and other effects, this error will not be significant over such a homogeneous area.Moreover, the high coefficients of determination (mean R 2 > 0.95) indicated strong correlations between the reference image and targeted images.

Training Samples Selection and Analysis
The training samples were selected from a false-color composite image and field sampling data.We used false-color composites of band 7 (2.08-2.35µm), band 5 (1.55-1.75µm) and band 1 (0.45-0.52 µm) (Figure 4) to differentiate different surface materials (in particular clay-rich and silt-rich surface).Since both bands 5 and 7 are located in the near-infrared region, they are sensitive to soil moisture and indicate water absorption regions.We identified four surface materials: A, B, C and D (Figure 4).
The field campaign included surface color analysis and sampling.Surface color analysis included the four documented types of surface materials based on false-color composite map (Figure 5).The colours of different surface deposits were determined using Munsell colours chart in the field: A is from 2.5 YR 8/2 to 5 YR 8/2 (light pink); B is from 2.5 YR 8/4 to 5 YR 8/4 (pink); C is 7.5 YR 6/3 to 7.5 YR 5/6 (light/strong brown).The samples were then analysed in terms of grain size with a Sympatec HELOS KR laser-diffraction particle sizer with a size range from 0.1-2000 µm.The results of the grain size analysis of field samples were used to validate the analysis of the false-color composite Landsat images on 11 November 2012 (Figure 6).We assigned the following names to the classes: A, salty surface; B, silt-rich surface; C, clay-rich surface; D, salt.Here, salty surface represents mixed sediments of salt and silt or clay with a light pink colour, while the silt-rich surface is coarser than the clay-rich surface.Salt indicated a thick, white salt crust on the surface.2) and the other for ground truth accuracy assessment.The training samples were analysed in terms of separability between classes with transformed divergence method [39].
where i and j represent the two signatures (classes) being compared, and T is the transposition function.
The index varies from 0 and 2. When the value is greater than 1.8, the compared pairs have good separability, whereas when it is less than 1, they should combine into one single class [40].The transformed divergence indices were greater than 1.97 for most of the classes (Table 3).These samples were then used as input for the MLC method.

Maximum Likelihood Classification
Evans [41] has proved that the prior knowledge of the relationships between input attributes and salinity in the classification of Landsat TM data is of importance for salinity mapping.We used a supervised classification, the MLC method, which employs discriminant functions to assign pixels to the class with the maximum likelihood [42].

Accuracy Assessment
The classification results were evaluated by different metrics including the overall accuracy, kappa coefficient, producer's accuracy and user's accuracy (Table 4).Due to the limited number of ground samples, the cross validation approach was implemented to make use of all the samples in both training the model and for the validation.The average of different runs was then used as the ultimate measure of the classification.The overall accuracy was 92.83% and the Kappa coefficient was 0.90.

Application to the Other Landsat Images Processing
The established spectral library based on the analysis of Landsat scene of 11 November 2012 was applied to other Landsat scenes for the period 1985-2011.We expanded the spectral libraries and randomly selected more than 1000 pixels with more than 10 polygons to training samples based on the classified map on 11 November 2012 (Figure 7).We used the toolbox "Endmember Collection" to perform classification on Landsat images from 1985-2011 with the remote sensing image processing software ENVI.

Areal Statistical Analysis of the Different Classes
Surface materials change over years in terms of the area that they cover (Figure 8).The clay-rich surface is the most dominated surface material with a few exceptions over the period 1985-2011 (Table 5) covering on average 40% of the study area.According to the standard deviations, the salt remained consistent over time (with standard deviation of 21.48 km 2 ) compared to the other surface materials in the study area.However, in the river terminus, all the three surfaces are characterized with large variations over time with standard deviations of 38.37 km 2 , 41.96 km 2 and 41.27 km 2 for salty, silt-rich and clay-rich surfaces, respectively.However, the classification results of the year 1996 showed very large discrepancies compared to the other years for the three surfaces which seemed to be abnormal in this case.We could not identify any specific reason for this discrepancy.

Interpretation of the Different Classes
In the river terminus, the silt-rich surface and the clay-rich surface tend to be associated with different sedimentary facies.For example, the silt-rich surface is distributed on crevasse splays and levees along channels, while the clay-rich surface can be found on the floodplain between channels.In addition, the silt-rich surface was generally distributed in the proximal part of the river terminus while the clay-rich surface was commonly seen in the distal part of the terminus.The clay-rich surface sometimes could be seen in the salt lake, probably attributed to peak floods or wind (Figure 9).During peak floods, the water is at its highest level in the salt lake and clay is transported to the salt lake and the edge of the water body.When the water evaporates and salt is precipitated, the clay in the deep water area is covered by salt but clay at the edge of the water body is mixed with salt and exposed.Depressed channels were also covered by the clay-rich surface, although sometimes salt could be found in the present channel (Figure 9).The salty surface was generally distributed in the channels and the distal part of the river terminus, while the salt was generally distributed at the end of the present channel and the salt lake.In combination with the annual precipitation, we analysed the areal changes of each class between 1985 and 2011.The silt-rich surface shows a weak correlation with the annual precipitation (Figure 10A) probably because although the annual precipitation is high, the daily precipitation all year around is low.However, some daily precipitation is high in those low annual precipitation periods and Li et al. [21] have demonstrated that high daily precipitation would lead to peak discharges, which result in large flow and sediments transport to the river terminus.During these short-lived and high magnitude of peak floods or multiple successive floods, channels experience massive over-spilling flow due to the downstream decrease in cross-section area on the low-gradient Río Colorado river system [21,43] and the spill-over flow reactivates crevasse splays and levees.For instance, the crevasse splay was partly reactivated in the rainy season of 1999-2000.However, it was completely reactivated in the rainy season of 2000-2001 because of the highest annual precipitation (406.4 mm) in the period from 1985 to 2011 (Figure 11).
By contrast, the salty surface was negatively correlated with the annual precipitation (Figure 10B), i.e., when the annual precipitation increased, the salty surface decreased.For the clay-rich surface, there is no correlation with the annual precipitation (Figure 10C).The area of river terminus, which includes salty surface, silt-and clay-prone surface, tended to be inversely proportional to the annual rainfall (Figure 10D).This is likely due to the fact that the higher the annual rainfall, the more water is draining to the lake, which results in high water levels in the lake.When the water evaporates, salt is precipitated on the surface.Therefore, the salt area will increase and the river terminus shrinks when the annual precipitation increases.

Geomorphological Change Detection
Geomorphological changes can be visualized with comparisons of classification maps.The positive correlation between silt-rich surface and crevasse splays enables us to characterise the expansion of crevasse splays over the observation years.Specifically, the development of a crevasse splay in the distal part of the river terminus was recorded through time-series classification maps.It started with a small crevasse splay and then grew into large splays during the subsequent floods (Figure 12).The expansion of the crevasse splay is in agreement with that by Li et al. [21], which also showed that this new crevasse splay was attributed to a compensational stacking pattern, whereby the new crevasse splay expanded in the topographic low between two adjacent existing splays.The expansion of crevasse splays also contributed to the increasing area of the silt-rich surface.The avulsion history was reconstructed by comparing classification maps (Figure 13).The observed avulsion was initiated with a crevasse splay and then developed into an avulsion with subsequent floods.Previous studies suggested that the splay-induced avulsions are attributed to downstream decrease in channel cross-section area in the dryland river system [21,43].Therefore, rivers in the low-gradient dryland region experience massive flow over-spilling onto the floodplain, which reinforces levee breaches during peak floods, and increases the chance of bank failure.In addition, the sediment-laden water would flow back to the trunk channel during the waning flood stage.The returning flow erodes the crevasse floor and deepens the crevasse channels; on the other hand, the process decreases the capacity of the trunk channel because sediments in the returning flow are deposited at the downstream side due to low energy of returning flow [43].These processes promote avulsions and simultaneously increase the silt-rich surface in the distal part of the river terminus.

Discussion
Due to salinity, the Time Domain Reflectometry (TDR), which is used to simultaneously obtain soil water content and electrical conductivity by a single probe with a minimal disturbance [44], did not work properly.However, field observation showed that the clay-rich surface is wetter than the silt-rich surface (Figure 14A).In addition, the density of desiccation cracks (Figure 14B) is also higher on the clay-rich surface than on the silt-rich surface.Both soil moisture and desiccation cracks exert a reduction on reflectance.In order to check the consistency of the surface reflectance for different years over the study area, we selected an area of 70 pixels located in the centre of a clay-rich surface on both 11 November 2012 and 30 September 2011 images (Figure 15).The tested area is in the centre of a floodplain area (Figure 9) which is marginally influenced by peak discharge.As seen in the figure, the extracted reflectances show good agreement between the two dates.However, some differences are still seen between the two images, partly due to the differences in the sensor configurations and partly due to the pre-processing of the images or changes in the area over time.In general, the absolute difference between the two reflectances is mostly lower than 3% (reflectance), which can be considered as a reasonable threshold for classification.
By comparing the reflectance of different classes (see Figure 7), it can be seen that the mean values are very different for each class, though they show similar trends.This might lead to the conclusion that classification methods based on angular mapping would work better in such a case.In order to demonstrate this, we have implemented other classification methods based on angular mapping like Spectral Angle Mapper (SAM) and as well as Support Vector Machine (SVM).However, these two methods showed lower accuracies (not shown) compared to the maximum likelihood method and failed in producing reasonable classification maps.Moreover, the classification results show the tested area is classified the same (i.e., silt-rich surface) in the two images using maximum likelihood method demonstrating that the good performance of the method.
We also suggest that post-processing such as majority/minority analysis, sieve classes should not be performed on classification maps.First of all, the inter-channel deposits are generally clay-rich and, therefore, their reflectance is low.In addition, the inter-channel area is sometimes small (e.g., crevasse channels) with the ground resolution of the Landsat TM images being 30 m by 30 m.Using these post-processing methods tends to blur these sedimentary phenomena and create errors.Therefore, the statistical analysis on playa morphodynamics was found the result of maximum likelihood classification to be better than that with post-processing.Surface reflectance from false-color composite images was used to differentiate surface materials, and these surface materials correspond to a particular color (Figure 5), but still with some errors.For example, some sampling locations apparently look pink and samples from these places are supposed to be silt-rich, but the grain size analysis shows that these samples belong to clay-rich materials.

Conclusions
This paper contributes to the understanding of playa surface dynamics using remotely sensed data over a period of 27 years, from 1985-2011.The study investigates potentials of multi-temporal Landsat surface reflectance in mapping coarse and fine sediments in the river terminus of the Salar de Uyuni, Bolivia.First, spectral libraries were generated based on the Landsat surface reflectance image of 2012 and field data analysis corresponding to the ground truth data for four different surface materials in the river terminus.These surface materials were salty, silt-rich, clay-rich and salt surface.Such spectral libraries were then used to classify the Landsat image using Maximum likelihood method.The cross validation approach was implemented to validate the classified map.The overall accuracy of 92.83% and the Kappa coefficient of 0.90 were achieved.Second, the consistency of the reflectance over time for different images was investigated.For instance, the absolute difference between reflectances of the images of 2011 and 2012 were found to be mostly lower than 3% (reflectance) which was considered as a reasonable threshold for classification.Finally, the spectral libraries were used to classify the other Landsat images from different years over the study area.
According to the classified maps, the clay-rich surface remained the most dominant surface material with a few exceptions over the entire period covering approximately on average 40% of the study area.The salt was relatively stable in terms of area over time (with standard deviation of 21.48 km 2 ) compared to the other surface materials in the study area.However, in the river terminus, all the three surfaces are characterized with large variations over time with standard deviations of 38.37 km 2 , 41.96 km 2 and 41.27 km 2 for salty, silt-rich and clay-rich surfaces, respectively.The silt-rich surface and the clay-rich surface are associated with different sedimentary facies in the river terminus.The silt-rich surface was widely distributed along channels and crevasse splays while the clay-rich surface was found in the floodplain and the channels.In addition, the proximal part of the river terminus was dominated by the silt-rich surface while the distal part of the terminus was mainly occupied by the clay-rich surface.The clay-rich surface occurs in the salt lake during floods or through wind-blown influence.The silt-rich surface has a weak correlation with the annual precipitation while the salty surface tends to be inversely proportional to the annual precipitation.High annual precipitations have better chances to reactivate rivers in the playa and lead to over-spilling flows onto levees and crevasse splays due to downstream decrease in bankfull width and depth.Geomorphological changes have been identified such as the formation of crevasse splays and avulsions.High annual precipitation-induced avulsions also led to increased silt-rich surfaces.

Figure 1 .
Figure 1.Map of the study area.(A,B) The location of Altiplano plateau in South America (modified from [29]).(C) A false-color image of Landsat Thematic Mapper (TM) with bands 7 (red), 5 (green) and 1 (blue) on 24 February 2006 indicating the study area.

Figure 3 .
Figure 3. Flow chart for the Landsat images processing.

Figure 4 .
Figure 4. False-color composite of bands 7 (red), 5 (green) and 1 (blue) of Landsat ETM+ on 11 November 2012.A, B, C and D represent different surface materials.The black strips are caused by the failure of the Scan Line Corrector (SLC) of Landsat 7.

Figure 5 .
Figure 5. Different types of surface materials in various locations in the river terminus on 11-Nov-2012.(A) Mixture of salt and sediment (salty surface); (B) Silt-rich surface; (C) Clay-rich surface; (D) Thick salt crust (pure salt).The locations for each surface type are indicated on Figure 4.

Figure 6 .
Figure 6.Grain size distribution of three types of sediment.Figures on the left-hand side (A1,B1,C1) are the grain size distributions of different surface materials (A: salty surface; B: silt-rich surface; C: clay-rich surface) while Figures on the right-hand side (A2,B2,C2) are their median grain size.

Figure 7 .
Figure 7. Statistics of the spectral library on 11 November 2012.A: Salty surface; B: Silt-rich surface; C: Clay-rich surface.Figures on the left-hand side (A1,B1,C1) are about the histogram of training samples (greater than 1000 pixels for each class) while on the right-hand side (A2,B2,C2) are the statistic boundary of each class.

Figure 8 .
Figure 8. Areal statistics of each class between 1985 and 2011.

Figure 10 .
Figure 10.The relationship between the annual precipitation and areas of different surface materials (A: Silt-rich surface; B: Salty surface; C: Clay-rich surface; D: River terminus) between 1985 and 2011.

Figure 11 .
Figure 11.Reactivation of a crevasse splay.The bold black lines indicate the trunk channel while the dashed lines show the expanded splays.See Figure 1 for the location of the splay.

Figure 12 .
Figure 12.The formation and expansion of a new crevasse splay.Figures on the left-hand side are classification maps, where the yellow area represents the crevasse splay surface, the green area the floodplain, and red area the salty surface.Figures on the right-hand side are schematic diagrams of the development of the new crevasse splay in a time series.See Figure 1 for the location of the crevasse splays.

Figure 13 .
Figure 13.Observed avulsion in the distal part of the river terminus.Figures on the left-hand side are classification maps, where yellow area represents silt-rich surface and green area is floodplain and red area indicates salty surface.Figures on the right-hand side are schematic diagrams of the development of an avulsion.A1 and B1 show the distal part of the river terminus.The red splay in A2 indicates where the avulsion happens.B2 shows a new river channel as a result of an avulsion.See Figure 1 for the location of the map.

Figure 14 .
Figure 14.Clay-rich surface with high soil moisture (A) and abundance of desiccation cracks (B).

Figure 15 .
Figure 15.Statistics of the tested area for a clay-rich surface.See the location in Figure 9.

Table 1 .
Days for which Landsat Climate Data Record (CDR) data are used.

Table 2 .
Region of interest (ROI) statistics of each class.

Table 3 .
Report of ROI separability for training pixels.

Table 4 .
Accuracy assessment of classification.

Table 5 .
Statistics of areas of each class.