## Abstract

The need to accurately document the spatiotemporal distribution of earthquake-generated strong ground motions is essential for evaluating the seismic vulnerability of sites of critical infrastructure. Understanding the threshold for maximum earthquake-induced ground motions at such sites provides valuable information to seismologists, earthquake engineers, local agencies, and policymakers when determining ground motion hazards of seismically sensitive infrastructures. In this context, fragile geologic features such as precariously balanced rocks (PBRs) serve as negative evidence for earthquake-induced ground motions and provide important physical constraints on the upper limits of ground motions. The three-dimensional (3D) shape of a PBR is a critical factor in determining its static stability and thus susceptibility to toppling during strong ground shaking events. Furthermore, the geomorphic settings of PBRs provide important controls on PBR exhumation histories that are interpreted from surface exposure dating methods. In this paper, we present PBRslenderness, a MATLAB-based program that evaluates the two-dimensional (2D) static stabilities of PBRs from unconstrained digital photographs. The program’s graphical user interface allows users to interactively digitize a PBR and calculates the 2D geometric parameters that define its static stability. A reproducibility study showed that our 2D calculations compare well against their counterparts that were computed in 3D (R^{2} = 0.77–0.98 for 22 samples). A sensitivity study for single-user and multiuser digitization routines further confirmed the reproducibility of PBRslenderness estimates (coefficients of variation *c*_{v} = 4.3%–6.5% for 100 runs; R^{2} = 0.87–0.99 for 20 PBRs). We used PBRslenderness to analyze 261 PBRs in a low-seismicity setting to investigate the local geomorphic controls on PBR stability and preservation. PBRslenderness showed that a PBR’s shape strongly controls its static stability and that there is no relationship between a PBR’s stability and its geomorphic location in a drainage basin. However, the geomorphic settings of PBRs control their preservation potential by restricting their formation to hillslope gradients <40° and the upper reaches of drainage basins. Such examples of our program’s utility have led to its use in archival efforts of PBRs in southern California and Nevada, USA.

## INTRODUCTION

Strong ground motions generated by large earthquakes pose significant hazards to human life, infrastructure, and the socioeconomic state of seismically active regions. Efforts to reduce seismic risk require knowledge of ground motions that are produced during seismic events. Models that describe these ground motions use prior knowledge of the distribution of seismic sources in space and time in order to provide forecasts of the spatiotemporal extent and magnitude of ground shaking. In some cases, direct evidence of shallow, surface-rupturing seismic events is manifested in the geologic and geomorphic records via fault scarps, offset stream channels, and coseismic mass wasting (Arrowsmith and Rhodes, 1994; Arrowsmith et al., 1998; Oldow and Singleton, 2008; Akçiz et al., 2010; Toké et al., 2011; Haddad et al., 2012). However, these records provide limited information about the spatiotemporal extent and magnitude of ground shaking that play key roles in building ground motion models for seismic hazard analyses and earthquake engineering studies. Due to the short record of seismic instrumentation and significant spatial variability in ground motions, empirical data for ground shaking are still insufficient. The high variability in peak ground accelerations has to be acknowledged in seismic hazard assessments and often results in higher estimates of seismic hazards (Bommer and Abrahamson, 2006). Moreover, the physical processes involved in generating high peak ground accelerations combine source-, path-, and site-related effects, making ground motion prediction more difficult (Strasser and Bommer, 2009). This high variability affects all parameters of ground motion and can change between earthquakes with similar magnitude, depth, mechanism, or even in the same event, when measured at different seismic stations (Bommer and Abrahamson, 2006). For example, the February 2011 Christchurch, New Zealand, earthquake (local magnitude, *M*_{L} 6.3) generated some of the highest recorded peak ground accelerations, which caused widespread liquefaction and, due to the close proximity of the epicenter to the town, caused severe damage and casualties (Bradley and Cubrinovsko, 2011). The observed peak ground accelerations exceeded previous records and exceeded those generated by the 2010 Darfield earthquake, which had an even higher magnitude (*M*_{L} 7.1), but was located ∼30 km west of Christchurch (Bradley and Cubrinovsko, 2011).

