Abstract

The advent of high-resolution, precise, back-pack portable terrestrial lidar scanners (TLS) provides a revolutionary new tool for obtaining quantitative, high-resolution (2-mm to 30-mm point spacing) measurements of landscape surface features. Moreover, data collected using these instruments allow observation of geomorphic processes in systems that can experience change on a daily basis. We have introduced TLS techniques in ongoing investigations of semiarid landscapes associated with weakly cemented sandstones along part of the Black Mesa escarpment of NE Arizona. Clay-cemented, Jurassic sandstones exposed along this escarpment are sensitive to moisture, and thus climate, via hydration-expansion weathering of interstitial clay. Sediment shed from weathered slopes has caused locally rapid valley floor aggradation and upper basin slope vertical denudation rates of 2–3 mm/yr over 10- to 100-yr timescales, as indicated by dendrochronology coupled with soil geomorphic analysis. These rates suggest rapid hillslope denudation rates. Employing the University of New Mexico Lidar Laboratory Optech Ilris 3D TLS, we are constructing a high-resolution model of two major basins along the escarpment. Focusing on a single, small (30 × 60 m) area of a mostly non-vegetated, steep slope (>35°), we demonstrate in this paper a method of comparative analysis of point-cloud data sets that can detect subcentimeter change resulting from a single season of monsoon precipitation along the escarpment. Using repeat scans can provide an empirical evaluation of single season erosion rates in the study site, and because our observations are geospatial in nature, we can also document the parts of the slopes that make the greatest contribution to local valley floor aggradation. In demonstrating the utility of this method, we expect that continued investigation of this site will provide insight to the key processes associated with soil-mantled versus bedrock-dominated slopes during modern escarpment retreat and hillslope modification, which, in turn, may further elucidate the impacts of Holocene climate change on this rapidly evolving landscape.

INTRODUCTION

Quantifying hillslope erosion rates and processes has the potential to provide insight on how hillslopes respond to both natural and human-induced perturbations. In turn this information can improve our understanding of how sediment is generated and transported from the hillslope to the active alluvial system (Yair, 1990; Gellis et al., 2004). Modeling approaches or the analysis of cosmogenic nuclides provide the most common way of determining geologic, millennial-scale, erosion rates (Brown et al., 1995; Schaller et al., 2001; Schaller et al., 2002), whereas modern erosion rates must be modeled (Govers et al., 2006) or measured empirically or by proxy. Traditional measuring methods employ a range of techniques, the most common of which are nail-washer lines or pins (Leopold et al., 1966), sediment monitoring within streams (Walling, 1991), or sediment trapping on the hillslope surface (e.g., Gerlach, 1967; Bryan, 1991; Larsen et al., 1995; Gellis, 1998; Gellis et al., 2004). Surveying (e.g., Govers and Poesen, 1988; Poesen et al., 1996) and digital survey techniques using laser theodolites (e.g., Belyaev et al., 2005) are also broadly employed to measure hillslope surfaces in situ. The application of reflectorless theodolites represents a significant improvement over previous methods because it allows for repeat measurements of the same system without disturbing the natural setting. It also provides for a multidimensional representation of hillslope change, whereas the other techniques are often only used to provide information about a single point (zero dimensional, 0D) or a vertical profile of a hillslope (a range of samples that sample an irregular line within a vertical 2D space). The greatest disadvantage to all of these measurement approaches is that they are time consuming and therefore limited in scope by the length of any reasonable investigation. Thus, it is difficult to collect a full spatial representation of hillslope shape from which process can be inferred.

To investigate the full scope of active hill-slope processes associated with a small drainage basin setting within a semiarid environment, we are incorporating the use of terrestrial lidar (light detection and ranging) scanners (TLS) technology to an ongoing investigation of punctuated hillslope modification and local basin floor aggradation (Tillery, 2003; Tillery et al., 2003; Burnett, 2004; McAuliffe et al., 2006). TLS approaches are capable of measuring 1000–20,000 points per second over a surface and depending on instrument, range, and conditions, can have subcentimeter precision and accuracy. Our principal purpose in this paper is to demonstrate this claim as well as to document a method for conducting a chronotopographic analysis of an active hillslope by comparing a time series of processed lidar point-cloud data.

The study area, which constitutes a much larger investigation, includes parts of the east- to southeast-facing escarpment of Black Mesa, which lies on the western flank of the Laramide-age Defiance Uplift (e.g., Kelley, 1958) in NE Arizona. The area scanned as part of this investigation includes Basins 1D and 1B (Fig. 2B) and focuses on an 800-m–long segment and 400-m–long segment, respectively, of the southeast-facing Black Mesa escarpment. We have also captured the valley floor to a distance of ∼500 m east of the escarpment base. Modern and millennial-scale degradation processes in this semiarid environment have produced a complex terrain and provide a unique opportunity to evaluate both hillslope erosional processes and localized aggradation of hillslope materials on the nearby valley floor. For this discussion, we are focusing on a small patch that forms a steep part of the Basin 1D escarpment (Fig. 2C). This particular location is representative of what is generally referred to as a detachment-limited slope formed upon soft, clay-cemented, Jurassic sandstones of the Morrison Formation and the uppermost San Rafael Group (Burnett, 2004).

