Fluvial channels in metamorphic core complexes are preferentially oriented parallel and perpendicular to the direction of tectonic extension. This pattern has been variably attributed to such causes as tectonic tilting during extension, channel elongation by slip along the range-bounding detachment fault, and the exploitation of extension-related joint sets during channel incision. In this paper we use field measurements, digital elevation model analyses, and numerical modeling to test hypotheses for the tectonic and structural control of fluvial channels in metamorphic core complexes, using the Catalina-Rincon core complex in southern Arizona, USA, as a type example. Field measurements and aerial photographic analyses indicate that channels of all sizes exploit steeply dipping joint sets during fluvial incision. As a consequence, channels become preferentially aligned along those joint sets. First and second Strahler-order channels preferentially exploit a joint set oriented perpendicular to the extension direction, while higher-order channels preferentially exploit a joint set oriented parallel to the extension direction. While these observations support the joint-exploitation hypothesis for structural control of drainage architecture, numerical modeling indicates that the spatial distribution of rock uplift during the initial phase of extension plays a crucial role by determining which joint set is preferentially exploited by channels of which Strahler orders. Numerical models indicate that higher-order channels exploit the joint set that is most closely aligned with the direction of initial tectonic tilting, even if that tilting is active for only a short period of time following the initiation of uplift. We conclude that the drainage architecture in the Catalina-Rincon core complex is the result of a combination of joint exploitation and tectonic tilting mechanisms. Structure also plays an important role in controlling the longitudinal profiles of channels in metamorphic core complexes. Channels in the Catalina-Rincon core complex are characterized by structurally controlled knickpoints with a wide distribution of heights and spacings. Field observations indicate that the occurrence of structurally controlled knickpoints and the resulting variability in longitudinal profile form is related to spatial variations in joint density. Numerical models that incorporate spatial variations in joint density using a stochastic bedrock erodibility coefficient are capable of reproducing the statistical properties of longitudinal profiles in the Catalina-Rincon core complex, including the power spectrum of longitudinal profiles and the frequency size distribution of structurally controlled knickpoints. The results of this study illustrate the important roles played by both jointing and the spatial distribution of rock uplift on the geomorphic evolution of metamorphic core complexes. More broadly, the study provides a recipe for how to incorporate joint-related structural controls into landscape evolution models.


A better understanding of the topographic evolution of mountain belts requires a better quantitative understanding of how bedrock channels behave. Early modeling work on the evolution of bedrock channels emphasized the relationships between drainage area, uplift rate, and longitudinal profile form in relatively simple cases of uniform tectonic uplift rate and rock type (e.g., Whipple and Tucker, 1999). More recently, the impact of variations in tectonic uplift rate, sediment supply, and substrate heterogeneity on the morphology of bedrock channels has been emphasized (Finnegan et al., 2005; Wobus et al., 2006; Amos and Burbank, 2007; Gasparini et al., 2007; Pazzaglia et al., 2007; Turowski et al., 2007; Whittaker et al., 2007). Studies have shown that bedrock channels can respond to differential uplift or variations in substrate erodibility in complex ways. Channels adjust both in width and slope to variations in uplift rate or resistance to erosion. As a result, positive or negative feedbacks can occur between erosion rate, channel width, and channel slope (Amos and Burbank, 2007; Whittaker et al., 2007; Attal et al., 2008). Channels can also pass a threshold channel gradient beyond which a reduction in bed shear stress occurs in spite of the bed steepness (Crosby et al., 2007). In such cases, tectonically and structurally controlled knickpoints can persist in the landscape long after they are formed. In this paper we focus on the role of substrate variability on the orientation and longitudinal profiles of bedrock and mixed bedrock-alluvial channels in a tectonically inactive metamorphic core complex in southern Arizona, USA.

Metamorphic core complexes form in areas of continental extension. Geologically, metamorphic core complexes are characterized by the exposure of rocks in the footwall of a low-angle detachment fault (Davis and Coney; 1979; Armstrong, 1982). Topographically, the complexes form spatially discontinuous mountain ranges with relief of as much as 3 km and widths of ~20 km. Fluvial drainages in the metamorphic core complexes of the southwestern U.S. are often aligned preferentially parallel and perpendicular to the extension direction (Spencer, 2000). This apparent tectonic and/or structural control on drainage architecture has puzzled geologists for at least two decades (Pain, 1985; Spencer, 2000). Nevertheless, despite their prominence in the Cordilleran landscape (and in landscapes elsewhere in the world) and Cordilleran tectonics, metamorphic core complexes have been the focus of few published geomorphologic studies. Fewer of these have been quantitative or extended our understanding beyond the conceptual stage. In this paper we address three problems: (1) a classic problem in metamorphic core complex geomorphology (the origin of linear valleys aligned parallel to the extension direction), (2) problems of general relevance to geomorphology in a range of settings (joint control on network geometry and separately on profile shape), and (3) a methodological problem (simulating anisotropic erodibility on a regular grid).

Pain (1985) was the first to describe the geomorphology of metamorphic core complexes and propose a model for how drainages are preferentially oriented in these systems. Focusing on metamorphic core complexes in southern Arizona and southeastern California, Pain (1985) noted locations where initially extension-parallel drainages had been captured by radially oriented drainages. Pain used these observations to propose a two-stage model for the tectonic control of drainage architecture in metamorphic core complexes. In Pain's model, extension-parallel drainages first form on a gently sloping ramp of footwall rocks exposed by tectonic denudation. Extension in metamorphic core complexes is commonly accompanied by orthogonal contraction and antiformal arching (Yin, 1991; Fletcher et al., 1995). Antiformal arching, Pain argued, diverts some of the linear drainages into a radial pattern while preserving the extension-parallel orientation of drainages uplifted on or near the major antiformal axes.

Spencer (2000) pointed out that extension-parallel drainages in metamorphic core complexes are not restricted to the major antiformal axes: as such, there must be some additional mechanism for controlling drainage architecture in metamorphic core complexes. Spencer (2000) argued that this additional mechanism was the repeated incision of channel segments into the footwall ramp after slip along the range-bounding detachment fault. Spencer (2000, p. 727) stated that “maintenance of hydraulic linkage during sequential fault-slip events will guide the lengthening stream down the fault ramp as the ramp is uncovered, and stream incision will form a progressively lengthening, extension-parallel, linear drainage segment.” In addition, Spencer argued that small extension-parallel grooves in the newly exposed fault surface could further guide channels along the extension direction. Spencer's (2000) model is problematic for three reasons. First, in the absence of structural control, channels generally align along the direction of water flow, i.e., the direction of steepest descent. The direction of steepest descent prior to fluvial dissection is determined by the spatial distribution of rock uplift (e.g., tilting) near the range front. That direction may or may not be parallel to the extension direction. Second, Spencer's (2000) model assumes that the channel is incising both upstream and downstream from the range-bounding detachment fault. In general, however, range-bounding faults coincide with the transition from erosion (upstream of the fault) to deposition (downstream of the fault). Variations in the ratio of water and sediment through time cause the transition from erosion to deposition to fluctuate, but it is not generally the case that the range-bounding fault is upstream from that transition, as Spencer's model assumes. Third, the models of Pain and Spencer do not address the occurrence of extension-perpendicular drainages in metamorphic core complexes, first documented by Miksa (1993) and later by Spencer (2000).