Slender human-made structures and fragile geologic features allow for the documentation of past earthquake-induced ground motions and have been of interest to geoscientists for more than a century (Mallet, 1862; Milne, 1881; Milne and Omori, 1893; Kirkpatrick, 1927; Housner, 1963). Such efforts were motivated by the prospect of inferring the intensity of earthquake ground shaking from overturned columns and water towers. The first attempt to use a precariously balanced rock (PBR; Fig. 1) to constrain estimates of peak ground motions was by Coombs et al. (1976), who used the Omak Rock PBR in northern Washington State to estimate peak ground accelerations 0.1–0.2 *g* generated by the moment magnitude, *M*_{w} 6.5–7 1872 Pacific Northwest earthquake. Advancement in seismic instrumentation and historical earthquake records enabled the documentation of the spatiotemporal distribution of seismicity over decadal to centennial times scales. Subsequent interest in the utility of PBRs was renewed in the early 1990s to extend ground motion constraints to millennial time scales for critical infrastructures such as the Yucca Mountain nuclear waste repository (e.g., Brune, 1992, 1993, 1994; Weichert, 1994). This was based on the assumption that PBRs have been in their present state for thousands of years, and that by linking their static stabilities with ground motion characteristics they may be used to constrain peak ground motions on the scale of multiple earthquake cycles. Later applications of PBRs were extended toward physically constraining probabilistic seismic hazard analyses for southern California by Brune (1996), who found that PBRs were located near large faults with modeled peak ground accelerations of >1.5 *g* in 50 yr with a 2% exceedance probability (∼2475 yr recurrence interval). These modeled peak ground accelerations were well over those needed to overturn the surveyed PBRs (0.2–0.3 *g*). Using the assumption that the PBRs have been precarious for thousands of years, Brune (1996) speculated that they must have survived multiple episodes of strong ground motions and that the existing probabilistic seismic hazard analyses for southern California overestimate the hazard posed by earthquake-induced ground shaking. A fundamental parameter in this assumption is the time since the contact between the PBR and its pedestal, and therefore the PBR’s balancing time, has been exposed (Fig. 1). Determining this time is vital to the utility of PBRs in physically constraining ground motion exceedance probabilities that refine probabilistic seismic hazard analyses (Purvance, 2005; O’Connell et al., 2007). As a result, applications of varnish microlamination and cosmogenic radionuclide surface exposure dating techniques in PBR research became widely used to determine the geomorphic history of PBRs (e.g., Bell et al., 1998; Stirling et al., 2002; Stirling and Anooshehpoor, 2006; Balco et al., 2011). Despite the seismological confidence in the utility of PBRs as natural seismometers, the processes that form them and control their static stabilities during their geomorphic histories are insufficiently understood. Interpretations of PBR exhumation histories from surface exposure ages are reconstructed without accounting for the PBR’s overall geomorphic situation in the landscape. Because the fundamental mechanistic hillslope transport and soil production laws are largely controlled by the local hillslope gradient and upslope contributing area in a landscape (Gilbert, 1877; Penck, 1953; Schumm, 1967; Kirkby, 1971), the geomorphic setting of a PBR may control its exhumation rate, residence time, and static stability. These lines of geomorphic reasoning have not been applied in PBR research; there is therefore a need to investigate the possible effects of local geomorphic settings of PBRs on their static stabilities.

In this paper, we present a simple and efficient method for determining the two-dimensional (2D) geometric and geomorphic parameters that control the static stabilities of PBRs in a low-seismicity region. We first describe the mathematical background for computing PBR stability parameters in 2D, then present the graphical user interface (GUI) implementation of a software and the workflow involved in estimating PBR stability parameters from unconstrained digital photographs. Because variations in the shape and rocking points of a PBR control its stability, we present reproducibility tests of our program’s stability estimates using digitization routines performed by multiple GUI users. We then compare our 2D estimates with those made by a three-dimensional (3D) parameter calculator that uses photogrammetrically generated 3D surface models of PBRs.

## METHODS

### PBR Geometric Parameters

An important parameter that is used in the equations of motion that describe earthquake-induced quasi-static accelerations of PBRs is slenderness (Fig. 2A; Housner, 1963; Shi et al., 1996). Slenderness is defined by the angle (α_{i}) between the PBR’s rocking arms (*R*_{i}) and a vertical reference passing through the center of mass; α_{i} is used to model the rocking response of a rigid body due to accelerating motion at its base as follows (Shi et al., 1996):