Past studies have demonstrated that as drainage area increases, the availability of short-term sediment storage sites increases, which results in a net decrease in sediment yield as basin size increases (Schumm, 1977; Walling, 1983; Trimble, 1990; Gellis et al., 2004). A key question within the study area concerns the evolution of the slopes and the dynamic interplay between erosion, sediment flux, and short-term sediment storage. Specifically, which slopes (soil mantled versus bedrock) within the basin provide the most sediment, and how is this sediment moving into the fluvial system? The work to date in Blue Gap study area has shown that slope aspect and topoclimate play a key role in the magnitude of weathering (Burnett, 2004; Burnet et al., 2008), and Tillery (2003) established the relative importance of a thin soil mantle on some slopes in delivering sediment to the valley floor. TLS examination of both modern hillslope denudation and coupled valley floor aggradation may provide a new perspective on processes related to this episodic transport of sediments into active alluvial systems. More importantly, it may elucidate the geospatial distribution of key processes and how they relate to aspect and different slope types. A great deal of work has demonstrated the importance of diffusive processes where creep of material along the slope delivers sediment from the crest of the hill to the valley floor (Kirkby, 1971; Nash, 1984; McKean et al., 1993; Roering et al., 1999; Nash, 2005). There seems to be little consensus regarding processes on bedrock-dominated slopes, where diffusion is an untenable process. In these settings, sediment transport may be dominated by more advective processes limited by failure strength of the rock and the presence of water to help deliver the sediment to the local valley floor.

Working in a semiarid location has several advantages—because exposed bedrock is common, rock types are sufficiently susceptible to erosion (Bull, 1991; McFadden and McAuliffe, 1997), and storm-generated surface runoff can result in measurable erosion in a single event. Climate change on multiannual timescales has the potential to induce significant changes in both the driving and resisting forces affecting hillslope change, e.g., through drought-induced reductions in vegetation cover. Seasonal variations in precipitation provide an opportunity to capture annual-scale variation in the behavior of the system. All of these factors are important as we try to understand the contrasts between millennial-scale variation in erosion rates from the day-to-day processes that actually do the work of erosion, aggradation, and transport.

STUDY AREA

The study area spans part of the eastern escarpment of Black Mesa in NE Arizona and is positioned almost entirely within the Blue Gap quadrangle, ∼31 km west-southwest of Chinle, Arizona (Fig. 1). Here, a 150-m–high escarpment is formed on Jurassic-age strata principally of the Salt Wash Member of the Morrison Formation and the Recapture and Bluff Sandstone Members of the San Rafael Group (Tillery et al., 2003; McAuliffe et al., 2006). Basins draining eastward from this escarpment have been subdivided into several second- and/or third-order tributary subbasins (0–5) (Tillery, 2003; Tillery et al., 2003). Detailed petrographic analysis has demonstrated that at this location, the Jurassic strata largely consist of sandstone cemented nearly entirely by kaolinite and smectite (Tillery, 2003; Burnett, 2004; Burnett et al., 2008). Thus, the sandstones are highly susceptible to granular disintegration by hydration and expansion of clay cement, which locally forms a mantle of weathered sandy detritus up to decimeters in thickness over the bedrock. On cliffs and drier southerly aspects, however, hard, unweathered sandstone is often exposed (Burnett et al., 2008). The lower half of the escarpment is composed of middle Jurassic San Rafael Group sediments including the eolian Bluff Sandstone Member, which locally interfingers with fluvial mudstones and siltstones of the Recapture Member. A prominent unconformity (the J-5, c.f. Anderson and Lucas, 1997) separates the San Rafael Group from the late Jurassic Salt Wash Member of the Morrison Formation. Calcite-cemented sandstone and conglomeratic concretions are found throughout the Jurassic-age sediments, which are much more resistant than the surrounding sandstones and siltstones and typically form caps for hoodoos in the escarpment. We also note that in the northern part of the field area, the Morrison Formation is locally truncated by an erosional unconformity, and the escarpment rim is underlain by a white, medium- to coarse-grained sandstone with large kaolinitized feldspars that appear to be equivalent to the Jackpile sandstone of north-central New Mexico. The Jackpile Member thickness varies with strike and is unconformably overlain by well-cemented, fine-grained sandstones and carbonaceous mudstones of the Cretaceous Dakota Formation. The Dakota Sandstone locally contains abundant, thick lateral accretion units that grade laterally into a 3-m–thick coal. Over most of the field area, the Jackpile and Dakota units have been erosionally removed, leaving the upper sandstones of the Salt Wash Member capping the escarpment.