Miksa (1993) documented field evidence that channels exploit extension-related joint sets in the Tanque Verde Ridge area of the Catalina-Rincon core complex. Joint exploitation is an appealing hypothesis for preferred drainage orientations in metamorphic core complexes because it is readily observed in the field and in aerial photographs, and because the structural control is continuously exerted as channels incise more deeply into bedrock, rather than being a relict of an earlier control localized in time (e.g., during an early period of tilting) or in space (e.g., at a channel crossing the detachment fault). Miksa (1993) documented a higher frequency of channels oriented both parallel and perpendicular to the extension direction compared to channels oriented in other directions in her Tanque Verde Ridge study area. Lower-order channels are preferentially oriented perpendicular to the extension direction; higher-order channels are preferentially oriented parallel to the extension direction. Joint alignment parallel and perpendicular to the extension direction is a general feature of all of the well-studied metamorphic core complexes of southern Arizona and southwestern California. Extension-parallel and extension-perpendicular joints developed in these ranges as tectonic denudation relieved the remnant stresses imposed on the rocks in Oligocene–early Miocene time. A close association between the orientation of joints and other extension-related structural elements (e.g., fracture zones, dikes) and either the extension-parallel or extension-perpendicular directions was documented by Peterson (1968) in the Catalina forerange, by Reynolds (1982) at South Mountain, by Richard (1983) in the Harquahala Mountains of west-central Arizona, and by Fletcher et al. (1995) in the central Mojave metamorphic core complex. The structural details of each of these ranges are, of course, complex and unique, but in all cases there are structural elements that provide extension-parallel and extension-perpendicular zones of weakness or heterogeneity that channels can exploit as they incise. Due to the general correlation between structural elements and extension-parallel and/or extension-perpendicular directions in metamorphic core complexes of the southwestern U.S., Miksa's (1993) hypothesis provides a mechanism for the structural control of drainage orientation that could be broadly applicable to many metamorphic core complexes, not just the Tanque Verde Ridge area that was the focus of her study. Miksa's hypothesis is also consistent with many papers that have documented joint control of fluvial channel morphology in other tectonic settings (e.g., Hobbs, 1905; Bannister, 1980; Miller, 1991; Eyles et al., 1997; Ehlen and Wohl, 2002; Wohl, 2008). Several questions remain, however, regarding Miksa's (1993) hypothesis. First, it is unclear whether joint control is the sole mechanism for controlling drainage orientation in metamorphic core complexes, or whether joint exploitation acts in concert with some other mechanism. Second, it is unclear why lower Strahler-order channels preferentially exploit the extension-perpendicular joint set while higher Strahler-order channels preferentially exploit the extension-parallel joint set.

Figure 1 illustrates four areas with preferentially extension-parallel and extension-perpendicular drainages from southern and central Arizona. The extension direction in southern Arizona and southeastern California is regionally correlative (Wust, 1986) and is oriented west-southwest or approximately along a 240° azimuth. Of the four areas illustrated in Figure 1, three are considered classic Cordilleran metamorphic core complexes: the Catalina-Rincon core complex, the Tortolita Mountains (located adjacent to the Catalina-Rincon core complex to the west), and South Mountain. Higher-order drainages in all of these metamorphic core complexes are predominantly oriented parallel to the extension direction; lower-order drainages are oriented perpendicular to the extension direction. The Galiuro Mountains, located adjacent to the Catalina-Rincon core complex, was formed during the same late Oligocene–early Miocene extension as the Catalina-Rincon core complex, but underwent a smaller magnitude of extension. The Galiuro Mountains exhibit the same extension-parallel and extension-perpendicular drainage architecture as the other examples, but higher-order drainages are preferentially oriented perpendicular to the extension direction along a south-southeast–north-northwest orientation. The extension-perpendicular drainages in the Galiuro Mountains can be explained by north-northwest–striking normal faults, but the extension-parallel drainages in this range are not well understood.

In this study we use field measurements, digital elevation model (DEM) and aerial photographic analyses, and numerical modeling to test hypotheses for the tectonic and structural control of drainages in the Catalina-Rincon core complex near Tucson, Arizona. In the study area section, we describe the study sites and document the evidence for the exploitation of extension-parallel and extension-perpendicular joint sets by incising channels in the East forerange area of the Santa Catalina Mountains. The conclusions of that section compliment those of Miksa (1993) in nearby Tanque Verde Ridge. Moreover, where joint sets locally deviate from the extension-parallel and extension-perpendicular directions in the East Catalina forerange, the preferential channel orientations deviate similarly, providing clear evidence that channels are joint controlled and orient parallel and perpendicular to the extension direction only when those directions coincide with joint directions. DEM analyses provide further documentation of the correlation between channel orientation and joint direction, with lower Strahler-order channels primarily exploiting the extension-perpendicular joint set and higher Strahler-order joint sets primarily exploiting the extension-parallel set, as observed by Miksa (1993). The numerical modeling discussion presents the results of numerical model experiments designed to determine whether the joint-exploitation mechanism is sufficient to explain the observed drainage architecture in the Catalina-Rincon core complex, or whether some additional mechanism is required. Numerical experiments highlight the important role played by the spatial distribution of rock uplift in determining which joint set a particular channel or set of channels exploits. Higher-order channels in these experiments become preferentially oriented along the joint set that is most closely aligned with the direction of steepest descent, which, in turn, is controlled by the spatial distribution of rock uplift (i.e., south-draining channels exploit the joint set most nearly aligned along the north-south direction). This pattern is consistent with that observed in the Galiuro Mountains but is inconsistent with that of the Catalina-Rincon core complex and other metamorphic core complexes of Figure 1, where higher-order south-draining channels are primarily oriented along the west-southwest–east-northeast extension direction. An alternative numerical model, which includes both joint exploitation and an early phase of extension-parallel tilting, is found to reproduce the observed drainage orientation patterns in metamorphic core complexes. In this model, the early phase of tilting causes higher-order drainages to exploit the joint set parallel to the extension direction. Lower-order drainages exploit the orthogonal, extension-perpendicular joint set, thus leading to drainage patterns similar to those observed in the Catalina-Rincon core complex. Numerical modeling also suggests that Spencer's mechanism of channel orientation by slip along the detachment fault is not a viable mechanism for creating extension-parallel drainages. We conclude that a combination of Pain's (1985) tectonic tilting model and Miksa's (1993) joint-controlled model best explains the observed drainage architecture of the Catalina-Rincon core complex.

The numerical modeling discussion also explores an important, related feature of fluvial channels in metamorphic core complexes, i.e., their complex, stepwise longitudinal profiles. Longitudinal profiles exhibit structurally controlled knickpoints and alternating zones of convexity and concavity over a wide range of spatial scales. In the numerical modeling discussion we show how large-scale (i.e., 1–10 km) structurally controlled channel segments can result from small-scale (i.e., 1–10 m) variations in joint density and hence bedrock erodibility. The model we propose is capable of reproducing the statistical properties of actual longitudinal profiles from the Catalina forerange with respect to both their power spectra and the frequency size distribution of structurally controlled knickpoints. In addition to exploring the controls of joint orientation and anisotropic density on the geomorphology of metamorphic core complexes, this paper also provides a recipe for quantifying the role of anisotropic rock strength in landform evolution models more generally. Given the frequent occurrence of jointed bedrock in nature, the methods of this paper should provide a useful starting point for modeling structural controls on bedrock channel morphology in a broad range of other tectonic settings.



The Santa Catalina Mountains are bounded on the south by the low-angle Catalina detachment fault and on the west by the high-angle Pirate fault. Offset along these faults occurred in two separate intervals of deformation and uplift. In the first, late Oligocene–early Miocene, phase of extension and tectonic denudation, offset occurred primarily along the Catalina detachment fault. Extension along an ~240° azimuth was accompanied by tectonic tilting of an extension-parallel topographic ramp and by antiformal arching along a direction approximately orthogonal to extension (Dickinson, 1991). The antiformal or domal arching in metamorphic core complexes is generally attributed to the isostatic response to tectonic denudation (Spencer, 1984) and its enhancement by the flow of weak lower crust and/or magma emplacement beneath the metamorphic core complex (Block and Royden, 1990). In the second, post–mid-Miocene tectonic event, faulting occurred along the high-angle Pirate fault in a manner similar to the block faulting associated with classic Neogene Basin and Range extension (Davis et al., 2004; Wagner and Johnson, 2006). Extension in this latter tectonic phase was oriented more nearly east-west compared to the earlier phase of extension.

Figure 2 illustrates the topography and geology of the Santa Catalina Mountains. The focus of our study is the Catalina forerange, defined approximately as the mylonite-dominated southern flank of the Santa Catalina Mountains. Figure 2A identifies several examples of large, predominantly extension-parallel drainages, including Pima, Sabino, Bear, Soldier, and Molino Canyons, as well as examples of two large, predominantly extension-perpendicular drainages, Ventana and Esperero Canyons. The central portion of the Catalina forerange is characterized by a topographic low, referred to here as the Sabino saddle. The Sabino saddle coincides with a normal fault, and enhanced weathering of this fault zone is most likely responsible for the formation of the saddle. Smaller folds also exist in the East Catalina forerange (our primary field study area, illustrated in shaded relief in Fig. 2C), including the Sabino Canyon and East forerange folds. These folds do not have a clear topographic expression.