where *I*_{i} is the moment of inertia about the PBR rocking points, t is time, and *A*_{x}(*t*) and *A*_{y}(*t*) are the horizontal and vertical accelerations at the PBR base, respectively. The value of α_{i} and its corresponding *R*_{i} depend on the azimuthal orientation of an imaginary vertical plane that dissects the PBR through its center of mass and rocking points (Fig. 2B). The smallest value of α_{i} (α_{min}) defines the PBR’s potential for overturning and thus its fragility (Purvance, 2005; Purvance et al., 2008a, 2008b). The challenge is to find the vertical plane that contains the PBR’s absolute smallest α_{i} value (α_{min}) and corresponding *R*_{i}.

Existing methods for determining α_{i} and *R*_{i} use force meters, plumb bobs, measuring tape, and a trained eye (e.g., Anooshehpoor et al., 2004). Recent 3D estimates of α_{i} and *R*_{i} were determined using photogrammetric and terrestrial laser scanning techniques to capture the 3D geometry of PBRs at decimeter and millimeter scales, respectively (e.g., Anooshehpoor et al., 2004; Haddad et al., 2012). These methods produce detailed 3D surface models of PBRs on which 3D shaking tests are performed. However, specialized equipment and expensive post-processing software are needed to build 3D models of PBRs that may not be necessary to capture the important parameters that define PBR stabilities, thus making these methods impractical for the reconnaissance of large (i.e., several hundred) PBR populations to identify PBRs that warrant more detailed fragility analyses.

### 2D Parameter Estimation

_{i}and

*R*

_{i}can be made along a variety of photograph azimuths (Purvance, 2005). Determining α

_{i}and

*R*

_{i}of a PBR in two dimensions from a digital photograph requires calculating the PBR’s 2D center of mass (

*C*

_{x},

*C*

_{y}) given

*n*vertices (

*x*

_{0},

*y*

_{0}), (

*x*

_{1},

*y*

_{1})…(

*x*

_{n–1},

*y*

_{n–1}) and the PBR’s 2D area (

*A*) as follows:

We developed a MATLAB-based program, PBRslenderness, which implements the preceding approach into a GUI (Fig. 3). PBRslenderness allows users to interactively digitize the outline of the PBR, its rocking points, length scale, and vertical reference from an unconstrained digital photograph. PBRslenderness then computes the PBR’s 2D center of mass, which is used in combination with the user-defined rocking points, vertical reference, and length scale to determine α_{i} and *R*_{i} values for the PBR. Where multiple photographs are to be analyzed, PBRslenderness allows users to append results of a current run to previous runs, thereby facilitating the analysis of large libraries of PBRs. (PBRslenderness and a comprehensive video tutorial can be found in the Supplemental File^{1}.)

### Study Area and Sampling Protocol

Selection of an appropriate study area to investigate the geomorphic controls on PBR static stabilities was dictated by the seismotectonic and geomorphic settings of the PBRs. Our study area is located in the Granite Dells precarious rock zone of the low-seismicity Arizona transition zone between the southern Basin and Range province and the Colorado Plateau (Fig. 4). The PBRs are developed in the Proterozoic Dells Granite, which forms a prominent pediment that is dissected by angular joint-controlled drainage networks. With the exception of localized compositional variations, the Dells Granite occurs as a massive, medium- to coarse-grained granite of locally porphyritic texture with grain diameters between 3 and 8 mm. We chose to study PBRs in a low-seismicity region in order to reduce possible seismic contamination of the static stabilities of the PBRs since their formation. Because a population of PBRs may contain a geomorphic life cycle and a seismic life cycle, for any one PBR the geomorphic life cycle begins when the PBR pedestal is exhumed and ends when the PBR is no longer precarious (i.e., toppled). Each PBR in a population may thus be at a different stage of its geomorphic life cycle. Conversely, the seismic life cycle of a PBR begins when the last strong ground motion episode took place before the PBR became precarious, and ends when the next strong ground motion episode occurs (whether or not the PBR topples). Combined, the two life cycles form a complex dynamic system that increases the complexity of the seismogeomorphic life cycle of a PBR population (O’Connell et al., 2007). Therefore, we chose PBRs that formed in the low-seismicity Granite Dells precarious rock zone to eliminate possible seismic or tectonic geomorphic contamination from the seismogeomorphic life cycle of the PBRs and to understand the underlying geomorphic processes that control the static stabilities of PBRs in their simplest setting.