The segment of the escarpment east of the village of Blue Gap (Fig. 2) is positioned along the eastern flank of the Defiance Uplift, which is a Laramide-age monoclinal flexure that trends NNW-SSE (Fig. 1). Overall, the escarpment at this location trends ∼060°, which is subparallel to a regional fracture pattern measured by Burnett (2004); the pattern is nearly orthogonal to the trend of the monocline of the Defiance Uplift. Given the absence of any discernable dextral offset along the uplift (Lucas and Heckert, 2003) and the orthogonal trend of the joints with respect to the NNW strike of the monocline, it is likely that the observed regional joint pattern represents the orientation of syntectonic, mode I failure of strata along the broadly warped western flank of the uplift. This necessarily implies that these fractures are of comparable age to the uplift and therefore predate formation of the Black Mesa escarpment. More importantly, the escarpment is not the product of or influenced by younger, if not modern, tectonism that has affected other parts of the plateau. The comparable trends between these mode I fractures and the trend of the escarpment suggest that the regional jointing pattern is influencing the pattern of escarpment retreat. However, there are no reported indications that the study area has been influenced by any late Cenozoic tectonism (Burnett, 2004).

Hillslopes within the subbasins of the study area have been extensively studied in a series of investigations that have focused on coupled evolution between escarpment slopes and nearby broad valley floors (Tillery, 2003; Tillery et al., 2003, Burnett, 2004; McAuliffe et al., 2006; Burnett et al., 2008). Evidence for generation of surface on bedrock exposures is clear and often associated with stripping of the weathered sandy mantle and gullying on subjacent slopes. Although mass wasting of cliffs is probably relatively inactive in Holocene climates (e.g., Howard and Selby, 1994), some small recent rockfalls are nonetheless apparent. Focusing on Basins 0, 1D, and 1B, (Fig. 2B), McAuliffe et al. (2006) demonstrated a link between high, punctuated hillslope erosion rates and local valley floor aggradation. Fans and aprons of sandy alluvial sediment below these slopes grade into main-valley alluvial fills and show that while slow creep is likely occurring in the weathered mantle, episodic transport by surface runoff is the dominant erosional process.

Tree-ring analysis on the north-facing slope of Basin 1D and a south-facing slope in 1B has demonstrated a centennial-scale erosion rate in these basins of ∼2 mm/yr (r2 = 0.71; McAuliffe et al., 2006). These erosion rates were derived from measurements of tree root exposure and tree age. Piñon trees positioned on transport-limited hillslopes exhibit episodic, negative growth anomalies, suggesting that punctuated hillslope retreat has resulted in periodic root exposure, associated physiological stress, and diminished tree growth. Likewise, piñon and juniper trees within the nearby valley floor have been subjected to multiple aggradation events that resulted in a range of growth anomalies reflecting both burial and removal of the root system from the zone of effective soil moisture and sometimes local incision, which resulted in return of the root system to the zone of effective soil moisture and a concomitant increase in growth (McAuliffe et al., 2006). Given the clear relationship between pulses of hillslope erosion and local valley floor aggradation, it is apparent that the study area slopes are experiencing episodic pulses of rapid sediment transport that may peak well above the longer term average rate demonstrated by the treering data. We can test this hypothesis by using the TLS-derived, point-cloud data to measure seasonal and annual variation in vertical denudation and aggradation rates. Moreover, we can evaluate how these rates might vary as a function of aspect, slope, and depth of weathered and parent material.

METHODS

Ground-based (Nagihara et al., 2004; Bellian et al., 2005; Trinks et al., 2005) and airborne (Adams and Chandler, 2002; Gamba and Houshmand, 2002; Woolard and Colby, 2002; Charlton et al., 2003; Gibeaut, 2003; Haugerud et al., 2003; Storlazzi et al., 2003; White and Wang, 2003; Saye et al., 2005; Hilley and Arrowsmith, 2006; Young and Ashford, 2006; Heritage and Hetherington, 2007) digital survey techniques have been employed widely for the generation of high-precision representations of landforms and geologic outcrops. In most cases, an airborne instrument is used to generate high-resolution, digital elevation models (DEMs) that can, with the use of newer instruments, have submeter resolutions. DEMs based on data from TLS technology can measure 2–3 orders of magnitude finer resolution, whereas most scanners tend to have accuracy and precision of 2–30 mm. To date, such ultrahigh-resolution DEMs have been used principally in examining changes in dune position and shape (e.g., Nagihara et al., 2004) or to characterize features found within an outcrop (Bellian et al., 2005; Bellian et al., 2007). Our effort in this study is based on similar data sets, but rather than looking at the entire hillslope, we are focusing a subset of our data that represents a region that can be covered within a single scan. This simplified case is necessary because the method is extremely sensitive to errors introduced by merging scans during processing. In the end, this approach, when applied to a range of hillslope examples, should still give a reasonable representation of denudation (−δz) and aggradation (+δz) within this dynamic landscape.