Figure 2B illustrates the geology of the Santa Catalina Mountains (after Dickinson et al., 2001). The rocks of the Catalina forerange comprise the mylonitic Wilderness and Oracle granites (Twm and Yom) at low elevations and the nonmylonitic Wilderness granite (Tm) at higher elevations. The mylonitic rocks are composed of leucogranite sills in which the more felsic Wilderness granite has intruded into the more mafic Oracle granite to form alternating light and dark bands. The higher biotite concentration and smaller mineral grains of the dark bands most likely are responsible for their higher weathering rates, easily observed in outcrop (Laughlin, 1959; Sherwonit, 1974). Individual bands range from several centimeters to several meters in thickness (e.g., Fig. 3A). Joints in the Catalina fore-range are predominantly oriented in the extension-parallel and extension-perpendicular directions (Fig. 3B). Peterson (1968) surveyed the joints in the East Catalina forerange in detail and rose diagrams of the several thousand joints he measured are illustrated in Figure 2C.

Figure 2D plots the longitudinal profiles of two representative channels of the Catalina forerange. The profiles were extracted from a 1 m/pixel light detection and ranging (LIDAR) DEM obtained from the Pima County Flood Control District and hence are capable of resolving steep and short knickpoints along the channel. In the absence of structural controls or episodic tectonism, bedrock longitudinal profiles are typically straight or concave. As these profiles clearly illustrate, channels in the Catalina forerange are characterized by highly irregular, stepwise profiles with alternating convex and concave reaches. At relatively large scales (i.e., ≥1 km), the channels of the Catalina forerange are characterized by steep (>10% slope), bedrock-dominated reaches separated by more gently sloping (5%–7% slope), alluvially dominated reaches. Within each of the bedrock-dominated reaches there are multiple, steeply dipping (45° to vertical) bedrock steps with heights ranging from 1 to 10 m (e.g., inset plot in Fig. 2D). Pools commonly form at the base of these steps. These pools provide important habitats for the fauna of the Santa Catalina Mountains, including the frogs that give the mountains their native Tohono O'odham name, Babad Doag (Frog Mountain). The alluvial reaches separating the zones of bedrock knickpoints vary in spacing, but these reaches are on the order of several hundred meters to 1 km. In this paper we refer to the zones of bedrock knickpoints as megaknickpoint zones to distinguish them from the smaller, steeper knickpoints of which they are composed. Given that there has been no significant offset along the range-bounding Catalina detachment fault for ~20 Ma, it is likely that the stepwise nature of these channel profiles is due to structural control rather than episodic tectonism or base-level change.

Field Observations

We conducted detailed structural and geomorphic mapping along portions of Soldier and Molino Canyons in the East Catalina forerange in order to test the hypothesis of structural control of fluvial channel process and form at the channel reach scale. We observed that, in bedrock knick-points, one or both channel walls commonly coincide with a joint surface. Figure 4 illustrates four such examples from Soldier and Molino Canyons. Field surveys provided dozens of such examples along the ~5 km of channel reaches we surveyed. Abrupt downcutting without joint exploitation could also form a planar surface on one or both channel walls, but the fact that planar channel walls in our study area occur at dips and orientations similar to those of nearby joints indicates that planar channel walls are, in nearly all cases, joint surfaces. We also observed that bedrock knickpoints are generally associated with the occurrence of unusually thick leucogranite sills (i.e., light-toned bands). These sills tend to have a lower joint density than the Oracle granite (i.e., dark-toned bands) into which they are intruded. The relationship between band thickness and joint density we observed is qualitatively similar to the relationship between bed thickness and joint density well documented in sedimentary rocks (e.g., Gross et al., 1995). Upstream of these zones of low joint density, the channel is commonly backfilled with alluvium and has a gentle slope (~5‰–7%). Inward from these zones, the channel incises abruptly and steepens to slopes that are commonly >45° as it is forced to incise through locally more resistant rock units. Within more gently sloping (<45°) bedrock reaches, the channel was often observed to be currently occupying and incising several parallel or subparallel joint surfaces, partitioning flow between the joint surfaces. As incision proceeds, presumably one joint surface will incise faster than the others and will capture more of the upstream flow in a positive feedback until one joint surface becomes dominant. Within the steepest knickpoint reaches, there is evidence that the channel switches from plucking-dominated erosion, which is common on more gently dipping bedrock reaches, to saltation-abrasion–dominated erosion. Within saltation-abrasion–dominated reaches, channel walls are not planar, but instead have undulating walls typical of slot canyons (Wohl et al., 1999). We infer that, when the spacing between joints increases locally above a threshold value (locally equal to several meters), channels are unable to pluck bedrock from the bed and must erode by saltation abrasion and/or cavitation. Extreme-flood flow depths along Soldier and Molino Creeks were estimated to be 2–3 m based upon the height of the highest evidence of channel wall scour (at higher locations, rock surfaces are more highly weathered and evidence of varnish and staining are present). Given these flow depths, it is reasonable to conclude that the channel is unable to pluck blocks of several meters in length and width from its bed. In such locations, saltation abrasion is the sole mechanism of erosion and the channel must steepen to offset the lower effectiveness of this process compared to plucking.

Along each of the study reaches, measurements were made of channel slope, channel width, rock hardness, and joint orientation, dip, frequency, and spacing. Figure 5 illustrates the measurements made along the Soldier Creek study area, where the upstream drainage area is ~6 km2. Joint data are reported in Figure 5 using rose diagrams. Joints that were spaced by less than a few centimeters were not counted separately from neighboring joints.

The channel reaches along the Soldier Creek study area are strongly joint controlled, with every channel segment corresponding closely to either the west-southwest–east-northeast– or south-southeast–north-northwest–trending joint sets. Transitions from one channel segment to another are almost all precisely at right angles. The channel geometry is characterized by a strong inverse correlation between channel slope and width, with steep bedrock reaches corresponding to channel widths typically <1 m. Local channel slopes are generally higher and channel widths narrower when the channel is oriented along the west-southwest–east-northeast–trending joint set than along the south-southeast–north-northwest–trending joint set. The difference in mean spacing between the two joint sets is one possible reason for these differences in channel slope. Joint density was approximately twice as great along the south-southeast–north-northwest–trending joint set compared to the west-southwest–east-northeast–trending joint set, hence lower slopes are required along this direction in order to incise at similar rates. Variations in rock hardness were minimal along these reaches. Schmidt Hammer readings of 62–64 were measured, the value for each reach obtained by acquiring 10 measurements and then computing an average value for that reach. Despite the compositional variation between light and dark bands, the two band types clearly have similar hardness values for fresh, unweathered rock recently exposed by channel incision. The Schmidt Hammer data provide evidence that variations in bedrock channel erodibility in the Catalina forerange are caused primarily by variations in joint density, not by variations in composition.

Figure 6 illustrates the measurements made along the Molino Creek study area, where the upstream drainage area is 13 km2 (i.e., twice as large as the Soldier Canyon study area). We measured structural and channel morphology data along the three megaknickpoint zones of Molino Creek located in Figure 6A. Broad alluvial zones occurred between each of these zones, precluding any observation or measurement of the bedrock structure beneath the channel bed. Channel orientations in the Molino Creek study area correlated with joint sets, although the correlation is not as strong as in the Soldier Creek study area. This is partly the result of the larger size of the Molino Canyon study area, which includes three megaknickpoint zones over an along-channel distance of >4 km. Also, a north-south–trending joint set was present in the southern portion of the study area and locally controlled the channel orientation along reaches in that area. The overall channel orientation is southwest (with a west-southwest orientation in the northern end of the study area and a nearly south orientation in the southern end of the study area). Along Soldier Creek, the steepest bedrock reaches are oriented along the west-southwest–east-northeast–trending joint set. Along Molino Creek, in contrast, two of the three megaknickpoint zones are associated with channel orientations along the south-southeast–north-northwest–trending joint set. This difference between the two study areas was not readily explained by the joint spacing data, which were broadly similar between the two areas. As Miksa (1993) pointed out, joint-controlled erosion may also be influenced by joint length or continuity in addition to spacing. As such, subtle differences in joint continuity could lead to significant differences in erodibility and hence channel slope. Overall, the longitudinal profile of Molino Canyon through the study area illustrated in Figure 6 is convex. This large-scale downslope-steepening trend could be the result of a more southward orientation to the channel in the southern end of the study area relative to the more southwestward orientation of the northern end of the study area (which, therefore, might be more capable of exploiting the primary west-southwest–east-northeast–trending joint set compared to the southern end). The farthest downslope megaknickpoint zone could also be influenced by the east forerange fold, the axis of which nearly coincides with the upstream boundary of that zone. However, any effects of the fold on channel morphology are subtle; there is no clearly identifiable influence of the fold on either channel slope or width through the study area.