We searched for PBRs by traversing along an ∼300-m wide segmented east-west transect (Fig. 4). Only those PBRs that were geomorphically developed in situ on bedrock pedestals were sampled. This assessment was made by ensuring that no discrepancies existed between the orientations of the sides and surrounding vertical joints of the PBRs, and by searching for fresh or broken surfaces that indicated damage or disturbance to the PBRs. In addition, only PBRs that were fully detached from their pedestals were surveyed in order to satisfy the assumption of a block that is free to rotate about its rocking points (Fig. 2A; Housner, 1963; Shi et al., 1996). We made this assessment by visually inspecting the PBR-pedestal contact for a throughgoing horizontal joint (Fig. 1). The location of each PBR was recorded using a Wide Area Augmentation System–enabled hand-held global positioning system (GPS) receiver (±2 m horizontal accuracy).

Our sampling protocol for each PBR was as follows: the PBR was selected according to the above-mentioned criteria and its location was recorded by the GPS receiver. The PBR’s maximum circumference and height were then measured using a graduated rope and measuring tape, respectively. A freely suspended plumb bob and scale of known length were placed next to the PBR. Digital photographs were then taken of the PBR, plumb bob, and scale from different azimuths using a Sony DSC-S600 digital camera. For each photograph, the scale was rotated such that it was perpendicular to the photograph’s azimuth (i.e., it was in the plane of the photograph). Each photograph was then processed using PBRslenderness to determined α_{i} and *R*_{i} for each PBR.

### PBR Geomorphic Parameters

^{−2}) point clouds in a global reference frame. These points were then used in an inverse distance weighting (IDW) local binning algorithm to generate submeter-resolution digital elevation models (DEMs) as follows (El-Sheimy et al., 2005):

*Z*

_{IDW}is the interpolated distance-weighted elevation computed for each grid node,

*n*is the total number of grid nodes,

*l*is the grid node index, and

*r*is the radius of the node-centered computational bin. Two geomorphic parameters were computed from these DEMs, and included local hillslope gradient and upslope contributing area. Local hillslope gradients were computed by fitting a plane to a 3 × 3 elevation grid window around a central cell node whose local slope was being computed. The maximum slope value of this plane was then assigned to the central cell and the computation window was moved to the adjacent central cell (DeMers, 2002). This process was performed for all raster cells in the DEM and produced gridded slope maps at the resolution of the DEM. Stream channels were defined as grid cells with an upslope contributing area >100 m

^{2}computed using the D∞ flow routing algorithm of Tarboton (1997). The D∞ algorithm routes flow into an infinite number of directions by using a triangular-faceted computation cell. The computation cell is divided into eight triangular facets and the flow direction for each facet is computed. The algorithm then searches for the facet with the steepest slope and assigns the flow direction. This in essence subdivides the computation cell such that a more realistic flow direction raster is produced and from which upslope contributing areas and channel networks may be extracted (Tarboton, 1997). Each PBR was assigned its corresponding slope and area values, which were extracted from the local hillslope gradient and upslope contributing area rasters, respectively.

## RESULTS AND DISCUSSION

### Parameter Reproducibility: Single-User and Multiuser Comparisons

A possible source of uncertainty in the α_{i} and *R*_{i} values estimated by PBRslenderness is the reproducibility of a user’s selection of rocking points, vertical reference, scale, and PBR outline during the digitization process. To test for this sensitivity, we performed 100 runs of PBRslenderness on a sample PBR photograph (Fig. 5). Estimates of α_{i} and *R*_{i} were centered closely about their means, with variation coefficients *c*_{v} 4.5%, 6.5%, 4.4%, and 4.3% for α_{1}, α_{2,}*R*_{1}, and *R*_{2}, respectively. This indicates that the single-user estimates of α_{i} and *R*_{i} by PBRslenderness are reproducible to within ∼5% from run to run.

Another source of uncertainty is the dispersion of α_{i} and *R*_{i} estimates between different PBRslenderness users. To assess this uncertainty, a population of 20 PBR photographs was processed by three of us using PBRslenderness (Fig. 6). Dispersions of α_{1} and α_{2} had a maximum of 40% (80% of the estimates were within the 0%–20% dispersion range) and 22.5%, respectively. Dispersions of *R*_{1} and *R*_{2} were as much as 9% and 8%, respectively. The large dispersions in α_{i} relative to *R*_{i} estimates may have been caused by how long each user spent digitizing PBR outlines, different choice of rocking points, plumb bob locations, or scale locations (Purvance, 2005).

### Parameter Reproducibility: 2D versus 3D Parameter Comparisons