We are employing the University of New Mexico (UNM) Lidar Laboratory's Ilris 3D, grayscale intensity, TLS system, wherein unlike newer scanners, color values of each point are not collected. With this system, intensity of laser reflectivity is collected and gained to a grayscale value between 50 and 255 for visual display. The raw-intensity data can also be extracted for processing purposes. The Ilris 3D is a lightweight scanner (∼13 kg) that provides subcentimeter scanning resolutions optimally from 3 to 300 m. Under optimal conditions (e.g., dry, clear, still air; 30%–40% reflective target surface, etc.), the scanner can collect subdecimeter data at ranges up to ∼800 to 1000 m. The unit employs a cylindrical beam design using a 905-nm wavelength laser source that yields modest spot sizes (13.7–182 mm, at ranges from 10 to 1000, respectively) under any light conditions. Minimum step size is 0.026 mrad, which yields a minimum spot spacing of 0.26–26.0 mm over the same functional range. Data are collected in a compressed format but can be exported as an American Standard Code for Information Interchange (ASCII) file, within which reflected points are reported in x, y, z coordinates, and a normalized intensity of reflection is reported as a 55–255 grayscale value. First-order variables that influence target reflectivity are distance and bulk reflectivity of the material. Variations in air density, humidity, and particulate load can induce error and/or reduce the intensity of the reflections, which implies that scanning should not be done during windy, moist, or dusty conditions.

Reflections are collected as a complete waveform, and the operating system can be programmed to collect the first or the secondary peaks. These reflections represent the primary and secondary strike on the target, respectively, where the secondary strike may represent an object directly behind the primary target (i.e., the ground behind a blade of grass). Given that it is an eye-safe laser, by design the beam is absorbed by water. The full kit including 4–6 batteries, cables, and necessary accessories (not including the tripod) weighs 25–30 kg and will function for ∼8 h of scanning activity. It can be deployed and used by a single user. The unit is self contained and operates in a range of environmental conditions. Optimal performance requires temperatures from 0 to 40 °C, and while the unit can function during light precipitation, beam absorption by water particles impacts reflectivity and creates some data loss. The beam is fully absorbed by accumulated snow.

Sources of error resulting from instrument design largely concern the configuration of proprietary electronic design. We are aware of some comparative studies to catalog such errors, which are currently under way, but to date no results have been made public. Optech cites a range of measurement error of no more than ±10 mm over a functional range of 3–500 m (B. Ysseldyk, 2005, personal commun.). Any other positional error is a function of spot size, where the spot represents a spatial Gaussian distribution, and the point measured is a function of the most reflective point within the spot with most of the signal returned within a radius that is roughly half of the spot size. In other words, in perfect conditions and spot diameter of 29 mm at 100 m, the maximum non-range error is ±14.5 mm. Such errors are typically overcome by a sampling step that is at least two-thirds of the spot size. At 100 m, the Ilris 3D can collect at a spot spacing of 2.6 mm with a spot size of 29 mm.

Field-Based Data Acquisition

The approach taken during data acquisition of Basins 1D and 1B required a minimum of 14 stations with 3–7 scans completed at each station. Each individual scan was calibrated for desired spacing 5–35 mm at the mean distance of the selected scan area. In some cases, high-resolution scans were desired over much smaller areas (e.g., Station 2, Fig. 2B). The desired spacing for these scans is generally 2–6 mm depending on distance and purpose of the scan. Collection of scan data in the field requires an organized approach—in our case, the data were collected from left to right and sequential scans were collected with a minimum of 2–5° of overlap. The overlap produced during scanning provides the basis for shape-recognition surface registration, a process that is handled in post-collection processing. Each station is referenced using global positioning system (GPS), and the orientation of the scanner is measured in the field using a Brunton compass. Both measurements aid in the merging process but are generally ineffective in porting the data accurately into a geospatial projection. This method has the added benefit of not forcing us to disturb the hillslope putting in survey targets. However, by measuring the trend and plunge of the edges of the base, we know the approximate orientations of both the y and x axes of the scan data, and they provide a reasonable reference when aligning and merging scans from different scan stations. High-resolution digital photographs were also collected of each scan. A single day of scanning results in 12–25 million data points that must be processed into a common reference frame, which can be geospatially referenced if targets of known position are captured in the scan. Alternatively, the data can be merged with an existing DEM. In this part of the study, we are only concerned with a ∼600 m2 area of a bedrock-dominated slope where we have conducted two repeat scans that were collected in June 2006 and October 2006 at a mean range of 35 m.