DEM and Aerial Photographic Analysis

DEM and aerial photographic analyses were performed in order to further document the relationships between channel frequency and orientation for channels of different sizes in the Catalina forerange. Figure 7A illustrates the channel frequency versus orientation data of Miksa (1993) from the Tanque Verde Ridge study area. Miksa digitized channel segments from topographic contour maps and plotted the frequency of channel segments in the Tanque Verde Ridge study area as a function of channel orientation, and found that channels of low Strahler order (1 and 2) are preferentially oriented along the joint set perpendicular to the extension direction, while channels of high order (3 and 4) are preferentially oriented parallel to the extension direction.

In this study we used an automated technique to perform a similar channel frequency versus orientation analysis in the Catalina forerange. Rather than digitize channel segments by hand as Miksa did, we used an automated procedure based on DEM analysis. In our analysis, the contributing area of each pixel was first calculated using the D8 (steepest descent) routing method. Second, the local channel orientation was determined by identifying the lowest pixel among the next-nearest neighbors of each pixel. Each pixel in a DEM has 8 nearest neighbors (including diagonals) and 16 next-nearest neighbors. Channel segments defined using nearest neighbors are restricted to multiples of 45°. By using next-nearest neighbors, the angular resolution of the analysis is improved but the results are still based on a discrete set. Nevertheless, this approach has two advantages over traditional map-based methods: ease of use and the ability to include tens to hundreds of thousands of channel segments in the analysis. Figure 7B illustrates rose diagrams of channel frequency for three sets of drainage areas: 10−3–10−2 km2, 10−1–1 km2, and >10 km2. These drainage area groups do not correspond precisely to particular Strahler orders, but they were chosen because they overlap closely with Strahler orders 1, orders 2–3, and order 4, respectively, based on comparing maps of drainage area with Strahler order. The rose diagrams shown in Figure 7B show results similar to those of Miksa (1993) in Figure 7A, i.e., channels of low Strahler order (1 and 2) are preferentially oriented along the joint set perpendicular to the extension direction, while channels of higher order (3 and 4) are preferentially oriented parallel to the extension direction. The patterns in the Catalina forerange are somewhat less defined than those Miksa found in the Tanque Verde Ridge area. This more diffuse pattern has three likely causes. First, the correlation between channel and joint orientations is less strong in alluvial reaches, where alluvial backfilling creates a broad floodplain for the channel to occupy, rather than a series of closely spaced joint surfaces. Because Tanque Verde Ridge has fewer zones of alluvial storage, there are more joint surfaces exposed that can act as structural controls on channel orientation. Second, the Santa Catalina Mountains have a larger area that includes more variability in structure compared to Tanque Verde Ridge. Third, the digital method itself could be partly responsible for the less defined pattern in the Catalina forerange.

The automated nature of DEM analyses makes it straightforward to test for other correlations between channel morphology and orientation. In addition to channel frequency, we also calculated the average channel slope as a function of channel orientation and drainage area only. The results (not shown) indicate that channel slopes are slightly less steep, on average, where channels are oriented parallel and perpendicular to the extension direction compared to when they are oriented in other directions (e.g., directly north-south). This result was surprising, at first, since we expected that channels would require significantly lower slopes to incise if they were taking advantage of one of the primary vertically dipping joint sets along the west-southwest–east-northeast and south-southeast–north-northwest directions. There are two likely explanations for this weak correlation between channel slope and orientation. First, the large majority of channel segments do exploit joint sets, either one of the two primary joint sets (i.e., west-southwest–east-northeast or south-southeast–north-northwest) or a secondary joint set (e.g., north-south). In instances where the channel is not oriented along the two primary joint sets, the channel is often aligned along a minor joint set that, while not as common as the two primary joint sets, locally plays a similar role in enabling the channel to incise at relatively low slopes. Second, field observations indicate that channels respond to lower joint densities through changes in both channel slope and width. When local joint densities decrease below a threshold value of several meters, channels often narrow by as much as an order of magnitude or more (e.g., from several meters to <1 m) in order to focus their stream power sufficiently to incise by saltation abrasion. The fact that channel slope is not very sensitive to orientation suggests that variations in channel width are an important mechanism, equal to or greater than variations in slope, by which channels locally achieve erosion rates comparable to those of nearby reaches despite their relatively low joint densities.

Aerial photographic analysis, while not as precise as field mapping, provides corroborating data that support joint-controlled channel incision in the Catalina forerange. Figure 8 illustrates three areas in the Catalina forerange where the predominant joint orientation can be inferred from 1 m/pixel aerial photographs. Vertically oriented joints provide conduits for runoff and hence undergo enhanced weathering and the formation of preferential planes for slope failure. Over time, these zones of enhanced weathering can be observed in aerial photographs and be used to infer the presence of joints. In each image of Figure 8, primary joints sets were identified by measuring the orientations of bedrock fractures visible in the photographs. White lines in Figure 8 highlight the locations where the local channel orientation closely aligns with one of the two primary joint sets. In Figures 8A and 8B, the orientation of the primary joint sets is within 15° of the regionally averaged joint orientations, which, in turn, closely align with the extension-parallel and extension-perpendicular directions. In Figure 8C, however, the primary joint sets are displaced ~45° counterclockwise from the regionally averaged orientations. In this area, the drainages maintain their alignment with the local joint sets, not the regional extension direction. The results of Figure 8C clearly show that channel orientations follow joint orientations locally within the Catalina forerange. It is only in those areas where joint orientations happen to be parallel and perpendicular to the extension direction that channel orientations are correlated with extension direction.


Three-Dimensional Model with Orientation-Dependent Bedrock Erodibility

The numerical models of this paper are based on the stream-power model of bedrock erosion, which in its basic form can be written as:  
where h is elevation, t is time, U is the rock uplift rate, K is the coefficient of bedrock erodibility, A is the contributing area, x is the along-channel distance, and m and n are empirical coefficients (Whipple and Tucker, 1999). The stream-power model is most applicable to bedrock channels dominated by plucking, because it is in those channels where erosion rates can be related most directly with contributing area, a proxy for peak discharge (Whipple et al., 2000). In addition to the stream-power model, bedrock landform evolution models often include a hillslope erosion component. In this paper we do not include a hillslope erosion component because the pixel sizes of the model (100 m) are larger than the average distance between divides and channel heads in the study area. As such, bedrock channel erosion is the dominant process in the model down to the scale of individual pixels. The contributing area A in the model is calculated at each time step using the multiple-flow-direction method of Freeman (1991). This algorithm is implemented by first ranking all pixels in the basin from highest to lowest in elevation. Starting with the highest grid point in the basin, which receives only local runoff, flow is distributed to all of the neighboring downslope pixels, weighted by slope. Routing is next performed for the second-highest grid point in the basin, then proceeding in rank order down to the lowest grid point. This method ensures that all incoming flow has been accounted for before the flow is distributed downstream.
Joint exploitation is included in the model by prescribing the bedrock erodibility, K, to be a function of channel orientation:  
and K0 is the maximum erodibility value, ε1 and ε2 are parameters that control the magnitude of the anisotropy, θ and ϕ are the local channel orientation and strike of the primary joint set, respectively. Equation 3 is a two-parameter function that mimics the form of the rose diagram of joint frequencies observed in the Catalina-Rincon core complex. Equation 3 also assumes that the principal and secondary joint sets are orthogonal. If they are not, a second joint orientation can be introduced into equation 3 in place of the term ϕ – π/2 in order to represent the secondary joint set. The parameter ε1 controls the strength of the anisotropy in K between channels aligned parallel to the primary and secondary joint sets and those not aligned parallel to these joint sets. Specifically, the ratio of the maximum to the minimum erodibility is given by 1/ε1. The parameter ε2 controls the magnitude of the difference between the erodibility values for channels aligned along the primary versus secondary joint sets. For example, a value of ε2 = 0.5 means that the bedrock erodibility coefficient along the secondary joint set is equal to half its value along the primary joint set. As an example, the rose diagram in Figure 9B plots the orientation dependence of K for ε1 = 2.5, ε2 = 0. As a caveat, equation 3 is only one possible approach to incorporating the effects of joints into bedrock-channel landscape evolution models. A more sophisticated model that explicitly incorporates data on joint spacing, length, and continuity would likely be the most ideal approach to this problem. Nevertheless, in this study we adopted a more empirical approach, effectively linking joint frequencies directly to bedrock erodibility values.

