Abstract
Spatially explicit data for soil properties governing plant water availability are needed to understand mechanisms influencing plant species distributions and predict plant responses to changing climate. This is especially important for arid and semiarid regions. Spatial data representing surrogates for soil forming factors are becoming widely available (e.g., spectral and terrain layers). However, field-based training data remain a limiting factor, particularly across remote and extensive drylands. We present a method to map soils with Landsat ETM+ imagery and high-resolution (5 m) terrain (IFSAR) data based on statistical properties of the input data layers that do not rely on field training data. We then characterize soil classes mapped using this semiautomated technique. The method distinguished spectrally distinct soil classes that differed in subsurface rather than surface properties. Field evaluations of the soil classification in conjunction with analysis of long-term vegetation dynamics indicate the approach was successful in mapping areas with similar soil properties and ecological potential.
1. Introduction
Soil properties are critical for understanding patterns of vegetation community composition and primary productivity in arid and semiarid ecosystems globally [1, 2]. The relationship between incoming precipitation and water availability for plant growth in these water-limited systems is modified by vegetation, topography, and soil properties that affect surface redistribution, infiltration, and water retention [3–5]. Therefore, spatially explicit information for relevant soil properties is needed to understand the mechanisms governing plant species distributions and vegetation dynamics, particularly in efforts to forecast responses to changing climate. High-resolution imagery and terrain data sets are increasingly available globally, enhancing opportunities for digital soil mapping opportunities to yield accurate spatially explicit soil maps that capture many of these important hydrologic properties with reduced effort.
Although soil mapping is a long-standing science, broad scale mapping efforts often lack sufficient detail for understanding mechanisms regulating vegetation patterns and dynamics. This is particularly relevant for spatially extensive arid and semiarid regions (i.e., drylands). For example, in the United States, soil mapping efforts were initially focused on lands suitable for cultivated agriculture and fewer resources were allocated to mapping other areas, including drylands. Map units outside of valuable agricultural lands typically contain several different soils that can function quite differently from an ecological perspective. Therefore, in many drylands, updated soil maps with greater resolution and accuracy are needed to understand plant community patterns and dynamics. Digital, raster-based maps of soil properties are ideally suited for such analyses.
Enhanced availability of satellite imagery with increasing spectral, spatial, and temporal resolutions provide ample opportunities for predictive soil mapping at different levels of detail across a range of spatial extents [6–8]. McBratney et al. [9] proposed a framework for predictive digital soil mapping (DSM) that generalized Jenny’s [10] five soil forming factors (climate, organisms, relief, parent material, and time) to also consider spatial position and allow for interactions between soil forming factors to predict either spatially-explicit soil classes or discrete soil properties. The McBratney et al. [9] approach capitalizes on the availability of computing power and the ever-increasing wealth of remotely sensed data sources that serve as environmental covariates. The choice of satellite imagery for soil mapping should be based on cost and logistical (e.g., storage/computing) constraints in conjunction with the requisite detail of mapped result.
Land management and planning for spatially extensive drylands require prudent evaluation of the costs and benefits of available resources to meet the ever-increasing need for digital soil map layers. For example, in the USA, while detailed soil survey data are available in places, mapped soil units are often grouped together into general types based on the type and amount of potential vegetation and the site’s ability to respond to management activities based on soil-vegetation feedbacks and properties (i.e., NRCS Ecological Site, see [11]). For such applications, a lower spatial resolution product that can be easily obtained over large areas, such as with Landsat imagery, serves a broader purpose more effectively than more detailed, labor-intensive soil map products. Furthermore, archive and contemporary Landsat imagery provides an easily assessable source of data commensurate with landscape features that coincide with the rangeland monitoring and land-cover mapping [12, 13].
Spectral information regarding soil surface conditions and vegetation indices as surrogates for vegetation cover have been combined with high-resolution terrain models to improve DSM efforts [7, 14]. Boettinger et al. [6] effectively demonstrate the utility of remotely sensed imagery (i.e., Landsat) for characterizing soil surface features in drylands with modest vegetation cover. In another study focused on automated soil mapping with Landsat imagery and terrain layers, Saunders and Boettinger [15] combine unsupervised classification techniques with classification trees to evaluate utility of this combined approach to classify soils compared to a classification approach based on expert knowledge and field survey data.
Physical properties of the land surface relevant to soil forming factors are provided by satellite imagery and topographic features derived from digital elevation data [6, 16]. The availability of environmental covariates in digital format along with computing power and integration with local knowledge of change and degradation are key components to a worldwide effort to map soils for land management and carbon storage planning [17]. Despite advances in the availability of digital data and modeling algorithms (e.g., [14]), predictive digital mapping of soil properties at broad spatial scales is commonly hampered by a lack of supporting ancillary or training data [18]. This is especially so in spatially extensive and often remote dryland ecosystems. We present a statistically-based approach to derive spectral signatures for classifying soils without prior extensive field sampling and then characterize the soil types and properties of the mapped soil classes to facilitate interpretation.
We set out to isolate the effects of soil properties on dynamic changes in shrub and grass vegetation over 71 years in an arid grassland in southern New Mexico, USA. We assume soil properties delineated in our contemporary mapping effort represent those over the 71 year period. Requisite digital soils map data commensurate with patch-level dynamics were not readily available. The objectives of this study were to (1) delineate distinct soil classes within the study area using spatially explicit data layers for topographic and spectral features and (2) characterize derived soil classes to facilitate interpretation of observed vegetation dynamics.
2. Materials and Methods
2.1. Study Area
The study was conducted across a 150 ha landscape within the Chihuahuan Desert Rangeland Research Center in the northern Chihuahuan Desert of southern New Mexico, USA (Figure 1). The soil mapping endeavor was conducted as part of a larger study to discern the influence of soils in an analysis of long-term vegetation dynamics. The climate is characterized by a warm dry spring, hot wet summer, and cold dry winter. Long-term (1930 to 2008) average annual rainfall is 232.0 mm. Annual pan evaporation rates far exceed rainfall, with a measured annual average of 2,204.1 mm (1953 to 1979) [19]. The soil temperature regime for the area is Thermic, and the soil moisture regime is Ustic to Typic Aridic. The study area, with a long history of livestock grazing [20], occurs at 1,324 m elevation.
The area is dominated by sandy soils that are part of the broad alluvial plain of the ancestral Rio Grande River. Deposition of sediments by the Rio Grande on the alluvial plain ended approximately 1.6 million years ago [21], providing time for substantial pedogenic development including formation of argillic horizons and thick petrocalcic horizons. However, the area is now a mosaic of primarily sandy soil types due to postdepositional geologic and geomorphic processes, mainly tectonic uplift and reworking of the sediments by wind [22]. In areas receiving recent eolian deposition, surface textures are typically very coarse (>90% sand) and there is little clay or carbonate accumulation within the top meter. In eroded areas, the depth to the petrocalcic horizon can be relatively shallow (<50 cm) and the surface textures finer due to loss of surface soils.
The area was historically dominated by grasses but has transitioned to a shrub-dominated ecosystem [23, 24]. Shrub cover, dominated by honey mesquite (Prosopis glandulosa Torr.), increased from 1% in 1937 to ca.16% in 2008 [25]; this shift in mesquite dominance occurred amidst declines in grass cover, constituting black grama (Bouteloua eriopoda), tobosa (Pleuraphis mutica), and dropseed species (Sporobolus spp.) from ca. 19% in 1937 to 1% in 2008. Proliferation of shrubs and the decline in grasses were highly heterogeneous across space and time; interactions between soils and climate (i.e., drought) are expected to explain some of the observed heterogeneity.
2.2. Mapping Soils
A combination of unsupervised and supervised classification techniques were used to delineate distinct soil classes using spectral and topographic variables across the 150 ha study area. We chose to use Landsat 7 Enhanced Thematic Mapper (ETM+) imagery over that from other sensors for its accessibility (i.e., freely available) to land managers and soil scientists through the GLOVIS portal administered by the USA Geological Survey [26]. Spectral information was derived from a 24 February 2002 Landsat 7 ETM+ Level 1T image product that was precision and terrain corrected [27]. The winter image was selected from the Landsat archive to maximize the reflected radiant energy from the soil surface in winter, when vegetation cover is low and during a period of below-average rainfall. We assume that spectral properties reflected exposed soil surfaces based on a preponderance of evidence from long-term ecological research (LTER) site data sources collected at nearby study sites. Field-estimates of vegetation biomass collected 12–18 February 2002 indicate that production of annual and perennial grasses, forbs, and shrubs was near zero [28]. In addition, monthly observations of plant phenology corroborate the supposition that perennial grasses and shrubs were dormant when the image was acquired.
2.2.1. Input Layers
Topographic variables were derived from a digital terrain model (DTM) acquired with an airborne interfermetric synthetic aperture radar (IFSAR) sensor in 2006 (Intermap Technologies Inc., Englewood, Colo. USA). The native spatial resolution of the DTM was 5 m. Topographic derivatives evaluated were slope, slope shape (planar and profile curvature), and the Topographic Wetness Index (TWI). TWI provides a relative index of whether a point is in a landscape position likely to receive run-in water by taking the natural log of the specific catchment area divided by the local slope [29, 30]. Higher TWI values indicate a point with greater contributing area and/or lower slopes; lower values indicate points with less contributing area and/or higher slopes. Slope, slope shape, and TWI were calculated using Model Builder and Spatial Analyst in ArcGIS (Environmental Science Research Institute, Redlands, Calif, USA, v 9.3).
The pixel-based parametric classification scheme requires all input layers share a common spatial resolution. We chose to summarize the fine-scale terrain data within a neighborhood coinciding with the larger footprint of Landsat image pixels rather than to increase the redundancy in the spectral data by resampling to a smaller pixel size. This decision was based on the need to avoid limitations related to collinearity in the signature derivation process (see Section 2.2.2). We summarized topographic derivatives (using mean, standard deviation, minimum, and maximum values) within a focal pixel window corresponding to the pixel size of the ETM+ image.
The ETM+ image was radiometrically calibrated, and spectral reflectance values were adjusted to remove the effects of attenuation and scattering of photons due to atmospheric interference using the COST model [31]. The atmospherically corrected reflectance values for the six multispectral bands ranged in value from 0 to 1; spectral index calculations were based on atmospherically corrected reflectance values. Short-wave infrared ETM+ bands (5 and 7) are sensitive to surface soil moisture [32] and are commonly used in digital soils mapping [6]. Normalized difference ratios were calculated for established band combinations (e.g., gypsic soil index [16] using ETM+ short wave infrared bands 5 and 7 and iron oxide index [33] using bands ETM+ bands 3 (red) and 1 (blue)). Our combination of spectral bands and topographic derivatives as inputs for the classification process represents environmental covariates commonly used for digital soil mapping (e.g., [7]).
2.2.2. Assessment of Spectral and Topographic Input Layers
Reflectance values for image bands representing distinct, yet overlapping, portions of the electromagnetic spectrum are correlated. Collinearity in image band values can prohibit the derivation of spectrally distinct signatures [34]. Our initial attempts to employ a statistically driven assessment of signature separability with the entire set of 20 spectral and topographic image layers were unsuccessful because redundancy in the dataset prevented the creation of invertible covariance matrices in ERDAS Imagine. Therefore, we implemented a combination of methods to reduce the number of image bands to drive the unsupervised classification. In an effort to facilitate data interpretation, we did not use principle components analysis to reduce dimensionality of the dataset. Instead, we evaluated correlation matrices for image bands to identify redundant bands, quantified the optimal index factor (OIF) [35] to identify the subset () of spectral bands that maximize the variance within the Landsat ETM+ scene while minimizing duplication, and examined input image bands with the coregistered digital landform data layer. OIF is calculated for each 3-band combination as the sum of the standard deviations for the six Landsat image bands divided by the sum of the two pairwise correlations between the three candidate image bands [22].
2.2.3. Unsupervised Classification and Signature Derivation
The first step in the image classification process was to perform an unsupervised classification of the entire study area. Unsupervised classification is a computer-driven process that applies user-specified input parameters (such as convergence threshold and number of output classes) to initialize thematic classes and to identify clusters of pixels with similar spectral characteristics. The unsupervised classification algorithm supported by ERDAS Imagine (v.9.3, ERDAS, Inc., Norcross, Ga, USA) is the Iterative Self-Organizing Data Analysis Technique (ISODATA). Because the initial clusters in the algorithm are based, in part, by the number of output classes specified by the user, it is important to double the number of classes you expect to retrieve (based on available data).
Existing soils maps were used to determine the number of output classes to specify for the ISODATA algorithm. In our case, the existing 3rd-order soil survey conducted in 1980 (mapped at 1 : 48,000 scale, largely mapped as associations or complexes) delineated two soil map units across the study area: Berino-Bucklebar association (25% Typic Haploargid and 60% Typic Calciargid) and Wink-Harrisburg association (35% Typic Haplocalcid and 25% Typic Petrocalcid) [36]. A spatially distinct, ephemerally flooded 2.8 ha playa occurs within the study area; the playa is characterized by heavy clay soils and was delineated previously [22]. In the classification process, we manually assigned image pixels associated with the playa to a separate soil class. Combining information from the 1980 soil survey with a recent landform map of the area that included four landforms (Alluvial Plain Uplifted, Alluvial Plain Eroded, Alluvial Plain Wind Worked, and Playa; [22]), we estimated a maximum of five distinct soil types in the study area. Doubling the number of potential soil classes to allow for good separation in unsupervised classification process, we specified 10 output classes from the unsupervised classification of the entire study area.
The ISODATA algorithm was implemented with the following parameters. We specified 10 output thematic classes, selected a default convergence threshold of 0.950, initialized clusters using statistics encompassing 95% of the data along the principal axis, and allowed a maximum of 10 iterations. The ISODATA algorithm assigns pixels to clusters based on the minimum spectral distance from the cluster centroid. The centroid is recalculated after each class assignment with the convergence threshold serving as a measure of “classification completion.” This threshold represents the proportion of pixels that do not change classes from iteration to iteration. The ISODATA algorithm runs iteratively until either 95% of the pixels do not change clusters or the maximum number of iterations is completed. The 10 spectral signatures representing distinct soil classes of spectrally similar data were output for subsequent analysis and evaluation.
Transformed Divergence (TD) is a commonly used measure of signature separability in image classification [35]. The TD metric represents the spectral distance between two signatures based on the covariances between the signatures for a specific combination of spectral bands. Transformed Divergence is computed as follows:
where and = the two signature classes being compared, is the covariance matrix of signature , is the mean vector of signature , tr is the trace function, and is the transpose function [36, 37]. Mausel et al. [38] found that TD outperformed other statistical measures of separability in selecting the optimum subset of multispectral bands to distinguish between crop types in Hildago County, TX, USA. We used the TD measure to achieve two primary objectives: (1) determine the number statistically and spectrally separable soil classes occurring within the study area and (2) identify the subset of data layers that best discriminates between the distinct soil classes. The TD measure of signature separability ranges from 0 to 2,000 for a given class pair. We used a minimum value of 1,900 for all class pairs as the threshold to indicate spectrally distinct classes [35].
We used the TD measure in the ERDAS Signature Editor to determine the number of distinct soil classes based on spectral and topographic properties. We invoked the TD utility to calculate divergence values for all possible pairwise comparisons of 10 classes ( class-class comparisons). The TD utility outputs two matrices; the first matrix contains the class pairs and the second matrix contains the TD measure of separability for each associated class pair in the first matrix. For each iteration, we identified the class pair with the lowest TD value, merged the two specified classes in the ERDAS Signature Editor, deleted the original two signatures, and recomputed the TD metric. This process was performed iteratively until all class pairs had a TD value of 1,900 or higher, indicating spectrally distinct soil classes.
2.2.4. Supervised Classification
The second step in the classification process was to apply soil signatures derived from spectral and topographic input layers in the previous exercise back to atmospherically corrected ETM+ image (containing relevant spectral and topographic image bands) to perform the supervised classification. Supervised classification with parametric spectral signatures uses the statistical properties of the signatures to drive the class membership decision rule. We used the Maximum Likelihood decision rule, which does not employ a separate rule for overlap. This signifies that each pixel in the image will be applied to one of the spectral classes in the signature file. It is worth noting that supervised classification is traditionally used to classify multispectral imagery using field-collected or otherwise derived known vegetation or land surface data. In the absence of field-sampled training data for soils, we derived spectral signatures based solely on the Transformed Divergence metric to map distinct soils classes across the study area.
Small pixel clusters in the classified images may result from misclassification or aberrant reflectance data [39]. Further, while small patches of 1- and 2-pixel clumps could possibly represent soil inclusions, the number of small patches was low and would have complicated stratification of the site for field assessment. Therefore, we generalized the classified images slightly by removing 1- and 2-pixel clumps. These values were replaced with values of surrounding pixels.
2.3. Field Assessment of Mapped Soil Classes
To characterize soil classes mapped with this semiautomated classification procedure, we formulated a stratified random sampling strategy based on within-soil class topography. Depth to petrocalcic horizon is one of the key distinguishing soil properties in the study area. Much of the variability in the depth to petrocalcic in this area is due to reworking of the overlying sandy horizons by wind and water. To ensure our samples captured the representative range in deposition (likely deep) and erosional (likely shallow) areas we used stratified-random approach based on mean planar curvature within a moving window. We selected eight locations for field sampling within each of the three spectrally distinct soil classes (see Section 3.1, Digital Soil Mapping) for a total of 24 field samples.
For each sample location, we recorded horizon morphology and collected surface and subsurface samples for textural analysis. Soil was excavated using a hand auger to a depth of 100 cm or until a petrocalcic horizon was encountered. Horizon morphology information included descriptions of clay and carbonate accumulation in the top 100 cm. Each sample was classified to soil taxonomic Suborder following Soil Taxonomy [40]. At locations where a petrocalcic horizon occurred within 100 cm of the soil surface, we classified the soil to Great Group (e.g., Petrocalcid). Soil samples were collected from the surface (0–5 cm) and the horizon encountered with the maximum clay content (as estimated in the field by hand) for later particle size analysis in the lab [41]. For each soil class, average soil surface and subsurface texture and frequency of each taxonomic class were calculated using the probability of inclusion from the stratified sampling (PROC MEANS; SAS version 9.2, Cary, NC, USA). Analysis of variance was conducted on surface and subsurface texture to test if classes differed in surface or subsurface sand, silt, and clay (PROC GLM; SAS version 9.2).
3. Results
3.1. Digital Soil Mapping
3.1.1. Assessment of Spectral and Topographic Variables/Input Layers
The three-step process to assess input bands for the unsupervised classification led to the selection of four spectral and two topographic input bands. The first step was to calculate the optimal index factor (OIF) on the six spectral bands to identify the three-band combination that maximized within-ETM+ scene variance with the least redundancy. The 1, 5, 7 band combination yielded an OIF value of 0.701. The second step was to evaluate the correlation between spectral data (bands and indices; Table 1) and topographic derivatives (Table 2). The correlation matrix for the spectral data revealed that the gypsic index exhibited the least correlation with other spectral bands when compared with the other two normalized indices and that band 1 was not strongly correlated with bands 5 and 7 (Table 1). The short-wave infrared bands (5 and 7) were highly correlated (); the TD metric was used to determine which of the two short-wave infrared bands were more informative to characterizing the spectral signatures (i.e., Band 7, see Section 3.1.2). In general, topographic derivatives were less correlated with one another than the spectral bands (Table 2). In this case, visual inspection of candidate image bands (step 3) was most effective for identifying topographic derivatives for classifying soils. We evaluated maps for the 10 topographic derivatives with the existing (albeit general) 3rd-order soil survey [42] to determine that mean slope and maximum planar curvature most effectively represented relevant topographic features (i.e., relief and curvature). Selecting candidate image bands for an array of options to drive the image classification, while subjective, is a common procedure best accomplished in conjunction with available ancillary data sources [6, 15].
3.1.2. Spectral Signature Derivation
Starting with 10 signatures associated with classes from the unsupervised classification, we computed the Transformed Divergence (TD) metric eight times using all seven input bands (Band 1, Band 5, Band 7, B5−B7/B5+B7, B3−B1/ B3+B1, mean slope, and maximum planar curvature) to achieve a minimum TD value of at least 1,900 for each class pair (Table 3). The minimum TD value ranged from 1,511 when evaluating 45 class pairs between 10 clusters to 1,940 when evaluating three class pairs between three clusters. There was a clear break in minimum TD values between the six pairwise comparisons for four (1,715) and the three pairwise comparisons for three clusters (1,940) (Table 3). Using the three spectrally distinct signatures, we then calculated the TD metric using the best six (of seven) input bands to determine that ETM+ Band 5 was the least informative of the seven image bands. The subsequent supervised classification was then performed using six input bands (i.e., Band 1, Band 7, gypsic index (B5−B7/B5+B7), iron oxide index (B3−B1)/(B3+B1), mean slope, and maximum planar curvature).
Spectral response curves for the three spectrally distinct signatures revealed a low dynamic range in values for the spectral bands and the normalized ratios; in contrast, mean slope (summarized within windows) demonstrated the largest difference between the three classes (Figure 2(a)). Relativized topographic derivatives were used in the image classification; absolute values are presented (Figures 2(b) and 2(c)) to illustrate that mean percent slope for soil class 1 was 3.4 and was considerably higher than mean percent slope of 2.1 for soil class 2 (Figure 2(b)). There were no distinguishable differences in mean maximum planar curvature for the three soil classes (Figure 2(c)).
(a)
(b)
(c)
3.1.3. Supervised Classification
Pixels that coincided with the ephemerally flooded playa were manually assigned to a separate soil class 4, constituting 1.6% of study area. This resulted in 41.2% of the study area classified as soil class 1, 50.1% mapped as soil class 2, and 7.1% mapped as soil class 6; soil class 6 was interspersed between class 1 and class 2 pixels (Figure 3(a); Table 4). Field assessment of the mapped soil classes was based on 24 field sites.
(a)
(b)
3.2. Field Assessment of Soil Classes
Field diagnostics and lab analyses indicated that soil classes 1 and 2 exhibited different subsurface properties and that soil class 6 is a transitional soil more similar to class 2 than 1 (Table 4, Figure 3(a)). Of the eight sites surveyed in soil class 1, none exhibited a petrocalcic layer within the top 100 cm. Five of the sites (68.5% of the area), with very little clay or carbonate accumulation in the top 100 cm, were classified as Cambids; the remaining three sites had illuvated clays and were classified as Argids. Carbonate accumulation in the three Argid locations was very slight, with only a few carbonate filaments occurring at depth in most locations. In contrast, all eight field sites in soil class 2 had a petrocalcic horizon within the top 100 cm, distinguishing them as Petrocalcids. Average depth to a petrocalcic horizon in soil class 2 was fairly shallow (51.5 cm). All eight sites in soil class 6 exhibited some carbonate or clay accumulation; the majority of sites were Petrocalcids (63.6%) with the remaining sites classified as Argids (36.4%). Argids in soil class 2 had more visible carbonate accumulation with depth than did soil class 1; most soil class 1 sites had common carbonate filaments and nodules. No significant differences were detected in surface (0-5 cm) textures, although soil class 1 was slightly coarser (with more sand and less silt and clay) than soil classes 2 and 6 (Table 4). Subsurface textures of soil classes 1 and 2 were significantly different with class 2 having more clay (13.5 versus 9.2%) and silt (15.5 versus 11.7%) and less sand (71.0 versus 79.1%) than class 1. Soil class 6 had significantly less sand and more silt than 1 but did not differ in percent clay.
4. Discussion
4.1. Implications for Digital Soils Mapping
We demonstrate that combined use of unsupervised and supervised image classification methods with a semiautomated approach using the Transformed-Divergence metric to derive distinct signatures based on spectral and topographic features successfully delineated two soils classes with distinct subsurface soil properties (soils 1 and 2) and one transitional class (soil 6). We explore reasons why subsurface (rather than surface) soil properties were distinguishable with our mapping approach in subsequent sections.
4.1.1. Mapping in the Absence of Field Data
Derivation of spectrally distinct signatures for supervised classification using the Transformed Divergence metric circumvented the need for guesswork in the soil mapping process. Yet, only with a field sampling effort to characterize the mapped soil classes across the study area were we able to interpret the final product. This approach not only yielded an interpretable map of soil classes, but the delineated unique classes also provide an effective means to stratify the landscape for other types of interpolation (e.g., kriging). Such a stratification approach has some advantages over other sampling approaches designed to mimic the natural distribution of soils (e.g., [43]) in that it is effective in capturing rare but complex soil units (e.g., soil class 6), where a higher number of samples was needed for optimal characterization.
The combined use of unsupervised and supervised classifications in digital soil mapping is not necessarily novel. Saunders and Boettinger [15] combine unsupervised classification techniques with classification trees to compare a semiautomated image classification approach based on expert knowledge and field survey data. The novelty of our approach was in using a statistical metric to delineate distinct soil classes in the absence of field data or expert knowledge to train the supervised classification and characterizing the mapped soil classes as a means to assess the performance of this approach.
Field evaluation of mapped soil classes in this study indicates this approach can accurately map the distribution of soil types within a landscape without field samples to drive the image classification. We were surprised by the relatively small contribution of image bands in the visible spectrum and the TWI to distinguishing soil classes. We expected that shallow calcic and petrocalcic horizons would display different levels of brightness near the soil surface. In addition, selecting a Landsat image acquired in winter during a period of below-average rainfall to maximize reflected soils properties may have dampened the benefit of short-wave infrared bands (Band 7 and the normalized ratio of Bands 5 and 7) which are generally sensitive to near surface soil moisture. Although Mean TWI (calculated within window) did reflect a higher mean index value for the ephemerally flooded playa (Figure 4), it did not relate well to ancillary soil layers in the data screening process. It is possible that the ideal catchment area size for TWI calculations in this system is larger (or smaller). Future work will examine TWI at a range of spatial scales (i.e., window sizes). In contrast, the slope layer was an effective covariate in distinguishing soil classes. This is due to the improved accuracy and detail (5 m) afforded by the IFSAR sensor and the prominence of landform in soil formation at this site (see Section 4.1.2).
The increased resolution of terrain derivatives contributed to the effectiveness of this mapping effort. We contend that the enhanced accuracy of the IFSAR digital terrain model effectively captured slope and slope shape and that these properties were maintained in the summarized depiction of the derived layers. If our analysis had been based on a 30 m resolution slope layer derived from a widely available digital elevation model, we feel that subtleties of slope shape contributing to overland water flow would have been missed. Previous attempts to incorporate slope derived from the previous standard 30 m elevation product to map vegetation were unsuccessful [34]. This was particularly important for distinguishing soil class 6 which was juxtaposed between soil classes 1 and 2 (Figure 3) and may not have been particularly discernable using a native 30 m DEM product.
Spectral response of bare soils has been shown to demonstrate a linear relationship between reflectance in red and near-infrared image bands [44], where soil moisture and surface roughness are the primary factors which determine where upon the line individual soils occur [45]. Vegetation within a given image pixel will contaminate the pure soil response, and this principle has been applied to depict vegetation contributions to reflectance [46, 47]. We demonstrate that reflectance values in the 23 Feb 2002 ETM+ image in this analysis were not detectably influenced by vegetation due to general conformity of all pixel values to the soil line (Figure 5). As noted earlier, soil moisture was very low for all soils during this time period, thereby circumventing its utility in distinguishing between soil types.
4.1.2. Mapped Soils in the Context of Landform
The success in mapping these subsurface properties was likely due to the correlation between slope, slope shape, and geomorphic processes that affect soil formation. Soil class 1 partially corresponds with the Alluvial Plain Uplift geomorphic area mapped by Monger et al. [22] (Figure 3(b)). The study area slopes gently to the north and is just below the ridgeline of this gentle uplift. The predominant wind direction in the study area is from the southwest to northeast [19]. Hence, soil class 1 is immediately on the leeward side of this ridge and likely has been accumulating eolian sediments since the uplift occurred, resulting in the younger soils with less clay and carbonate accumulation in the soil class. In contrast, soil class 2 is on a flatter, more stable portion of the landscape, where deposition of eolian sands is less pronounced and the petrocalcic horizon is within 100 cm of the soil surface.
Slope shape (maximum planar curvature), which in this system is likely an indicator of reworking of sediments by wind, was similar for soil classes 1 and 2. Although broad-scale deposition was not likely occurring across soil class 2, reworking of surface sediments by wind is common in this geomorphic surface [22, 48]. Soil class 6 is intermediate to classes 1 and 2 both in terms of soil properties and topography (Table 4, Figures 2(b) and 2(c)). Soil class 6 slope is intermediate to the flat class 2 and sloping class 1 and slope shape of class 6 was smoother (i.e., fewer undulations) than that of classes 1 and 2. This class primarily occurs at the toe-slope of soil class 1, and in the transition from soil class 2 to the playa (soil class 4). The shallow depth to petrocalcic and commonly observed surface carbonate fragments in the field indicate that soil class 6 represents areas dominated by erosion of surface horizons, due either to wind or water.
4.2. Enhanced Interpretation of Long-Term Vegetation Dynamics
We set out to generate a spatially explicit depiction of ecologically relevant soil properties to evaluate the influence of the topoedaphic template on our long-term analysis of vegetation dynamics. An understanding of the interactions between soils, vegetation, and landscape treatments is strongly desired by land managers and decision makers. The digital soil map produced by this effort captured important differences in soil properties that help to explain the divergent shrub and grass vegetation dynamics over the study area. Our analysis of changes in shrub and grass cover from 1937 to 2008 revealed that grasses were distributed over both dominant sandy soils in 1937, but grass loss and shrub proliferation diverged following the 1950s drought [49]. Soil class 1 was characterized by low water holding capacity in the top 100 cm with little clay or carbonate accumulation with depth. Soils with these properties have been correlated with low grass resilience to drought due to inability to retain water during dry periods [3, 50]. Such low near-surface water holding capacity is expected to benefit mesquite (Prosopis glandulosa, the dominant woody invader in this system) since this species develops extensive, deep rooting systems soon after plants are established [51, 52]. In contrast, soil class 2 captured soils with slightly higher clay amounts and much greater near-surface carbonate accumulation (in the form of petrocalcic horizons) than soil class 1. Such shallow petrocalcic soils are known to promote resilience of grasses during drought due to their ability to retain water at plant available tensions for nearly a year following rain events [3, 40]. This high near-surface water holding capacity would likely lessen the competitive advantage of deep-rooted mesquite plants relative to shallow-rooted grasses by retaining a larger proportion of infiltrated water where it is available to both shallow- and deep-rooted species. Soil class 1 does also have a slightly steeper average slope (3.4%) and soil class 1 (0.9%, Figure 2). Soils with steeper slopes are likely to have a greater fraction of precipitation be lost to runoff than soils with shallower slopes. However, coarse surface textures (Table 4), modest differences in slope, and the importance of petrocalcic horizons for water dynamics [3] and grass persistence during drought [50] indicate that the subsurface soil properties were likely more important factor than the topography governing vegetation dynamics.
We demonstrate that it is possible to extract ecologically meaningful information about soil properties from a remotely sensed perspective. Extensive field sampling and knowledge of ground conditions is not required a priori, but the method does require post-hoc allocation of effort to characterize soil class maps generated. This study used established image processing techniques with a semiautomated method to derive soil class signatures to ultimately distinguish differences in site potential (i.e., shrub dominated versus shrub-invaded grassland; [49]). Field evaluations of the resulting soil classification and analysis of long-term vegetation dynamics among soil classes indicate the approach was successful in mapping areas with similar soil properties and ecological potential. As such, the method provides a basis for mapping soil classes across landscapes or for effectively stratifying sites to make the most efficient use of human resources.
Acknowledgments
This research was funded by the USA Department of Agriculture, Agriculture Research Service, and the National Science Foundation Long-Term Ecological Research Program, Jornada Basin LTER V: Landscape Linkages in Arid and Semiarid Ecosystems. M. Mattocks made substantial contributions to the field assessments and lab analyses to characterize soil properties. Derek Bailey facilitated access to the Chihuahuan Desert Rangeland Research Center for field analysis and sampling.