Sources of error during data collection begin with the instrument capabilities and are compounded by range distance, target obliquity and roughness, presence of water on the target, reflective intensity, atmospheric conditions, and ground motion. The latter can be avoided by setting up the scanner on a stable substrate. The other factors all serve to either reduce intensity of the return signal from the target or by beam deflection, which likewise impacts the return signal intensity. These errors are variable where range, water, and atmospheric conditions seem to have the most influence. Additional errors can be introduced by transient factors such as objects moving through the scan scene during data acquisition (e.g., animals, rain drops, etc.). Such types of errors are generally easy to identify and remove manually from the data set using any basic point-cloud processing software.

Post-Collection Processing

Once the data are collected, they must be processed to complete three basic tasks: (1) removal of measurement errors or bad points/scans; (2) merging of acceptable scans into a single reference frame; (3) data output into portable format or model generation. Data points representing unwanted target strikes can be easily trimmed from the data. However, depending on the objectives, it may be necessary to remove vegetation to map the underlying topography. Data representing vegetation can be much harder to remove and requires the development of filtering algorithms based on simplified DEM surfaces or assumptions regarding the thickness of the canopy above the topography. For this paper, we are presenting scans of essentially bare soil or rock surfaces.

Merging consists of two steps. First, scans collected from a single station must be merged into a common reference frame. Using I-Site Studio v. 3, scans are merged by generating a temporary surface for each scan and manually translating and or rotating the scan that is to be merged until it approximately overlaps the scan that is to remain fixed. In most cases, scans from a single station only need to be rotated. The same result can also be achieved by rotating the data to orient the x and y axes of the individual scans to the orientations measured in the field. Executing a best-fit algorithm for overlapping shapes completes the final merge between two scans. Errors in the final merge are reported as root mean square (RMS) values. For scans where the overlap spans area of cleanly exposed bedrock, RMS values range between ±3 and 25 mm. Because wind can cause plants to move, scans containing abundant vegetation often yield RMS values that range between ±15 and 45 mm. If errors exceed these values, merging is reprocessed using only clean bedrock exposures to reduce error, or the scan pair is flagged so that these scans are not used to merge scans from different stations.

The second step is to merge scan stations into a common reference frame. As before, a single station is identified as the fixed station, and the other stations are brought into the same project and visually aligned such that visible regions of overlap approximately coincide. Ideally, auto registration is applied only to the fixed scans of each station. During data collection, an effort is always made to include some ideal target in each station that can be used for registration in this step. Once the multiple stations have been merged, most of the multiple scans of the registration targets can be discarded. Typical RMS errors for this second step range between ±7 and 45 mm, and most full data sets contain multiple registration sites to evaluate the possibility of compounding merging errors. If desirable, the fully merged project can then be translated to a UTM (Universal Transverse Mercator) coordinate reference frame or rotated such that reference lines measured in the field are correctly oriented with respect to geographic north. The data can also be exported as ASCII format point data, which can be loaded into any surface modeling software capable of handling the size of the data files generated during scanning. It should be noted that point-cloud data require a voxel (a volume element within a gridded, three-dimensional space) or pseudovoxel rendering and processing environment. Unlike airborne lidar, which has a nadir (birds-eye) view of the ground and therefore generally only produces one z value for any given x-y position, a typical TLS-derived, point-cloud data file can contain multiple x-y pairs that have different z values. This is especially true when a vertical cliff face is captured within a scan. Within a pixel-based processing environment, all but one of the redundant x-y pairs is ignored or must be deleted. Most affordable, commercially and academically readily available software is only capable of fitting surface models to point-cloud data that contains unique X-Y data pairs (e.g., MATLAB, ArcGIS, Fledermaus Pro, RiverTools v. 3.01, Surfer v. 8, etc.). The work presented in this paper uses a combination of I-Site Studio v. 3 (beta release 8) and MATLAB to clean and merge scan data and to process the differences between two data sets that span the summer 2006, monsoon season, respectively.

To evaluate the differences between two point clouds, we choose to first demonstrate the method by comparing two scans from Station 2 along the northern end of the escarpment along Basin 1D (Fig. 2B). In June 2006, we captured an ultra-high resolution scan (2–5 mm) of an exposure of emergent and fully formed hoodoos positioned just below the escarpment rim (Fig. 3). Another scan (5- to 10-mm resolution) was collected in late September 2006. The difference in resolution between the two scans should only impact the lower limit of change detection between the two scans based on the resolution of the September scans, and assuming minimal instrument error, we are limited to only capturing 1–2 cm of change with this particular data set. Using the digital surface models (DSM) based only on data representing the hoodoo caps (the resistant heads of the hoodoos), we fully merged the June and September data sets into a single reference frame (Fig. 4). The merged data were then exported into ASCII format and loaded into MATLAB for comparative analysis.