Introducing an orientation dependence into a landscape evolution model raises a problem with regard to the underlying grid structure of the model. In landscape evolution models that operate on square pixels, the local channel slope and orientation are usually computed using the eight nearest neighbors to each pixel (including diagonals). This becomes problematic, however, when introducing an orientation dependence into K because channel orientations computed in this way are restricted to multiples of 45°. This fundamental discreteness resulting from a regular grid can be mitigated in two ways. First, the channel orientation and slope can be computed at a larger spatial scale using neighbors separated by one or more pixels between them (i.e., next-nearest neighbors or next- to next-nearest neighbors). Alternatively, the D ∞ algorithm can be used (Tarboton, 1997). The D ∞ algorithm identifies the two adjacent pixels amongst the nearest neighbors that have the lowest elevations. The elevations of those two pixels are then used to identify a channel orientation that can take on any continuous value between these two lowest adjacent neighbors. In the model of this paper we use the first method by calculating the channel slope and orientation using the next-nearest neighbors. This approach reduces, but does not eliminate the discreteness in channel orientation computed by the model. For the purposes of this paper, however, this method was found to work well.

The results of three model experiments designed to illustrate the effects of orientation-dependent bedrock erodibility are shown in Figure 9. Model results are shown for U = 1 m/ka, K = 10−2 ka−1, m = 0.5, n = 1, L = 40 km, and Δx = 100 m, where L is the width and length of the domain and Δx is the pixel size. At the beginning of each model, the initial topography is a flat surface at h = 0 with a white noise of 1 m root-mean-squared amplitude superimposed. This white noise provides small-scale variations in flow pathways that represent the small-scale variations in topography or erodibility present in any area subject to incipient uplift. Each of the model experiments in this paper was run for 2 Ma. This duration was found to be sufficient for channel incision to propagate from the range front to the headwaters of the model domain near the center of the grid. Slope-area relationships for bedrock rivers indicate that m/n ≈0.5 (Whipple and Tucker, 1999; Snyder et al., 2000). Here we choose the linear case of the stream power model (resulting in m = 0.5 and n = 1), but we verified that similar results were obtained for both linear and nonlinear model experiments. The values of U and K in the model control the elevation attained by the model for a given duration (i.e., higher values of K for a given U and model duration result in lower peak elevations). The values used in these experiments resulted in topography with a maximum elevation of ~2 km above base level, which is comparable to the relief of metamorphic core complexes of the southwestern U.S. The value of U can be increased or decreased, and the resulting model topography is essentially unaffected. Only the time required to propagate knickpoints is affected by the absolute value of U as long as the ratio of U and K remains constant.

In the control case with ε1 = 0 and ε2 = 0, illustrated in Figure 9A, the drainages exhibit no structural control. The topography illustrated in Figure 9A is qualitatively similar to the output of stream-power-based landscape evolution models that have been described in the literature since the early 1990s (e.g., Willgoose et al., 1991; Howard, 1994). In contrast, Figure 9B illustrates the results of the model with ε1 = 2.5, ε2 = 0, and ϕ = 240°. These model parameters correspond to two equally erosive joint sets oriented at ϕ = 240° (west-southwest–east-northeast) and ϕ = 150° (south-southeast–north-northwest). Figure 9B clearly shows that channels tend to become aligned along the directions of peak erodibility values (which represent joint sets in the model). It is also apparent from Figure 9B that the choice of which joint set a particular set of channels exploits also depends on the channel orientation as established by the spatial pattern of uplift and/or the shape of the model domain. For channels draining to the bottom and top of the model domain, for example, the lower-order channels are predominantly aligned along the west-southwest–east-northeast joint set while the high-order channels are predominantly aligned along the south-southeast–north-northwest joint set. In the higher-order channels draining the Catalina forerange (e.g., Pima, Ventana, Esperero, Sabino, Bear, Soldier, and Molino Canyons), the majority are aligned along the west-southwest direction, not the south-southeast direction. We conclude from this that joint control alone is insufficient to produce the drainage architecture observed in the Catalina forerange. Increasing or decreasing the value of ε1 in the model increases or decreases the apparent structural control somewhat, but does not change the results qualitatively.

In addition to the anisotropy in bedrock erodibility caused by preferential joint orientations, there is another, more subtle, geomorphically relevant anisotropy that develops in metamorphic more complexes. As described in the previous section, mylonitic rocks are composed of alternating bands (i.e., sills and the country rock into which they intrude) that strike along the extension direction. If joint densities differ between bands and folding occurs, elongated zones of relatively high and low joint densities aligned parallel to the extension direction are exposed in plan-form. This joint-density effect is separate from, but complimentary to, the joint-orientation effect represented in equation 3. This effect is likely enhanced by the presence of corrugations. Corrugations are large-scale folds (wavelengths of kilometers to tens of kilometers) that develop during exhumation. If the folding occurs late in exhumation, it can promote extension-parallel drainages via the tectonic effect of antiformal arching described in Pain's (1985) conceptual model. Folding also causes an additional structural effect in which mylonitic bands (i.e., laterally continuous zones of similar joint density) become exposed in planform, where they can be more readily exploited by incising channels. In the absence of corrugations, rivers will erode vertically through strong and weak bands, and this effect will likely be less significant in terms of generating preferential drainage orientations.

In order to model this effect, we created a random noise for K with a mean value of 10−2 ka−1 uniformly distributed between 0.5 × 10−2 ka−1 and 2 × 10−2 ka−1 with spatial correlations in the extension-parallel direction and no correlations in the extension-perpendicular direction (Fig. 9C). This random field was then input into the stream-power model, which produced the output shown in Figure 9D. As in the model output illustrated in Figure 9B, higher-order, south-draining channels tend to become aligned along the south-southeast–north-northwest–directed joint set, not the west-southwest–east-northeast joint set, as observed in the Catalina forerange. As such, this effect does not appear to be sufficient to reproduce the drainage patterns observed in the Catalina forerange. The degree of anisotropy produced by this model depends on the standard deviation of K prescribed when generating the random noise; greater variation in K causes more anisotropy. One way to identify this effect separately from the joint-orientation effect represented by equation 3 is to note that bedrock erodibility is a function of channel orientation and strike direction, i.e., K(θ, ϕ), in equation 3 while in the case of Figures 9C and 9D it a function of space, i.e., K(x, y) and the channel orientation plays no explicit role. It is unclear how to distinguish these two joint-controlled using the topography or field observations alone, hence more research is needed. In the subsequent model experiments, we include the effects of jointing using equation 3 only, bearing in mind that joint control can be exerted in two different ways.

Figures 10B and 10D illustrate the results of two numerical experiments designed to test the tectonically driven hypotheses of Spencer (2000) and Pain (1985), respectively. Both model results use the same model parameters and forcings as Figure 9A, except where noted. To test Spencer's hypothesis, we imposed spatially uniform uplift on the range, as in the control experiment of Figure 9A, but we simultaneously extended the portion of the model domain subjected to uplift along a southwestward direction over time, as illustrated schematically in Figure 10A. Clearly, spatially uniform uplift is a simplification that is rarely, if ever, achieved in nature, and it is unlikely to be appropriate for any metamorphic core complex. Nevertheless, it is useful to consider simplified end-member scenarios in order to isolate the effects of specific end-member mechanisms (e.g., mountain-front elongation versus tectonic tilting). Joint-controlled erosion with ε1 = 2.5 and ε2 = 0 was also incorporated in this model, as in Figure 9B. The results for this model scenario are shown in Figure 10B as a grayscale map of the topography at the end of the simulation (t = 2 Ma). The model results clearly illustrate that channels do not align in the extension-parallel direction. Instead, the model develops drainages that are perpendicular to the mountain front (i.e., south and west). This was also verified by performing numerical experiments with other mountain-front geometries; drainage networks develop perpendicular to the mountain front in all cases. In contrast, most of the higher-order drainages of the Catalina forerange are oriented along a west-southwest–east-northeast direction even though the mountain front is oriented in a west-northwest–east-southeast direction. The results of this numerical experiment indicate that elongation of the mountain front does not, as a general rule, create a preferred drainage orientation parallel to the direction of elongation or extension.