To compare the 2D estimates of α_{i} and *R*_{i} made by PBRslenderness with calculations derived from 3D surface models of PBRs, we used a 3D geometric calculator that employs photogrammetrically generated triangular mesh models of real PBRs (Anooshehpoor et al., 2004). The 3D calculator assumes that each triangular facet forms a closed convex surface. Each PBR surface model is thus divided into convex elements whose mass-weighted centroids are computed and combined into a single 3D center of mass. The PBR’s absolute smallest α value (α_{min}) and *R*_{i} lengths are then computed given user-defined rocking points. Using six 3D PBR models and a total of 22 viewing azimuths, we compared the computed α_{i} and *R*_{i} values from the 3D parameter calculator with those estimated in 2D by PBRslenderness (Fig. 7). Our 2D estimates compared well with their 3D counterparts, attaining correlation coefficients R^{2} of 0.86, 0.85, 0.97, and 0.98, and root mean square errors of 4.83°, 5.16°, 0.04 m, and 0.03 m for α_{1}, α_{2}, *R*_{1}, and *R*_{2}, respectively. As observed in the single-user and multiuser sensitivity analyses, the relatively large scatters between the 2D and 3D α_{i} calculations are likely due to variabilities in our selections of rocking points (*RP*_{i}; Purvance, 2005). Similarly, our estimated 2D centers of mass compared well with the calculated 3D centers of mass for each of the 22 viewing azimuths (Fig. 8). However, PBRslenderness appeared to systematically overestimate the heights of the centers of mass, especially for tall and slender PBRs, by a maximum of 8.8% of the distance between the centers of mass and the PBR bases. Despite this overestimation and given the complex geometrical configuration of the 3D PBR models, our reproducibility results demonstrate that PBRslenderness estimates of center of mass, α_{i}, and *R*_{i} compared sufficiently well with the 3D calculations. This demonstrates that PBRslenderness can be effectively used as an efficient α_{i} and *R*_{i} estimator.

### Geomorphic Settings of PBRs

The mean α_{min} value of the PBRs surveyed in the Granite Dells precarious rock zone is 29° with a 38% variation coefficient (Fig. 9). The relatively wide range of slenderness values suggests that the PBRs surveyed in the Granite Dells precarious rock zone have a large variety of stabilities that may be controlled by several geologic and geomorphic factors. These include the local hillslope gradients and upslope contributing areas in which the PBRs reside, local joint densities in the bedrock within which the PBRs formed, and PBR shape. To explore the degree of influence these controls have on PBR stability, we used PBRslenderness to evaluate α_{min} values for the surveyed PBRs as a function of local hillslope gradient, upslope contributing area, and the shapes of the PBRs.

Figure 10A plots slenderness as a function of local hillslope gradient computed over a 3 m lidar-generated DEM. There appears to be no correlation between local hillslope gradient and PBR slenderness, suggesting that PBR stability is not controlled by the steepness of the local hillslope on which it is balanced. This likely reflects that the stability of a PBR is controlled by its local geometry and pedestal arrangement, thus implying that the local geomorphic setting of a PBR does not affect its stability. However, few PBRs appear to be preserved on hillslopes steeper than ∼40° (Fig. 10A), which suggests that steep hillslopes do not have geomorphic conditions that are conducive to forming and preserving PBRs due to unfavorable slope-dependent soil and corestone production rates. As a result, PBRs that formed on hillslopes steeper than ∼40° early in their life cycles may have toppled at a time closer to their formation compared to PBRs that are still preserved. Figure 10B plots slenderness as a function of upslope contributing area, which is a geomorphic proxy for a PBR situation in a drainage basin. The stability of a PBR does not appear to be controlled by its geomorphic location in a drainage basin; PBRs near drainage divides (small contributing areas) are equally as stable as those near the outlets of drainage basins (large contributing areas). However, PBRs appear to form in upslope contributing areas between ∼5 and ∼30 m^{2} per unit contour length, further indicating that geomorphic conditions in the upper reaches of drainage basins are conducive to forming and preserving PBRs. Figure 10C is a plot of slenderness as a function of PBR height to diameter aspect ratio. The closer the aspect ratio is to 0, the more precarious the PBR becomes, and vice versa. For the surveyed PBRs, α_{min} increases with increasing aspect ratio, confirming that slender PBRs are more precarious than their cubic counterparts and indicating that the 3D geometry of a PBR is an important parameter that controls its static stability.

*h*is the slope-normal soil thickness (secant of the slope angle), is the rate of bedrock to soil conversion, ∇ is the divergence vector , is the downslope soil mass flux vector, and ρ

_{s}and ρ