Comparison of lidar data sets can be completed in two ways: (1) gridding individual scans and comparing values of z or intensity as modeled at regularly spaced grid nodes, or (2) comparing individual reference points to the nearest neighbors in the new data set. Since most programs cannot effectively generate true 3D DSMs, we are using a method defined in this study to compare the scans by a method we term Spot Comparison Analysis. A MATLAB script was prepared to calculate a range of statistical parameters regarding the difference between the June 2006 reference data set and the post-monsoon data set collected in October 2006 (Appendix 1). For each point in the reference data set, the script searches the new points data file to find any points within a user-defined search window defined in x-y space. Points found within the search window are then assigned to a temporary matrix, and for each new point graphic, the magnitude of the difference vector, δd, is calculated where the reference point graphic, using the following relation:  
formula

The smallest of the δd values (δdmin) derived by this calculation is a metric that represents the magnitude of a vector sometimes referred to as the Housdorff distance (Rote, 1991), which in this context is the shortest distance between the reference point and the nearest of the new points found in the search window. Although this is a useful metric for some analyses, we are ultimately more interested in the mean value of δd for each reference point.

The absolute difference in elevation (δz) is also calculated:  
formula

For any given point, if δz is less than zero, δd is assigned a negative value to represent denudation, which is an assumption that will not always be true when the data sets span a surface that has not undergone denudation or aggradation. The negative δz value could reflect instrument error rather than real change. When true denudation exceeds the range of error, this assumption should be valid.

For each reference point, the following parameters are recorded:  
formula
which are the positional (x,y,z), and intensity (i) values of reflection data for the reference point and the minimum δdmin, or shortest distance to any given point found in the search window around the point of reference, and the maximum δdmax or the longest distance for any given point found in the search window around the point of reference.
Other statistical moments of the data are also output. These include the mean of the difference magnitudes  
formula
where n is the number of points in the search window around the reference point, and δdi is the difference between the reference point and a single point in the search window. The mean difference in elevation for a given number of points (n) found in the search window;  
formula
where δzi is the difference between the reference point and a single point found in the search window. Finally, the standard deviation of the difference in elevation for the points found in the search window was calculated using:  
formula

The resulting matrix can then be used to construct files that can be imported back into I-Site Studio for visual rendering or, if desired, DSM generation. This requires that the parameter of interest be dynamically scaled to positive values ranging between 1.4135e104 and 6.5535e104 (following the scaling format for intensity of laser reflectivity in a Cyrex scanner format as defined in I-Site Studio) and exported as a space delimited, column formatted ASCII file containing what is equivalent to x, y, and z positional data associated with a fourth column containing the scaled Cyrex format intensity value. Mean denudation and aggradation rates can be calculated directly using functions built into MATLAB.

Sources of error can be generated in every step of the construction of merged data sets. Compounding errors are introduced when building a montage of scans from different positions. Careful processing and limiting the amount of merging is the only way to minimize such errors. At present, complete models of the entire escarpment cannot be completed because compounding errors of >1 m have been found over the 800-m length of the Basin 1D escarpment. Such an error is insignificant when compared to a 10-m DEM; however, it is unacceptably large for this analysis, which is the primary reason we are developing this approach and applying it to an area covered by a single scan.

RESULTS

To test our methods, we generated models of the entire escarpment of Basin 1D (Fig. 5; see Fig. 2B for location) and of a region of ∼600 m2 positioned just below the escarpment of Basin 1D (area 2 of Figs. 2B and 2C). Two topographic samples of Station 2 were collected ∼3 months apart (June 2006 and September 2006) and were successfully registered into a common reference frame (Fig. 4) with a best-surface fit of hoodoo caps with RMS error of 0.301e-3 (±3.01 mm). It is assumed that the calcite-cemented hoodoo caps are sufficiently resistant to not have changed during a single year. These two scan data sets span the southwestern USA 2006 summer monsoon season. A small part of the area 2 scans from both vintages were used to generate two DSM surfaces to graphically show the detectable variation between these two data sets (Fig. 6), representing changes that are visually obvious for this part of the Basin 1D slope (Fig. 3). The statistical results based on a 10-cm search window and the range of data spanned by these two scans is summarized in 01Table 1. Using MATLAB, the mean difference between each reference point and the new points was gained to a scaled intensity value based on the graphic range of –0.409 to 0.399 m, which represents a graphic span about the mean (graphic) of mean values (graphic) (−0.005 m), which were calculated using the following:  
formula
 
formula

The scaled intensity values and positional data were then imported back into I-Site Studio to display the range of lidar-detectable modification of the area 2 slope (Fig. 7). We hypothesize that the mean δz between the reference scan and the scan set taken three months later is representative of net slope denudation of 0.005 m and that the instrument is capable of mapping changes as small as 1–2 cm. Also, from the processed images, it is clear that the preliminary results of spot-comparison analysis provides a geospatial representation of change on the hillslope (Fig. 7) that agrees well with visual evidence of localized stripping of the slope during the 2006 summer monsoon (Fig. 3).