In the second tectonically driven model scenario, designed to test the hypothesis of Pain (1985), the model domain was first subjected to a short (i.e., 0.2 Ma) interval of tilting down to the west-southwest (240° azimuth). Following the initial phase of tilt-style uplift, a second, longer (1.8 Ma) interval of antiformal arching was imposed on the model domain with the arch axis following the same west-southwest–east-northeast orientation as the initial tilting. In nature, extension and antiformal arching in metamorphic core complexes are essentially coeval. In these numerical experiments, however, we separated extension-parallel tilting and antiformal arching into two separate phases in order to simplify the experiments, and because, in nature, extension-parallel tilting most likely dominates the pattern of rock uplift early in metamorphic core complex development. Contour maps of uplift rates for this model scenario are illustrated in Figure 10C. An orientation-dependent bedrock erodibility was also included in this model scenario to model the effect of Pain's model in combination with joint-controlled erosion. The results of this model, illustrated as a grayscale map of topography in Figure 10D, shows that drainages align parallel to the direction of initial tilting despite the fact that that tilting is only active for a small fraction (0.1 or 10%) of the model duration. The results of this model suggest that the initial spatial distribution of uplift plays an important role in controlling the subsequent drainage architecture because it is during this initial phase of rock uplift that incised channels become a permanent part of the landscape. Subsequent (or coeval) antiformal arching causes some degree of channel reorganization into a radial pattern, but the extent of that reorganization is not sufficient to counteract the effect of the initial tilt-style uplift. We experimented with varying the duration of the initial phase of tilt-style uplift in the model by varying the duration of the initial tilt phase from just a few percent of the total model duration to as much as 50%. We found that as the duration of the initial tilt phase falls below 10% of the model duration (i.e., 200 m maximum relief generation given maximum uplift rates of U = 1 m/ka), the effects of initial tilting become less significant relative to the effects of subsequent antiformal arching. As such, the initial phase of tilting must be sufficiently long to cause channels to incise to some minimum depth, below which they can be reoriented by a subsequent change in the spatial distribution of uplift. Once incision takes place beyond that minimum depth, however, the model results suggest that it is difficult for subsequent tectonic adjustments to substantially reorganize the drainage pattern.

Figure 11 illustrates the results of applying the different structural and tectonic model scenarios to a model domain equal in size and shape to the actual Santa Catalina Mountains. It is useful to consider the model predictions in a domain similar to that of the study area because the results can then be more directly compared to the modern topography of the Catalinas (shown in Fig. 11A using the same scale as the model results in Figs. 11B and 11D) and because the shape of the boundaries of landscape evolution models can, in some cases, control the model results in complex and unexpected ways. The 1100 m contour was chosen as the domain boundary for these models based on the fact that this elevation closely approximates the bedrock-alluvial contact at the mountain front of the actual Santa Catalina Mountains. However, is only an approximation to the actual location of the mountain front, because the mountain front does not follow contour precisely. Figure 11B illustrates the results of the model with joint control only (assuming ε1 = 2.5, ε2 = 0, and ϕ = 240°). In these figures, the actual basin topography is shown for elevations below 1100 m (i.e., downstream from the base-level boundary condition imposed in the model). Uplift is assumed to be spatially uniform in this model scenario. The results of the model are qualitatively similar to those of the same model performed on a square model domain (Fig. 9B). That is, in the area of the Catalina fore-range, higher-order channels are preferentially oriented along the south-southeast–north-northwest joint set while lower-order channels are preferentially oriented along the orthogonal west-southwest–east-northeast joint set. In Figure 11D, the effects of joint-controlled incision are combined with the tectonic model of Pain (1985), illustrated schematically in Figure 11C. Higher-order channels draining the southern and western flanks of the range are more preferentially aligned along the observed extension direction in this case as a result of the additional effect of the initial phase of tectonic tilting (not present in the model of Fig. 11B). Low-order channels in the model are preferentially oriented along the south-southeast–north-northwest joint set, as illustrated by the rose diagram of channel orientation illustrated in Figure 7C. These rose diagrams closely match the rose diagrams of channel orientations from the Tanque Verde Ridge and Catalina forerange, also illustrated in Figure 7A.

The model results illustrated in Figure 11 greatly simplify the tectonic, structural, and geomorphic processes of the range and do not represent a comprehensive model for the evolution of the Santa Catalina Mountains. For example, the topography of the Santa Catalina Mountains is controlled by high-angle faulting along the Pirate fault in ways that are difficult to determine and that we have not attempted to include in the model. Despite the model simplifications, however, the results illustrated in Figures 9–11 strongly support the conclusion that drainages are oriented parallel and perpendicular to the extension direction in the Catalina-Rincon core complex due to a combination of joint-controlled bedrock channel erosion and an early phase of tectonic tilting that imposed an extension-parallel drainage architecture on the higher-order channels draining the southern flank of this metamorphic core complex.

Two-Dimensional Model with Random Spatial Variations in Erodibility

Field observations (described in the study area discussion) indicate that knickpoints in the East Catalina forerange are controlled by spatial variations in bedrock erodibility related principally to local variations in joint density. Field observations indicate that knickpoint formation tends to be associated with the exposure of unusually thick (e.g., 1–3 m) leucogranite sills. It is these sills or light-toned bands that tend to have the lowest joint density and hence can be expected to have the lowest bedrock erodibility or K values. When incising channels encounter these bands, the channel must steepen locally relative to nearby channel segments in order to maintain comparable erosion rates along the longitudinal profile. Local channel steepening is required to maintain uniform erosion rates because higher values of stream power are required to pluck the larger blocks associated with channel segments of low joint density compared to nearby channel segments of higher joint density. If the local joint density falls below a certain threshold value, the channel may also be forced to erode by saltation abrasion, a process that is less efficient than plucking (in the field, evidence for channel erosion by saltation abrasion is restricted to the steepest channel reaches, hence plucking is likely the dominant process in all but those reaches). In this subsection we explore this conceptual model for the formation of structurally controlled longitudinal profile development using a simple mathematical model that treats bedrock erodibility as a random variable. The purpose of this model is not to predict precisely where knickpoints will occur in a given study area, but rather to understand the statistical properties of variations in longitudinal profile form. The longitudinal profiles of channels in the Catalina forerange all exhibit variations over a wide range of scales (Fig. 12C). Despite this apparent complexity in form, the statistical properties of structurally controlled longitudinal profiles in the Catalina forerange show remarkable similarities that point to a common underlying mechanism for structural control in these channels.

The stream power model with spatially variable erodibility can be written as:  
using m = 0.5 and n = 1. In this model, we represent spatial variations in joint density and the effect of those variations on erodibility using a stochastic model for K(x). In this model, K(x) is a random variable with values sampled at a prescribed interval Δx along the channel from a prescribed probability density function (pdf) with minimum and maximum values Kmin and Kmax, respectively. For example, K(x) = K0 for 0 < x ≤Δx, K(x) = K1for Δx < x ≤ 2Δx, where K0 and K1 are samples from the prescribed pdf. We assume, for simplicity, that the values of K are uniformly distributed between the minimum and maximum values. If we further assume that erosion and isostatic rock uplift are both uniform and approximately in balance, then equation 4 can then be integrated to give  
where E is a constant. The model leading to equation 5 is based solely on the stream-power model for bedrock channel erosion. Channels in the Catalina forerange, however, are of mixed type, including both bedrock and alluvial reaches. Bedrock reaches in the Catalina forerange have a range of slopes from 0.1 (10%) to vertical, while alluvial reaches have lower slopes of between 0.05 and 0.07 (5%–7%). In order to include the effects of local alluvial storage in the model, we took the output of the stream-power model from equation 5 and backfilled all reaches where the slope was below Smin, a prescribed minimum value for alluvial storage. The resulting model is defined by five parameters: the ratio E/A1/2, which has units of one over time, Kmin and Kmax, which also have units of one over time, Δx, which has units of length, and Smin, which is dimensionless.