_{r}are the bulk soil and rock densities, respectively. Early field observations by Davis (1892) and Gilbert (1877) suggested that the origin of convex landforms by sediment transport via mass flux is linearly related to slope:

*K*is the diffusion coefficient and

*z*is the ground surface (Fig. 11). At steady-state conditions where soil thickness is constant throughout the hillslope (∂

*h*/∂

*t*= 0), Equation 6 may be simplified to produce the following relationship between soil production and sediment divergence (Heimsath et al., 2001b):

## CONCLUSIONS

We presented a simple and efficient method for determining the 2D static stabilities of PBRs from unconstrained digital photographs. Our method was implemented as a GUI-enabled MATLAB program, PBRslenderness, which is available in the Supplemental File (see footnote 1). However, the methods presented here can easily be implemented in other programming languages that include basic image processing functionalities.

We tested our program using single-user and multiuser analyses, both of which showed that the reproducibility of 2D static stabilities estimated by PBRslenderness are within 5% for single users and 20% for multiple users. We also tested the program’s 2D calculations of PBR stabilities against a 3D PBR stability calculator that uses 3D surface models of PBRs. Comparison analyses between our 2D method and 3D PBR models showed that PBRslenderness estimates the static stabilities of PBRs to within ∼5° and 0.03 m error for α_{i} and *R*_{i}, respectively, while overestimating the height of the center of mass by only ∼9% of the distance from the PBR base. PBRslenderness showed that there is no correlation between the static stability of a PBR and the hillslope angle and drainage area from which it was exhumed. Conversely, PBRslenderness showed that the 3D shape of a PBR strongly controls its static stability. These results demonstrate that PBRslenderness is a useful tool for understanding the geomorphic processes that control PBR formation and preservation in a drainage basin.

There is an ongoing effort to archive PBR images and testing results of ∼15,000 digital images of PBRs that have been acquired over 18 years of PBR research (Biasi and Brune, 2010). In this effort, PBRslenderness has been used to analyze α_{i} and *R*_{i} parameters of >1700 PBRs that have been screened as valuable candidates for constraining earthquake-generated ground motions in southern California. Furthermore, independent validation of PBRslenderness was used by Anderson et al. (2011) with field tests of 41 PBRs in southern California; their 2D static stability analyses made using PBRslenderness and digital photographs of PBRs produced α_{i} and *R*_{i} estimates that were comparable to those computed from PBRs that have been tested in the field. Both uses of PBRslenderness demonstrate its utility in the rapid and accurate estimation of the static stabilities of large PBR populations. Given the relatively simple GUI and minimal computing requirements of PBRslenderness, future PBR investigators are encouraged to use a tablet or hand-held computer and a consumer-grade digital camera to estimate the quasi-static stabilities of large PBR populations directly in the field and prior to performing full fragility analyses on select PBRs.

We appreciate insightful discussions with Rasool Anooshehpoor, Glenn Biasi, Jim Brune, and Kelin Whipple about our two-dimensional method and its validation. Special thanks to Thomas Hanks for his insights into the geomorphology problem of precariously balanced rocks. Virginia Seaver generously provided access to the precariously balanced rocks. This project was supported by the Southern California Earthquake Center (SCEC) Extreme Ground Motion Special Group. The SCEC is funded by a National Science Foundation (NSF) Cooperative Agreement EAR-0529922 and the U.S. Geological Survey Cooperative Agreement 07HQAG0008. This is SCEC contribution 1632. Haddad was supported by a grant that was awarded jointly by the Arizona State University (ASU) Graduate and Professional Student Association, the ASU Office of the Vice President of Research and Economic Affairs, and the ASU Graduate College. The Granite Dells airborne lidar data were collected by the NSF-funded National Center for Airborne Laser Mapping (NCALM) for a Seed Grant awarded to Haddad. We appreciate comments provided by Carol Frost, Antonio Rodriquez, Dylan Rood, Mark Stirling, and an anonymous person, whose reviews served to improve this paper.

^{1}Supplemental File. Zipped file that contains: (1) “PBRslenderness.p” file, which is the program that is presented in this paper and requires MATLAB to run. (2) “PBRslenderness_user_manual.docx” file, which is the user manual for the program. (3) “PBR_sample_photo.jpg” file, which is a sample photograph of a PBR that the reader can use to experiment with. (4) “PBRslenderness_video_tutorial.wmv” file, which is a video demonstration of the program. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00788.S1 or the full-text article on www.gsapubs.org to view the Supplemental File.