DISCUSSION

Although our results and interpretation are preliminary, it is clear that the spot comparison approach is a viable method for processing point-cloud data for change detection at resolutions approaching the measurement limitations of the instrument. Also, significant limitations still exist for using these data to generate large DEMs or DSMs to evaluate seasonal erosion and aggradation rates. Similarly, the lower detection limit needs to be further addressed, and it is likely that such limits need to be calibrated for a range of scanning conditions and individual instruments. The nature of this approach provides a spatial perspective on the day-to-day processes that influence hillslope modification in a semiarid environment, which is important in the context of how variation in the modern climate can influence erosional and aggradational processes.

The lower limit of detectable change using this approach is mostly a function of the instrument capabilities and environmental conditions during scanning. If the targets are not highly reflective, or it is windy, which is a product of density variations in the air and motion of targets on the ground, data quality is diminished and errors for individual measurements can increase beyond design tolerances (i.e., high winds or heat shimmer produce a much noisier data set that is not well suited for this type of analysis). Under ideal conditions, the instrument design limitations have the greatest impact. The two key properties that influence the detail that the instrument can capture are the size of the spot produced by the laser on the target and minimum spacing allowed by the instrument. Spot size is largely dictated by laser design, and it represents an area of reflection from which the peak in reflectivity is detected. Minimum spacing is a function of distance, and it is limited by the smallest step increment of the mirrors directing the laser source. The Ilris 3D used in this study uses a non-focused (cylindrical) beam from a fixed laser source where the laser spot size is given by:  
formula
and the minimum spacing is given by:  
formula
where D is spot diameter, R is the range to the target, and S is minimum spot center spacing (2003 Optech Design Specifications). From these relations, it is clear that data can be collected with overlapping spots (Fig. 8), which serves to provide multiple and overlapping data points that serve to constrain the shape of the target with sub-centimeter precision. This is not the case for coarsely spaced data where sampling distance is larger than the spot size.

Another factor that impacts the lower limit of detectible changes is the sampling resolution. Theoretical considerations dictate that the sampling resolution having a value d provides a series of points that sample the terrain but do not capture the full frequency of variation in that landscape. It has been shown that for a given value of d, the smallest frequency of variation that can be represented by the sample data set is 2d (Peucker and Chrisman, 1975). In other words, a 90-m DEM only portrays features in the landscape that have topographic variation expressed over a horizontal distance of 180 m. Finer spaced features simply are not sampled enough to appear accurately in the model representation. An example would be low-order stream channels spaced at less than 180 m apart. In our example, the sample resolution of the reference data set is ∼2–5 mm, which means the finest detail we can capture is 4- to 10-mm–scale features. Combined with an Ilris 3D ranging error of ±10 mm and the merging error of ±3 mm, it would seem that for these scans, the finest resolution of change we can hope to capture is on the order of 2–3 cm (under ideal environmental conditions). Inspection of the data seems to confirm that this may be the lower limit of detection for this particular data set, which has a mean range of 30 m with a corresponding spot size of ∼17 mm. By increasing the resolution of future scans, we may be able to push our detection limits to <2 cm, but this will require scanning at subcentimeter sample spacing. Moreover, it seems that comparing a high-resolution data set to a lower resolution data set results in detection limits dictated by the higher resolution scan.

On the scale of Basin 1D, using ultrahigh- resolution point-cloud data to generate a DSM of the escarpment is problematic. The scan merging process introduces 0.3- to 4.0-cm errors when merging the data, and these errors are propagated and cumulative as scans from a single station are merged. The misfit errors are larger still when combining multiple stations into a single reference frame. Seemingly well-merged data sets from two different stations have shown up to 20–40 cm of misfit for targets that are more than 600 m from a target used to merge the stations. Having multiple registration targets within a scan study area can minimize these types of errors. In our effort, we have initiated a reprocessing effort to develop a full, basin-wide data set of Basin 1D where such errors are minimized. Most of this work is still being done in I-Site Studio because the 3D complexity of these data goes well beyond the capacity of any pixel-based processing tools currently available.

At larger scales (smaller areas), merging chronotopographic scan data into a common reference frame is relatively straightforward, if some objects or targets within the scans have not moved. In this study, we merged the ultra-high-resolution scans from area 2 using part of a vertical face and several hoodoo heads to determine a best fit based on modeled surfaces generated from the two data sets (Fig. 4). We achieved a best-fit RMS error of 0.00301. This equates to a spherical error of ±3 mm when comparing reference data points to individual points collected in the follow-up scanning effort. The biggest limitation to this treatment of the data is that it enables the analysis of only a small piece of the entire slope. Although we are capturing a wide range in slope detail, we need to compile a series of these sites in different slope settings to gain a more representative sample of the variation of slope evolution along the Basin 1D escarpment.