Figures 12 and 13 illustrate the results of the model for different values of the input parameters, along with comparisons of the model to actual profiles in the Catalina forerange. Figure 12A plots the values of K obtained from sampling a uniform distribution between a minimum value Kmin = 0.01 ka−1 and a maximum value Kmax = 1 ka−1 using a horizontal step size of Δx = 10 m. The form of equation 5 indicates that the local channel slope is controlled by the inverse of K, not by K, so in Figure 12B we plotted K−1 for the same values of K shown in Figure 12A. When the values of K−1 are plotted, it is clear that, even though the values of K have a simple (uniform) distribution, the values of K−1 have a skewed distribution capable of producing values of K−1 that are much larger than its mean value. It is these extreme values of K−1 that cause the channel to hang locally and form structurally controlled knickpoints. Figure 12D plots the results of the model assuming E/A1/2 = 0.02 ka−1. In the limit where Kmin equals Kmax, the model predicts a channel segment of constant slope given by S = E/(A1/2K). As the difference between Kmin and Kmax increases in the model, the variability in longitudinal profile form and local slope also increases, and the profiles become more stepwise. The model reproduces all of the qualitative features observed in the channels of the Catalina forerange, including steep, bedrock-dominated reaches characterized by a wide range of step sizes, alternating with gently sloping reaches with alluvial storage. For comparison, Figure 12C illustrates the longitudinal profiles of all channels >10 km in length draining the Catalina forerange. Alternating reaches of convexity and concavity and a stepwise character over a range of spatial scales is a universal characteristic of these channels.

Equation 5 states that the channel longitudinal profile is the integral of a random function. The simplest example of such a function is the random walk. In a simple random walk, the elevation can increase or decrease with equal probability by a prescribed amount Δh for each horizontal increment Δx. Because the probability of an increase or decrease in elevation is equal for each horizontal increment, the average elevation of the random walk is always zero. Due to the probabilistic nature of the model, however, the random walk will exhibit excursions from its mean value. In the stream-power model with spatially random erodibility, megaknickpoint zones (and their complement, the intervening reaches with alluvial storage) are the equivalent of the excursions of a random walk. Conceptually, mega-knickpoint zones form in the model when the cumulative bedrock erodibility falls below the mean value, which causes the channel to locally steepen above the mean channel slope. Conversely, when the cumulative bedrock erodibility is higher than average, the channel decreases in slope and becomes locally concave. If the channel slope falls below the threshold Smin, alluvial backfilling occurs. The exact locations of these excursions cannot be predicted, but their statistical properties can be predicted. A key goal of the modeling in this subsection is to predict the heights and spacings of the megaknickpoint zones as a function of the model parameters, in order to better understand the controls on the geometry of structurally controlled longitudinal profiles in our study area and in similar areas with locally variable joint density.

If the probabilities of increase and decrease in a random walk are unequal, a linear trend will be superimposed on the random walk, causing the walk to trend up or down depending on which direction is favored probabilistically. In an analogous way, the elevation always decreases with increasing distance downstream in our model. To compute the average slope of the longitudinal profile predicted by the model, it is necessary to integrate the local slope over all possible values of K, weighted by the probability density function. This gives:  
Using the parameters chosen for the example given in Figure 12D, (6) gives Sav = 0.093034. This value is comparable to that of channel profiles in the Catalina forerange, which rise ~100 m for every 1 km of horizontal distance.

Several examples of the model output are presented in Figure 13B for the chosen parameters Δx = 10 m, Kmin = 0.01 ka−1, Kmax = 1 ka−1, E/A1/2 = 0.02 ka−1, and Smin = 0.07. Actual longitudinal profiles from the Catalina forerange are plotted in Figure 13A for comparison. The profiles plotted in Figure 13B are each offset by 100 m so that they can be more easily distinguished in the figure. The stepwise character of the observed and modeled profiles can be highlighted by plotting a moving average (with a 200-m-wide window) of the slope of each profile (Figs. 13C, 13D). Note that the scale on the y axis of Figures 13C and 13D is divided so that all of the plots can be shown on the same curve. Both the actual and the model profiles have megaknickpoint zones (zones where the average slope exceeds 10% over length scales of 200 m) with typical spacings of ~0.5–1 km. The megaknickpoint zones are not precisely periodic in the model or in nature, but they do have a characteristic distance between them that is approximately several hundred meters to 1 km.

The height of the largest knickpoint, hmax, in the model is controlled by the segment with the smallest value of the erodibility coefficient, Kmin:  
Using the parameters chosen for the example given in Figure 13, equation 7 gives hmax = 20 m. This value is larger than, but broadly consistent with, the largest knickpoints observed along in the Soldier and Molino Canyon study sites, which are ~10 m in height. Equation 7 can also be used to calculate the height of every bedrock step in the model by replacing Kmax with the local value of K. Figure 14 plots the cumulative distribution of bedrock step heights for Molino and Soldier Creeks (extracted from the 1 m/pixel LIDAR DEM; Fig. 14A) along with the same distribution for the model (Fig. 14B). Both the real data and the model closely follow an exponential distribution (i.e., a straight line on log-linear scales) for large step sizes. At smaller step sizes, alluvial backfilling causes the steep spike in the model at low slopes because backfilling prevents any slopes lower than Smin from occurring. In the real data, alluvial backfilling occurs over a range of slopes rather than a single value, hence the spike at low slopes is more spread out compared to that in the model. It is difficult to compare the results of the model and the actual data with respect to the distribution of individual step heights precisely, because the model operates on a 10 m horizontal scale (i.e., Δx = 10 m), while the DEM data have a resolution of 1 m. It is not possible to simply run the model at the same scale as the DEM resolution, however, because the value of Δx in the model is not just the resolution of the model output, it also represents the scale above which there is assumed to be no correlation in K values and below which K values are uniform. In the field, individual sills or light-toned bands exceed several meters in thickness, and field observations suggest that several thick sills can cluster together to form resistant rock units that can be as much as 10 m long laterally along the channel direction in the study area. This observation suggests that Δx = 10 m is an appropriate value for the model, but it is difficult to be more precise. As a result of the difference in model and DEM resolutions, however, some ambiguity will necessarily exist in comparing the model predictions with DEM data with regard to individual step heights. That being said, it should be noted that step heights are the only measure of the model that depends on Δx. All of the other measures of longitudinal profile form we consider in this section (e.g., average slope, power spectrum, distribution of distances between megaknickpoint zones) are independent of Δx, and therefore may be directly compared readily with similar measures extracted from DEM data.
In order to compute the characteristic distance between megaknick-point zones predicted by the model, it is useful to introduce some additional concepts from the study of random walks. One of the fundamental properties of a random walk is that the root-mean-squared displacement of the walk increases with the square root of the horizontal distance, x. Specifically, the size of the average excursion of the walk is quantified using the root-mean-squared displacement, 〈h21/2, given by Gallager (1996) for any random walk as:  
where the brackets denote an average value and Sstd is the standard deviation of the slope at the scale of Δx. In the stream-power model with spatially random erodibility, Sstd can be calculated as:  
Using the parameters chosen for the example given in Figure 12D, equation 9) gives Sstd = 0.220579. The average trend of the model above the minimum slope threshold set by Smin is  
The characteristic spacing, λ c, between megaknickpoint zones is obtained by setting 〈h21/2 equal to hav and solving for x = λc. At scales smaller than this characteristic distance, the random walk variations are, on average, greater than the linear trend in the model. At larger scales, the random walk variations are, on average, smaller than the linear trend. The characteristic spacing between megaknickpoint zones coincides with the crossover between the random walk and deterministic linear trends, and is given by:  
Using the parameters chosen for the example given in Figure 12D, equation 11 gives λc = 0.415 km.
Predictions of the model can be compared in detail to the actual longitudinal profiles observed in the Catalina forerange in two principal ways. First, the power spectrum of a function quantifies the relative variability of that function at multiple scales. For a simple random walk, the power spectrum is a power-law function of wave number, k, with an exponent of −2:  
Given the similarities between the stream-power model with spatially random erodibility and the simple random walk, it is reasonable to expect that the longitudinal profiles of the model and of the actual profiles observed in the Catalina forerange will also follow equation 12. The power spectrum is a powerful tool for analysis because, while random walk–like functions can look very different in the spatial domain, their power spectra all collapse to the same power-law dependence given by equation 12.

Figure 15A plots the power spectrum of the model along with the average power spectrum of the actual profiles of the Catalina forerange. In each case, the linear trend was first removed from both the model and actual data sets. Detrending is a standard preprocessing step in time series analysis and has the advantage that the power spectra will not depend on the magnitude of Sav, which varies depending on the model parameters. Figure 15A plots the power spectra, S(k), as a function of wave number, k, on logarithmic scales for the model (light gray curve, results obtained by averaging the spectra of 1000 independent simulations), for all of the channel profiles shown in Figure 12C (profiles extracted from the 10 m/pixel U.S. Geological Survey DEM) (dark gray curve), and for profiles of Molino and Soldier Canyons extracted from the 1 m/pixel LIDAR DEMs (dark line). Also shown (straight black line) is the result for a simple random walk. The data sets and model closely follow the power spectrum of a simple random walk. This is not surprising considering that the model profiles (given by equation 5) are constructed from a model that is conceptually very similar to a random walk.