Establishing a seasonal perspective on hill-slope process reveals that erosion and aggradation events are strongly coupled to existing soil conditions (mantled versus non-mantled slopes) and levels of precipitation necessary to transport hillslope materials. Over a three-month period, we appear to have documented a mean rate of change of −0.005 m (20 mm/yr), which, if accurate, would represent a short period of high-denudation rates that is much higher than the 10- to 100-yr timescale rate of 2 mm/yr established by McAuliffe et al. (2006). Such high rates are permissible with measurements over short timescales (Gardner et al., 1987), wherein it is possible to capture lower probability, larger magnitude but geomorphically significant events.

Based on regional records of rainfall for Arizona, we appear to have captured the impact of a major monsoon season following seven years of extended drought (Fig. 9), which may have been similar in some respects to the events that terminated a period of vegetation-reducing drought ending in 1906. The 1907 event was associated with a major hillslope degradation and local valley bottom aggradation event along the study escarpment (McAuliffe et al., 2006). It is possible that our TLS scans have in part captured the effects of unusually strong monsoon-related rainfall that caused sufficient runoff and locally some erosion along drainages on the slopes of Basin 1D, despite the presence of a piñon woodlands and at least some intercanopy vegetation. If this is the case, we may have captured the main erosional event within the past 4–8 yr. Averaging over this time period reduces the observed denudation rate to 2–5 mm/yr, an estimate that agrees well with the longer 10- to 100-yr timescale rate documented by McAuliffe et al. (2006). Continued monitoring of this site should help determine whether this is the case.

Finally, it is clear that this analysis provides an error range on the rate estimate of ±14 mm, which is also a function of temporal and spatial scale. Since the slope we analyzed has areas of aggradation and degradation, the error range should be large because the mean determination includes parts of the slope that have increased and decreased elevation. A low error range would imply that the slope only degrades or aggrades, or simply does not change at all. In this example, it is clear that there are large areas of aggradation. If we had not captured this particular event, then the denudation estimate would have been much larger. This highlights the importance of sampling scale and event frequency when attempting to directly measure changes in hill-slope morphology. Given the small study site with a maximum length sampled of ∼60–65 m, we are only adequately capturing hillslope erosion tracks that are spaced less than ∼30 m apart. Larger events such as the block fall (Figs. 3 and 7) may be exceedingly rare, or they may occur with significant frequency in the study area; further surveys along the entire escarpment may help to answer this question.

CONCLUSIONS

Chronotopographic analysis directly from point-cloud data appears to adequately detect changes on the hillslopes that exceed 1–3 cm along the Z axis, consistent with visual confirmation of stripping of hillslope mantle between June and September of 2006. This lower limit of detection is consistent with instrument and merging errors that are roughly an order of magnitude less than the observed change. The observed rate from point-cloud data is ∼20 mm/yr, which is higher than the 10- to 100-yr timescale rate determinations of ∼2 mm/yr (McAuliffe et al., 2006). This high rate may reflect the fact that we captured some unusually high rates in a single part of the slope, or we captured the effects of an exceptional range of precipitation events during the summer of 2006. In either case, continued scanning of the site and efforts to improve our processing approach should provide additional insight.

The results indicate localized erosional losses due to storm-generated runoff and erosion mostly by channelized flow, consistent with observational evidence for the dominance of these processes in Basin 1D and other parts of the escarpment. Erosion by creep may be locally significant on some of the more vegetated and mantled slopes, over longer timescales, but capturing this process at these time and spatial scales may be difficult. Application of TLS-derived, point-cloud data combined with construction of bedrock-anchored monuments and analysis of tree trunk shape should at least help us assess, if not quantify, the relative importance and spatial distribution of these processes throughout the mapped escarpment.

The authors gratefully acknowledge the comments provided by an anonymous reviewer and by Christopher Crosby, whose careful reviews served to greatly improve the initial manuscript. Monetary and material support for this study was provided by the UNM Lidar Laboratory, the UNM Center for Rapid Environmental Assessment and Terrain Evaluation (CREATE) (National Aeronautics and Space Administration Grant NN504AB25G to Louis Scuderi), and the UNM Department of Earth and Planetary Sciences. The lead author thanks J. Frechette, G. Weissman, and J. Galewsky for many fruitful discussions regarding the processing of TLS-derived, point-cloud data. This study was conducted on lands of the Navajo Nation, and we thank the Navajo Nation for the opportunity to conduct this research. Fieldwork was conducted under a permit from the Navajo Nation Minerals Department, and any persons wishing to conduct geological investigations on the Navajo Nation must first apply for, and receive, a permit from the Navajo Nation Minerals Department, P.O. Box 1910, Window Rock, Arizona 86515 (telephone: 928-871-6587). The authors also acknowledge Maptek and KRJA Systems for providing access to I-Site Studio v. 3 for processing of lidar point data.

Supplementary data