In addition to the power spectrum, the statistical distribution of distances between megaknickpoint zones provides another basis of comparison for the model and actual profiles, as well as a means to test the model prediction for λc in equation 11. Equation 11 represents the characteristic distance between megaknickpoint zones. Megaknickpoint zones are not separated by a single value, but instead have a distribution of values. The distribution of distances between knickpoint zones can be calculated using the crossing statistics of a random walk. For a random walk, the distribution of intervals between successive crossings of the origin is given by a normalized exponential distribution (Gallager, 1996):  
where λc is given by equation 11. Figure 15B plots the normalized probability density function, f(λ), of the intervals between megaknickpoint zones for all of the channel profiles shown in Figure 12C (black histogram) and for the model (dark gray histogram). In order to compute the frequency distribution of distances between megaknickpoint zones, it is necessary to define precisely megaknickpoint for the purposes of the analysis. This definition is, to some extent, arbitrary because knickpoints have a range of sizes and distances between them. Nevertheless, we seek a definition of the megaknickpoint zone that recognizes that knickpoints tend to come in clusters that are separated by zones of alluvial storage or backfilling. For the purposes of this analysis, we defined a mega-knickpoint zone in the model and in the real data to be any portion of the profile where the moving average of slope exceeds 0.15 or 15%. Using different definitions (i.e., smaller averaging windows and/or different thresholds for average slope) results in only slight differences in the results. For comparison, the light gray curve represents the exponential distribution (the expected distribution for a random walk process), i.e., f(λ) = λce−λ/λc with λc = 0.5 km. The exponential distribution clearly represents a good fit to the real and model data. The value of λc predicted by equation 11 should not be exactly equal to the value obtained by fitting the exponential distribution to the data in Figure 15B, because the analysis that led to the data included nonunique spatial averaging of the slope and a threshold criterion for the identification of megaknickpoints. Nevertheless, the value of λc illustrated in Figure 15B (i.e., λc = 0.5 km) is in broad agreement with the theoretical value predicted by equation 11 (i.e., λc = 0.415 km).

The stream-power model with spatially random erodibility, despite its simplicity, is capable of reproducing the basic morphological features of channel longitudinal profiles in the Catalina forerange. More broadly, the model provides a starting point for quantifying the structurally controlled complexity that exists in many bedrock (or mixed bedrock-alluvial) channel longitudinal profiles. Structurally controlled knickpoints are common in nature, but they have generally been underemphasized in bedrock erosion studies, partly because of the difficulty of incorporating structural information explicitly into numerical models. The stream-power model with spatially random erodibility provides one approach to including structural heterogeneity into landform evolution models, and it illustrates how small-scale structural heterogeneity can give rise to large-scale variations in longitudinal profile form.

Several caveats of the model should be noted. First, the model explicitly incorporates lateral variations in bedrock erodibility only. In nature, variations in erodibility will also be encountered vertically as channels incise. As such, a more realistic model would include both lateral and vertical variations in bedrock erodibility. Such a model would not achieve a static, steady-state geometry, but instead would reach a dynamic steady-state condition. Such a model would yield longitudinal profiles that exhibit random walk–like behavior for snapshots in time, but knickpoints would shift laterally over time as new heterogeneities are exhumed. This point is particularly clear when considering the nature of alluvial storage in the model and in actual profiles. Zones of alluvial storage in eroding mountain belts are necessarily ephemeral. In the model of this paper, the steady-state solution predicts zones of alluvial storage but makes no explicit prediction regarding how long these zones persist over time. A model that explicitly incorporates vertical changes in rock strength would be capable of addressing the lateral migration of knickpoints and zones of alluvial storage through time. Second, the model assumes that banding is spatially uncorrelated below the resolution of the model (defined by the value of Δx). The fact that the power spectrum of the model matches the power spectrum of real topographic profiles from the Catalina forerange suggests that this assumption is correct to first order. However, no simple statistical model is likely to precisely honor the complexity of small-scale structural elements in real mountain belts. For example, leucogranite sills come in a variety of thicknesses, and field observations suggest that thick sills can occur in clusters. These complexities make it difficult to define a single, unique value for Δx.

Several additional caveats should be noted regarding the application of the stream-power model to channels of the Santa Catalina Mountains. First, the stream-power model is most accurately applied to channels undergoing plucking-dominated erosion. Field evidence clearly indicates, however, that bedrock channels in the study area undergo both plucking-dominated and saltation-abrasion–dominated erosion, with the predominant process type depending on the local rock resistance to erosion. Modeling such process dependence on structure would require a more sophisticated model that includes stream-power–driven erosion, sediment-flux–driven erosion, and the transition between the two (e.g., Gasparini et al., 2007). Structural knickpoints may be enhanced as the dominant process switches from plucking to saltation abrasion within a structurally controlled knickpoint. The separation of flow from the channel bed may also enhance the persistence of knickpoints (Crosby et al., 2007). Also, channels respond to variations in rock type by varying both channel slope and channel width. Specifically, channels both narrow and steepen in response to more resistant rocks (Wohl and Merritt, 2001; Duvall et al., 2004; Amos and Burbank, 2007). Field evidence indicates that channels in the Catalina forerange respond similarly. The model of this paper does not explicitly vary the channel width in response to variations in rock strength or bedrock erodibility. This is an important limitation of the model that should be emphasized. It should also be noted, however, that bedrock erodibility coefficients in the model are empirical coefficients, which cannot (at this time) be calibrated directly based on quantitative measurements of joint spacing and/or rock hardness. Because K values are empirical, they can be defined to implicitly include the effects of channel width adjustment to variations in rock strength. In the absence of channel width variations, zones of unusually resistant bedrock would be characterized by even lower values of K than we assumed because the channel would have to steepen even more that it already does in order to erode through those resistant units. Therefore, by choosing an appropriate range of K values, the effects of channel width adjustment to variations in rock strength can be implicitly included so that the model predictions for variations in channel slope and longitudinal profile form are realistic.


Metamorphic core complexes commonly exhibit drainages that are preferentially oriented parallel and perpendicular to the direction of tectonic extension. In this paper we used field observations, DEM and aerial photographic analyses, and numerical modeling to test hypotheses for the structural and tectonic control of drainage architecture in metamorphic core complexes. Field observations clearly show that channels preferentially exploit steeply dipping joint sets in the Catalina-Rincon core complex. DEM analyses also show that drainage segments are more frequently oriented parallel to joint sets than along other directions. In order to test the joint-controlled hypothesis in more detail, we incorporated joint-controlled bedrock channel incision into a numerical model based on the stream-power model by using an orientation-dependent bedrock erodibility coefficient. The results of this numerical model indicate that joint-controlled bedrock channel incision alone is insufficient for producing the observed pattern of drainage architecture in metamorphic core complexes. When an initial phase of extension-parallel tectonic tilting is also included in the model, the model predicts drainage architectures very similar to those observed.

Field observations indicate that structurally controlled knickpoints in the Catalina forerange are related to the exhumation of unusually thick, felsic-rich leucogranite sills with relatively low joint densities. In order to better understand the relationship between spatial variations in joint density and variations in longitudinal profile form, we modified the numerical model using a spatially random bedrock erodibility coefficient. In topographic steady state, longitudinal profiles predicted by the model can be written as the integral of a random variable. The resulting model is similar to a random walk with a linear trend superimposed. The model illustrates how large-scale (1–10 km) steps in longitudinal profiles can develop from small-scale structural variations. As a stochastic model, the model cannot predict the exact configuration of structurally controlled knickpoints, but it can predict statistically the geometry of structurally controlled knickpoints and zones of alluvial storage. The model is capable of reproducing the basic features of longitudinal profiles of the Catalina forerange, including the power spectrum and the frequency-size distributions of bedrock steps and distances between megaknickpoint zones. This model should form the basis for future studies that seek to quantify how spatial variations in bedrock erodibility influence longitudinal profile form.

We thank the Department of Geosciences, University of Arizona, for financial support of this work through course Geos 650: Field Studies in Geomorphology. We also thank Ellen Wohl and Scott Miller, who provided helpful reviews that greatly improved the manuscript, Evan Canfield of the Pima County Flood Control District, who graciously provided the LIDAR data used in the study, and Paul Kapp, whose review and helpful suggestions greatly improved the manuscript.