Establishing chronologies for alluvial-fan sequences with analysis of high-resolution topographic data: San Luis Valley, Colorado, USA

Where channels emerge from confined catchments onto unconfined valley floors, conical accumulations of sediment termed “alluvial fans” often form, consisting of series of depositional surfaces truncated and inset by progressively younger surfaces that record a history of proximal sedimentation. Recognition of fan surfaces of differing age is critical for appropriately categorizing flood hazards (House, 2005; Pelletier et al., 2005; Lancaster et al., 2012), as fans that were long ago abandoned by active sedimentation should be associated with minimal flood- and debris-flow–related risk. Where they are cut by active faults, alluvial fans serve as an important datum in measuring rates of deformation to quantify seismic hazard. The presence of distinct, inset fan surfaces that appear to record alternating periods of fan construction and entrenchment has led past workers to question the role of climate and tectonics in driving cycles of erosion and deposition on fans, leading to hypotheses that often require a chronology of fan development to be tested (e.g., Bull, 1977; Wells et al., 1987; Abrams and Chadwick, 1994; Ritter et al., 1995; Reheis et al., 1996; Harvey et al., 1999; McDonald et al., 2003; DeLong et al., 2008).

Primary sedimentation processes on alluvial fans, i.e., those that carry sediment from source catchments and construct alluvial fans (e.g., debris flows), produce a rough morphology often composed of boulder berms and levees and lobes (Blair and McPherson, 2009). In the absence of primary, fan-constructing processes, other “secondary” processes reshape alluvial-fan surfaces by reworking or removing the sediment comprising the fan. Blair and McPherson (2009) noted that secondary processes include weathering and soil formation, bioturbation, eolian deposition and winnowing, and reworking by overland flow. In this paper, we take the “age” of an alluvial-fan surface to be the time at which constructional processes on that part of the fan ceased (e.g., due to a change in climate driving incision into the fan) and secondary processes began to dominate fan-surface evolution. Differences in the extent to which secondary processes, particularly weathering and soil development, have reshaped alluvial-fan surfaces have long been utilized to identify fan surfaces of differing ages (e.g., Machette, 1985; Wells et al., 1987; McFadden et al., 1989; Pazzaglia and Wells, 1990; Ritter et al., 1993), and, recently, high-resolution topographic data have enabled the quantification of the general smoothing of fan surface morphology as a function of age (Frankel and Dolan, 2007; Mushkin et al., 2014).

The progressive smoothing of rough topography produced by catastrophic processes (e.g., debris flows deposited on fan surfaces) and the expanding catalog of digital topographic data present an alluring possibility: Can we extract information about the “age” of a landform from topographic analysis alone? Given two ingredients, knowledge of the initial landscape morphology and a description of how those landscapes evolved, it is possible to establish a “morphologic age,” which is defined as a measure of the time since a landform was produced, given a certain efficiency of surface modification. This approach has been used on fault scarps through recognition of the rounding of these features from an initially sharp escarpment produced by surface rupture (e.g., Andrews and Hanks, 1985; Avouac, 1993). Similarly, Hsu and Pelletier (2004) established a morphologic age for fan incision through analysis of the escarpments produced at the margins of incised fan surfaces. The topographic evolution of meteor impact craters has been suggested to be governed by diffusive transport (Soderblom, 1970), an observation that in turn enables quantification of planetary chronologies on remote bodies (Fassett and Thomson, 2014). Correlations between the roughness of geomorphic surfaces and their age have proven useful in reconstructing landscape histories (Cerovski-Darriau et al., 2014). Recognizing the progressive smoothing of landslides with increasing durations of inactivity (Wieczorek, 1984), recent efforts have extended landslide chronologies determined from radiocarbon dating to undated features by comparing simulations of landscape evolution to observations of surface roughness (LaHusen et al., 2016; Booth et al., 2017).

We build upon these past efforts by developing a method to establish a morphologic age for the surfaces of alluvial fans. We describe the evolution of abandoned fan surfaces with a linear, slope-dependent transport rule (i.e., diffusive modification, as described in Culling, 1960). Using a solution for the diffusive modification of arbitrary topography that relies on a transformation of topographic data into frequency space, we directly compute morphologic ages from comparison of spectrally transformed data. While we focus on alluvial fans, the method we present may be applied to any landform for which we can measure a modern analog for the initial topography and for which the evolution is adequately described by linear, slope-dependent transport. In testing this method on alluvial fans, we illustrate how these two criteria can be evaluated from topographic measurements alone.

Using calculated morphologic ages, we developed a calibration between fan surface topography and age on a suite of alluvial fans from the San Luis Valley of Colorado, USA (Fig. 1). We demonstrate that the calculated morphologic ages agree with expectations from stratigraphic and geochronological constraints and suggest efficiencies of surface modification broadly in agreement with past estimates. However, some aspects of the topography of the fans studied here are inconsistent with our model of fan-surface evolution; the oldest fans are likely approaching or beyond the limit of the resolving power of the technique, and systematic variations in morphologic ages point to either problems in our definition of the initial morphology of the fans or our choice of landscape evolution model. While this limits the ability to interpret some of the morphologic ages calculated here, it also illustrates how this method can be adapted to define detection limits for morphologic age determination. Thus, when appropriately applied, this method for morphologic dating of fan surfaces provides an assessment of the history of fan deposition and abandonment that can be calibrated in terms of an absolute age and validated in terms of the model used to derive that age.

The Sangre de Cristo Mountains define the eastern margin of the San Luis Valley, in southern Colorado (Fig. 1). The Blanca massif is located at a prominent step in the range front, with peaks extending to elevations over 4300 m. Many of the valleys that dissect the Blanca massif show evidence of repeat glaciations, and series of inset alluvial fans are present at the outlets of these valleys downstream of glacial moraines (Figs. 1 and 2A). Here, we focus on the largest fans sourced from west-draining, glaciated catchments. Based on the previous geologic mapping of McCalpin (1982) and Thompson et al. (2015) and new geologic mapping conducted herein, we simplified these fan units into three broad characterizations, termed young (Qay), intermediate (Qai), and old (Qao).

Areas of fans near the mouths of valleys are typically characterized by apparently anastomosing channel morphologies with relatively high relief and sharp definition; we refer to these youngest surfaces as Qay. This primary depositional topography, or ridge-and-swale topography, is characterized by cobble- to boulder-mantled surfaces. From field observations, the relief between ridges and swales is on the order of 1 m, and ridges often occur in pairs separated by <10 m and elongated in the transport direction. The most recently active parts of fans have evidence of recent transport: disturbed vegetation, no soil cover, and evidence of scour into vegetated, soil-mantled surfaces (Fig. 3B). These surfaces typically have sparse vegetation and minimal accumulation of fine organic matter or silt, sand, and gravel between boulders (Fig. 3A). Reworking of fan deposits can begin shortly after fan surfaces are abandoned by active transport. This is evident from the observations that some fan surfaces set a few meters or less above these most active surfaces have smoother topography, accumulations of organic material and other fine sediment in the interstices of clasts, and commonly well-developed vegetation (Figs. 3B and 3C). Qay surfaces therefore contain some primary depositional topography and other topography that has been slightly modified by the secondary processes described by Blair and McPherson (2009). However, the ridge-and-swale landscape is still readily identifiable in the field, and individual ridge crests can still be paired with one another. Both the individual ridges and swales, and the relief developed from incision into surfaces with minimal soil mantle produce a rough texture to the fan’s surface that is readily identifiable in a 1-m-resolution, airborne laser swath mapping (ALSM) digital elevation model (DEM; Figs. 1 and 2).

Fan surfaces set farther above active channels show evidence of prolonged modification of primary depositional topography. Immediately above young fan surfaces, we identified “intermediate” surfaces (Qai) as retaining some ridge-and-swale topography. These surfaces display a uniform mantling of silt, sand, gravel, and organic material burying a majority of boulders and cobbles (Figs. 4A and 4B). In these settings, ridge-and-swale topography is more subdued than is evident on young surfaces, to the extent that it is often difficult to trace individual ridges for more than a few ridge widths upslope.

Above intermediate fan surfaces, we identified “old” fan surfaces (Qao) as having minimal to no identifiable primary depositional topography. On these surfaces, a continuous mantling of soil covers nearly all weathered and fragmented cobbles and boulders, such that the internal surface of these parts of fans is almost topographically featureless (Figs. 4C and 4D).

Maps of surface roughness highlight how texture varies on alluvial-fan surfaces (Fig. 2B). Topographic profiles further highlight the subdued variability in surface elevations on Qai and Qao surfaces with respect to young surfaces, particularly at long wavelengths (Figs. 5C and 5D).

Morphologic variability exists within each of these categories, suggesting that there is room for further subdivision of these units in more detailed studies. While this classification scheme shares similarities with other nearby geologic maps (McCalpin, 1982; Thompson et al., 2015), we do not suggest that the ages determined here are representative of similarly named surfaces in those maps. Our goal was to develop a topographic estimate for the timing of abandonment of individual fan surfaces. For this reason, we did not assume that Qai surfaces emanating from one catchment were coeval with those sourced from a neighboring catchment.

We sought a prediction for the evolution of fan surface topography as a function of the time since abandonment, so that we could work this problem in reverse: Given a modified fan surface, how much time has elapsed since it was abandoned? Rather than attempting to capture all the secondary processes that are important on fan surfaces (e.g., Blair and McPherson, 2009), we assumed that local transport from disturbances (e.g., bioturbation, salt growth, and the shrink-swell behavior of clays) would dominate the evolution of fan surfaces. The primary disadvantage of this approach is that it assumes that processes that are known to be important in the evolution of many fan surfaces, such as dust deposition (McFadden et al., 1987), the physical weathering of large clasts (e.g., McFadden et al., 2005; Mushkin et al., 2014), and the downslope transport of reworked sediments by overland flow (de Haas et al., 2014), have a negligible impact on fan surface morphologies. A primary advantage of assuming transport due to local disturbances is that it provides a simple mathematical description of fan evolution that has explicit and testable predictions and that relies on a single parameter (the morphologic age) to describe the smoothing of fan surfaces. In addition, local disturbance–driven transport can reproduce the general smoothing of fan surfaces (Frankel and Dolan, 2007; Mushkin et al., 2014), and it is consistent with differences in cosmogenic nuclide concentrations in ridge-and-swale sediments, highlighting the transport of material from ridges to swales at rates that decrease as fan surfaces are smoothed (Matmon et al., 2006).

Given A_i^t=0 and λ_i, the morphologic age of the surface can be calculated by minimizing the difference between the observed roughness (Eq. 11) and the predicted roughness (Eq. 12).

Because k and t are only present as a product, kt (the morphologic age), in the solution for the diffusion of topography (Eq. 7), we cannot solve directly for either one independent of the other. Equivalent morphologic ages can be associated with surfaces of different absolute age given different efficiencies of surface modification. A best-fitting k can be calibrated from the slope of a plot of morphologic ages as a function of absolute age constraints, and this k value can be used to transform morphologic ages into absolute ages. Recognizing that the morphologic age, which has units of area, is a model of how much diffusive modification has taken place, we will focus on morphologic ages, kt, in discussing computed chronologies.

There are practical and theoretical challenges to computing morphologic ages from comparison of spectral transformations. First, it is impossible to know the true A_i^t=0 for a modified fan surface. Here, we assumed that the morphology of a nearby, active fan substitutes as an effective proxy for the initial morphology of the fan in question. The chaotic nature of active fan surfaces (Figs. 2 and 3) suggests that the details of the surface morphology between active and ancient fans were no doubt different, but it is the typical amplitude of different component waveforms that is of primary importance in applying Equation 10. So, while the particular braided patterns of ridge-and-swale morphology on fans are unlikely to be constant with time, we are primarily interested in ensuring that the dimensions of, for example, ridge-and-swale topography are relatively consistent through time. If the distribution of debris-flow rheologies constructing the alluvial fans controls the roughness of active fans (Whipple and Dunne, 1992), characterizing A_i^t=0 with modern topography is a reasonable starting assumption for catchments of similar size, substrate, and climate.

Diffusion will eventually smooth topography beyond our ability to recognize features that represent primary depositional topography. Equation 10 clearly shows that for small wavelengths, amplitudes decay rapidly with increasing morphologic ages. For expected values of k (10^–4–10^–1 m²/yr; Hurst et al., 2013; Richardson, 2015), after ∼10 k.y., the amplitude of a feature with a wavelength of 3 m will be reduced by at least a factor of 100. In other words, if fan surfaces exhibited perfectly periodic rolling topography with a 3 m wavelength and 1 m amplitude, diffusion would reduce the amplitude of these features to 1 cm or less in 10 k.y. While rapid changes in the amplitude of short-wavelength features would be useful for constraining the evolution of alluvial-fan surfaces deposited after the Last Glacial Maximum, locally as recent as ca. 16 ka (Leonard et al., 2017), they pose a problem for those studies seeking to document older alluvial fans. Given that 1 cm is less than the reported precision for many airborne light detection and ranging (LiDAR) surveys (e.g., those in the U.S. Geological Survey [USGS] 3D Elevation Program; Heideman, 2014), this suggests that morphologic ages will saturate rapidly in the absence of longer-wavelength topographic variability.

Spatial variability in the morphology of active fans (Fig. 2) and partial preservation of older fan surfaces (Figs. 1 and 2) create challenges for defining an appropriate reference initial condition for fans. Past work has recognized down-fan variability in the characteristics of the sediment comprising alluvial fans, and noted the importance of this in establishing comparative criteria for the relative age of fan surfaces (McFadden et al., 1989). These differences may be the result of changes in the rheology of debris flows traversing different parts of the fan, with proximal coarse sediment–rich flows tending to build rougher topography than farther-traveled, coarse sediment–poor flows (Whipple and Dunne, 1992). To calculate a morphologic age, it is important to account for this down-fan decay in roughness when defining the initial roughness ⁠.

In addition, there are challenges related to the application of spectral transformations. First, fans are not infinitely continuous, periodic features. While undisturbed subsets of fan surfaces hundreds of meters across are present in some settings, in other instances, older fan surfaces are crosscut by younger surfaces to form irregular shapes, or they are dissected by drainage development that produces high-amplitude, long-wavelength topography unrelated to the diffusive modification of primary depositional fan morphologies (Fig. 1). Moreover, fans are characterized by their arcuate morphologies. The flows that construct the initially rough surface topography of fans are both responsible for and impacted by this morphology, which causes the dominant directionality of roughness to arc along with the landform.

The random variability in measured elevations (herein “noise”) inherent in even the highest-quality elevation data sets presents an added complication. Because this noise is present across both active and modified surfaces, it impacts the observed power spectra of any age. While the expectation is that the amplitudes of the power spectra will continuously decay from modern geometries (e.g., Eq. 10), noise will create some finite amplitude even on surfaces where any record of the original fan topography was long ago flattened (Fig. 5D).

To apply Equation 12 to establish morphologic ages for fan surfaces, we first identified the spatial extent of the best-preserved fan surfaces (Fig. 1). We defined the margins of fan surfaces using the truncation of relatively smooth surfaces by inset, rougher surfaces identifiable in 1-m-resolution gridded topographic data (Figs. 1 and 2). In an effort to honor the assumptions of Equation 2, we attempted to map around the margins of channels incised into alluvial-fan surfaces, including small tributaries that merged downstream into larger, incised valleys. This led to the embayed morphology evident in the mapped extent of the intermediate and older fan surfaces (Figs. 1 and 2). We buffered the margin on fan surfaces inward from where fan surfaces were truncated by incised valleys or younger fans to exclude convex margins that were likely the product of downslope transport.

Within mapped fan surfaces, we measured the amplitudes and wavelengths (A_i and λ_i of Eq. 7) derived from applying the type 2 discrete cosine transform (DCT; Press, 2007) to many topographic profiles collected within fan boundaries. The convolution solution to linear diffusion (Eq. 8), and thus the ability to compare “fresh” and “diffused” power spectra to estimate morphologic ages is easily transformed to two dimensions, allowing us to instead conduct this analysis with moving windows. We have chosen profiles here because they capture the primary roughness recognized on young fans, are simple to orient in the dominant roughness direction, and, as we show in benchmarking experiments, are able to reproduce expected morphologic ages. Topographic profiles were 150 m long and oriented approximately perpendicular to slope. To smooth out local topographic gradients (e.g., on the margins of ridges or swales), we measured the orientation of the slope here broadly, using a 40 m × 40 m centered difference kernel. We ignored any profiles that extended off the mapped margin of individual fan units and limited our analysis to wavelengths shorter than our profiles and larger than the grid spacing (i.e., λ < 150 m and λ > 1 m). Within each fan, we limited the analysis to 10,000 profiles randomly distributed throughout the fan and calculated a morphologic age for each profile. The Qao surface of the Zapata complex was an exception, because only 3129 profiles within the masked region were available. The choice of 10,000 profiles was arbitrary and was simply intended to give an accurate depiction of the topographic variability within the fan. Prior to performing the DCT on each profile, the best-fitting parabola was subtracted from the fan profile to remove the mean elevation and the long-wavelength variation that come from collecting a straight-line profile on a roughly conical landform.

To calculate a morphologic age, it is important to account for the observed variability in the morphology of active fans. We normalized for down-fan variability in older surfaces by defining a relationship between DCT amplitudes and along-fan distance based on analysis of an active, reference fan. For each wavelength of interest in the DCT, we performed a regression between the magnitude of the DCT amplitudes and the distance down fan (Fig. 6). Down-fan distance was measured from the point where the largest channel supplying an alluvial fan exited the confines of its mountain valley (apex point), and we used the same outlet position for each unit within a fan complex (Fig. 1). At the distance associated with each topographic profile, we calculated from predictions determined from the characterization of a reference fan. Specifically, we characterized the change in initial amplitudes with distance using a linear fit to observations from the Qay unit of the Pioneer complex (Fig. 6). Our choice of a linear fit was arbitrary, but it captures the general decay in amplitudes present in the data. While we defined a Qay unit at each fan complex, the Qay deposits in the Pioneer fan appeared to be the freshest, with expansive areas of bare cobbles and boulders (Figs. 3A and 3B). In contrast, large parts of the Zapata Qay unit appeared to have undergone some modification (Fig. 3C). We computed morphologic ages for all mapped fan surfaces. This included the Qay units from the Holbrook and Zapata complexes, where we expected morphologic ages to be near zero, but possibly higher due to regions of modified depositional topography (e.g., Fig. 3C). Additionally, we calculated a morphologic age on the Pioneer Qay unit. By definition, the mean morphologic age of the Pioneer Qay unit will be zero, but it is nonetheless useful to analyze because it provides an opportunity to assess the variability in morphologic ages within a unit due to variations in roughness within a fan unit (Fig. 2B).

To characterize L_n (noise threshold in Eq. 11), we followed past work quantifying surface roughness from spectral transformations (Brodsky et al., 2011; Mushkin et al., 2014). We measured L_n as the wavelength at an inflection in plots of A as a function of λ for each of our measured profiles (Fig. 7). Specifically, we sought the wavelength that divided two piecewise power-law regressions where the sum of the squared residuals between the log-transformed predictions and observations is minimized.

We tested the validity of this morphologic age approach with a set of benchmarking experiments on topography that has been forward-modeled with a two-dimensional version of Equation 2. We used a subset of the Qay unit of the Zapata fan complex as an initial condition (Fig. 1), because the broad extent of this unit in proximal regions allowed us to clip a large rectangular area in the region of the fan that displayed the greatest roughness (Fig. 2B). To ensure that reference topographies used for benchmarking our morphologic age determinations were produced in a manner unique from our morphologic age methodology, we relied on forward models of topographic evolution rather than the solution for topographic diffusion presented in Equation 7 and 8 herein. Forward models are integrated for varying morphologic ages using a forward-Euler approach, where the second derivative of elevation is approximated with a second-order, centered difference approach, and the time step is limited by the von Neumann stability criteria ⁠. We compared predicted morphologic ages to the expected morphologic ages derived from the forward model duration for varying profile lengths. To account for the influence of noise on measurements of topography, we added random, Gaussian noise with 0 mean and 10 cm standard deviation before making each calculation of kt values. A value of 10 cm was chosen because it was comparable to one of the quality criteria for “QL2” LiDAR, the most widely used standard in the USGS 3D Elevation Program (3DEP; Heideman, 2014). We present our benchmarking tests with and without this added noise.

We placed independent constraints on the age of fan surfaces in the study area using U-Th series dating of pedogenic carbonate at sites on the Qai (site ZP-01) and Qao (site ZP-02) surfaces of the Zapata fan (Fig. 1). U-Th series dating has previously provided reasonable ages for pedogenic carbonates to constrain surface stabilization in similar environments (Ludwig and Paces, 2002; Sharp et al., 2003). These ages are interpreted based on the inference that carbonate formation is synchronous with or postdates soil development, which would in turn postdate the last transport events that constructed the fan surface (i.e., a minimum-limiting age). We assumed that the variability in age within a single map unit was small enough that we could characterize the age of a broad fan surface with a set of samples from a single location. We sampled 8–10 individual clasts from >50 cm below the fan surface in alluvial gravels exposed either by drainage incision into the Qai surface, or in a hand-dug soil pit on the Qao surface. Sampled clasts had dense carbonate undercoatings (rinds), which were then cemented by finer matrix carbonate within a clear carbonate horizon within the soil profile (Supplemental Fig. 1A¹). Carbonate powders were sampled from cross-sectioned, polished carbonate undercoatings from each clast (supplemental Fig. 1B [see ^{footnote 1}]) and measured by thermal ionization mass spectrometry at the USGS Southwest Isotope Research Laboratory following established methods (see Supplemental Methods for full details [^{footnote 1}]).

The distributions of predicted morphologic ages (kt_pred) versus expected morphologic ages (kt_exp) from 1–500 m², using 50-, 100-, and 200-m-long profiles, are shown in Figure 8. The two panels show predicted morphologic ages with (Fig. 8B) and without (Fig. 8A) a 10 cm component of random noise. For a given expected morphologic age, there is considerable variation in predicted morphologic ages. However, the distributions of kt_pred are generally well centered about the expected values for low to intermediate morphologic ages. However, as the expected morphologic age increases, predicted morphologic ages begin to taper off, yielding a near-constant, median-predicted value for increasing expected morphologic age. In these experiments, the observed saturation in morphologic ages occurred at ∼10² m², which, for a value of k in the middle of the expected range (∼0.003 m²/yr), corresponds to an age of ∼30 k.y.

For low expected morphologic ages, there is considerable overlap in predicted morphologic ages for the range of profile lengths assessed here (Fig. 8). Differences between predicted morphologic ages emerge for varying profile lengths, as predicted morphologic ages begin to saturate at lower expected morphologic ages when using shorter profiles. Profile length also determines the maximum wavelengths resolvable by the DCT.

When noise is present, predicted morphologic ages are commonly lower than expected morphologic ages, particularly when using shorter profiles, despite our efforts to mask out small-wavelength, noise-dominated data (Fig. 7), although there is considerable overlap between predicted morphologic ages and the expected value, even with noise, particularly when measuring topography with longer profiles (Fig. 8B). The addition of noise also causes the apparent saturation in morphologic ages to occur at lower morphologic ages (Fig. 8B).

Measured DCT amplitudes across fan surfaces tend to increase log-linearly with wavelengths on both active and older fan surfaces. The magnitude of the amplitudes tend to decrease on older fans, similar to observations from landslides of increasing age (Booth et al., 2017). There is little variation in the slope of amplitudes as a function of wavelength in log-space, contrary to the expectation that diffusion should more quickly reduce the amplitude of short-wavelength topography (Fig. 9B). Figure 9 also shows the expected decay of amplitudes as a function of wavelength and morphologic age (i.e., Eq. 10) given an initial condition described by the best-fit power-law relationship shown in Figure 9A (defined where λ > L_n).

Median morphologic ages for fan surfaces ranged from 1 to 60 m², and distributions of morphologic ages observed on individual fan surfaces are shown in Figures 10 and 11 and summarized in Table 1, where we report the lower 5th and upper 95th percentiles of observed morphologic ages as a measure of uncertainty. In some cases, individual profiles had larger DCT amplitudes than were observed in the mean condition of the reference fan. These required negative morphologic ages in our formulation (Eq. 12; e.g., Holbrook young in Table 1), which we excluded when constructing the distributions in Figure 10 to facilitate plotting in log space. Distributions of kt have long tails, some of which include significant fractions of data (e.g., Zapata young and Holbrook young), and some of which occasionally show multimodal distributions (e.g., Zapata old).

We utilized two tests to assess how well observed morphologic ages met expectations derived from our choice of initial conditions and our use of linear diffusion to describe fan smoothing. First, we tested whether calculated morphologic ages corresponded to the expected stratigraphic order defined by map relations. We compared the distribution of morphologic ages observed in a fan unit to the distribution observed in the next-youngest unit within the same fan complex; that is, we tested that morphologic ages of the Qao surface of Zapata were older than the Qai surface of Zapata, and that the Qai surface of Zapata were older than the Qay surface of Zapata. For young fan units (Qay, Fig. 1), we compared morphologic ages to the reference fan that defined our initial condition (the Qay unit at Pioneer; Fig. 1). Specifically, we applied the Mann-Whitney U test to pairs of data from mapped regions, testing the null hypothesis that a randomly selected morphologic age from one fan is equally likely to be smaller or larger than a randomly selected morphologic age from another fan (Mann and Whitney, 1947). Rejecting the null hypothesis, which we do for p < 10^–3, requires accepting the alternative, that randomly selected morphologic ages from one fan are likely to be larger than those from another. All the distributions of morphologic ages obeyed this basic stratigraphic test (p << 10^–3; Table 1), with the exception of the old unit of the Pioneer fan, the morphologic age of which was not greater than the intermediate Pioneer fan unit.

Second, we tested whether morphologic ages were consistent along the fan in plots of kt as a function of distance (Fig. 12). Down-fan variability in morphologic ages was observed, which we summarized in Table 1 by reporting the best-fitting slope to linear regression and a p-value, p(kt slope), for the null hypothesis that the slope is zero. We considered down-fan variability in morphologic ages to be significant when p(kt slope) < 10^–3 (Table 1).

The minimum limiting surface ages provided by U-Th series dating allow independent constraints to be placed on the efficiency of surface modification on the studied fans that can be compared to previous work. We generated a total of 19 U-Th series ages on 16 clasts from the two sampling sites at the Zapata fan complex, which ranged from 19.8 ± 0.1 to 26.2 ± 1.2 ka for the Qai surface soil, and 19.0 ± 0.2 to 115.1 ± 1.2 ka for the Qao surface soil (errors are 2σ). One additional older age obtained from the Qao surface was rejected due to suspected open-system uranium loss (see Supplemental material [^{footnote 1}]). Samples from both sites showed a range of ages, which we interpreted to reflect sampling mixtures of carbonate formed semicontinuously between the time of soil formation and the present. Given this, we adopted the oldest credible age from each site to provide the closest constraint on surface stabilization (Table 1).

Figure 13 highlights the range in possible diffusion coefficients that could explain our observations of morphologic ages given available geochronologic constraints (k < 2 × 10^–3 m²). This range is simply the bound on values of k that would predict 95% of kt values observed within dated fans and predict ages of surface abandonment on U-series–dated surfaces older than the lower 2σ uncertainty bound. To simplify discussion, we also calculated a single preferred value for the diffusion coefficient, k_bf = 8 × 10^–4 m², where the subscript bf indicates best fit. To calculate k_bf, we assumed that the age and analytical uncertainty of U-series ages effectively characterized the age of fan-surface abandonment and the onset of diffusive modification. We then determined k_bf as the value that minimized the sum of the squared errors between the observations of kt within fans and the product of k_bf and the U-series age of the surface. Calibrated morphologic ages for fan surfaces are summarized in Table 1 and Figure 11.

Morphologic dating has long been useful in studies of fault scarps and scarp-like landforms (e.g., Andrews and Hanks, 1985; Hanks, 2000) and shows promise for other geologic hazards (e.g., LaHusen et al., 2016; Booth et al., 2017). In our benchmarking experiments, we show that we are able to effectively extract young morphologic ages when linear diffusion is responsible for topographic smoothing and the initial conditions are known (Fig. 8). However, these benchmarking experiments highlight the finding that internal fan variability in original morphology can lead to a large spread in calculated morphologic age, even in the absence of noise (Fig. 8A). Noise will produce persistent roughness in the topography at all morphologic ages, biasing morphologic ages toward smaller values, unless its contribution to roughness measurements is minimized (e.g., by ignoring those wavelengths dominated by noise; Fig. 7).

In applying our morphologic dating method to real examples, calculated morphologic ages showed even greater variability than observed in benchmarking experiments (Fig. 10; Table 1), with significant portions of morphologic ages for a single fan spanning orders of magnitude, indicating that any single morphologic age assigned to these fans would be associated with a large uncertainty. Observed distributions of kt revealed generally increasing morphologic ages with unit age (Figs. 10 and 11). This provides confidence that, at the least, our method of morphologic age determination provides an objective criterion for assigning relative ages to fan surfaces based on their topography alone.

Our best-fitting diffusivity of 8 × 10^–4 m² yr^–1 is within the range of diffusivity values estimated by past work (Hurst et al., 2013; Richardson, 2015), although it lies at the low extreme of these estimates, and a broad range of diffusivities is allowed by the data (Fig. 13). Perhaps this low diffusivity is the consequence of the mean climate in the San Luis Valley, as Richardson (2015) noted a positive correlation between diffusivities and the aridity index (Zomer et al., 2008), suggesting that more arid climates have lower diffusivities. Alternatively, as discussed below, the morphologic ages of the older fan surfaces may be saturated and suppressed, as benchmarking experiments show that calculated morphologic ages on these alluvial fans saturate at ∼10² m² (Fig. 8). Following this, the upper confidence bound on diffusivity (2 × 10^–3 m² yr^–1; Fig. 13) may be a better estimate and is closer to commonly observed values from past work (Hurst et al., 2013; Richardson, 2015).

Beyond the uncertainty in the modeled morphologic age and calibrated diffusivity, there is further uncertainty in whether the model itself is applicable. In the sections that follow, we attempt to address this latter question.

To assign a morphologic age, it is necessary to assume an initial topography. Here, we utilized active modern fan surfaces as a proxy for the initial conditions of old fan surfaces. Because active fans are smoother in their distal reaches than proximal to their source catchments (Fig. 2B), we assumed that the decay in roughness observed in an active fan (Fig. 6) may be applied when characterizing the initial roughness of older, modified fans. If this assumption is valid, and fan evolution could be perfectly described by linear diffusion, kt values would be constant moving down fans; however, that is not generally observed (Fig. 12; Table 1). Down-fan variability in kt was seen in all fans, except the youngest two units from the Holbrook complex (Table 1). Again, this provides a check on estimated morphologic ages that depends exclusively on topographic analysis. The down-slope variability in kt may arise from an inappropriate characterization of the initial conditions in fan surfaces; that is, the initial roughness of older fan surfaces may have been different from what is observed today. Differences in initial roughness could result from variability in the processes and properties (in space or time) of flows carrying sediment from source catchments and constructing fans (Whipple and Dunne, 1992). Alternatively, as we discuss below, postdepositional processes that modify fans may also be responsible for down-fan variability in kt.

Linear diffusion predicts that short-wavelength, high-frequency elements of the topography (e.g., individual debris-flow levees; Fig. 3) will be smoothed and flattened more quickly than long-wavelength landscape elements (e.g., avulsed and incised channels; Fig. 3). As a result of this effect (Eq. 10), we would expect the power-law scaling of amplitude as a function of wavelength to steepen with time as diffusion draws down the amplitudes of short wavelengths (e.g., Fig. 9A). However, that is not generally observed (Fig. 9B). The exponents of power-law regressions are relatively similar across fan units, similar to past findings on alluvial fans (Mushkin et al., 2014). The ability to recognize that the data (in this case, amplitude-wavelength relationships) do not fit the model (the wavelength-dependent decay in amplitudes predicted by linear diffusion) without independent geochronologic control is a strength of the approach to roughness dating we present here.

In support of linear diffusion, we observed processes that we would expect to produce diffusive modification of surfaces (e.g., tree-throw; Fig. 14), and the Blanca massif fans do generally get smoother, such that morphologic ages are consistent with map relationships (Figs. 10 and 11), and rates of smoothing (Fig. 13) are broadly consistent with expectations (Hurst et al., 2013; Richardson, 2015). However, as discussed above, our simple model for diffusive modification (Eq. 2) ignores a host of processes important for fan surface smoothing (e.g., dust deposition and overland flow). Indeed, some elements of our analysis of fan surface topography suggest that aspects of the evolution of these fan surfaces are missed by the assumptions and simplifications we made in order to apply Equations 2 and 10 to understand fan surface age. This is highlighted by the fractions of fans that failed at least one of our tests of consistency with expectations (Fig. 11; Table 1). Other efforts to morphologically date landforms may be less concerned about the impacts of dust input or transport by overland flow; however, we would still encourage future efforts to assess the validity of their initial conditions and characterization of diffusive transport. Spectral transformation provides one way of accomplishing that, because of the predictions of linear diffusion (Fig. 9), as has been shown for other hillslope transport rules (Doane et al., 2018).

Similarly, our simplification of this problem to one-dimensional diffusion in the across-fan direction ignores diffusive modification in the down-fan direction. However, given that this simplification accurately recovers morphologic ages in the young, benchmarked cases (where diffusion was forward modeled in two dimensions; Fig. 8), we believe this assumption is unlikely to have significantly limited our results and address likely complications in the sections that follow.

At the wavelength of our 150 m profiles, morphologic half-lives are ∼400 m², although for the longest wavelength we examined (100 m), kt_1/2 is 175 m². Approaching the wavelength of individual debris-flow deposits (i.e., 3–10 m), morphologic half-lives would be as small as 0.2–2 m². Thus, many of the topographic features we associate with active fans have experienced many half-lives in the time represented by the sets of fans studied here. While the longest wavelengths that we could capture may not have decayed to half of their original amplitude, the integrated impact of decaying the suite of wavelengths we observed seems to result in saturation of morphologic ages at ∼10² m² (Fig. 8). The actual saturation limit for morphologic age determination of different landforms will depend on the amplitudes and wavelengths of topographic elements comprising those features and the level of noise present in the elevation data. We suggest that future studies quantify the limits to their morphologic ages imposed by the spatial scale over which roughness is being measured and the initial topography being modified (e.g., Fig. 8). Given that morphologic ages of the oldest fans around Blanca massif have likely reached saturation and are therefore artificially low, they should be interpreted with some caution. This is not surprising in light of the fact that there is little-to-no observable topography visible on these surfaces in the field (Figs. 4C and 4D).

The application of Equation 2 implies a continuum behavior to land-surface change. This approach may break down at length scales that approach the displacement distance of random disturbances. Jyotsna and Haff (1997) noted that large-scale disturbances (e.g., slumps and the burrows of large animals) would actually produce roughness that would subsequently be smoothed by smaller-scale disturbances (e.g., rainsplash, shrink-swell processes, and smaller burrows). Given that processes of different scales occur at different rates, the relative concentration of larger-scale roughness elements can be used to estimate the efficacy of diffusive transport by a short-length-scale process (Jyotsna and Haff, 1997). Similarly, in our approach, length-scale–dependent diffusive processes could act to modify the component waveforms at different rates, making k wavelength-dependent.

If, as concluded by Jyotsna and Haff (1997), the majority of transport is accomplished by short-length-scale processes (i.e., short with respect to the wavelength of primary alluvial-fan morphology), then assuming a constant k across all wavelengths should still appropriately characterize the diffusive modification of fan surfaces. At least in the present climate, most of these fan surfaces are mantled by trees and shrubs. If the growth of vegetative root systems and the death, felling, and resultant uprooting of plants were responsible for a large portion of diffusive transport (Fig. 14A), as is likely the case in other well-vegetated regions (Roering et al., 2010), most of these disturbances likely would result in transport distances on the order of meters or less, i.e., length scales approaching or smaller than those for which we measured roughness (Fig. 7).

Past variations in climate likely impacted the efficacy of mass-transport processes such that k may also vary through time (Marshall et al., 2015; Booth et al., 2017). As highlighted by Doane et al. (2018), amplitude ratios would still follow the general decay predicted in Equations 10 and 13, but morphologic ages would not increase linearly with actual age. Demonstrating this effect would require a suite of calibrated alluvial-fan ages, to which we did not have access in this study. Microclimatic gradients, such as those that occur on slopes with different aspects, have also been shown to impact the efficiency of diffusion, and hence morphologic ages, on fault scarps (Pierce and Colman, 1987). We ignored this effect, but we also limited our analyses to west-facing fans, and therefore we do not suspect that this factor hampered our analysis. While our west-facing fans do contain north- and south-facing aspects, the aspect dependence of diffusion on fault scarps likely arises because these scarps have sufficient slopes and relief to shade the scarp face at certain times of day. For the fans studied here, the relatively minimal amplitude of primary fan topography (i.e., the ridge-and-swale topography on the order of ∼1 m tall) offers minimal shading on the surface of the fan. Moreover, if tree throw is an important diffusive process on these fans, as the density of trees in modern time suggests, the tree species present typically exceed the height of the primary fan topography (Figs. 3A, 3B, and 14A), and they are therefore unlikely to be significantly impacted by differences in topographic aspect.

From field investigations of the fans around Blanca massif, the two most obvious contributors to nondiffusive sediment transport are wind and overland flow. Great Sand Dunes National Park is located immediately northwest of the extent shown in Figure 1 (Fig. 14D), and some barchan dunes are even visible at the northwestern edge of the fans studied here (Figs. 1 and 11). If the fans studied here follow the patterns observed at other sites proximal to playa dust sources, such as the former basin of Lake Alamosa immediately west of these fans (Machette et al., 2013), then a steady or climate-modulated rate of dust input may drive dust deposition (Pelletier, 2007). As a consequence, while we do not document this here, it is likely that our fan surfaces of different ages record differences in average rates of dust accumulation (Reheis et al., 1995). The net and variable input of material by a process other than local transport invalidates the mass conservation assumption used to transform Equation 1 to Equation 2. However, in assessing the importance of dust accumulation on smoothing effects of inactive fan surfaces, the issue at hand is the way in which a net accumulation of dust influences surface morphology. If dust preferentially accumulates in low-lying areas, then diffusion alone would underpredict the rate of fan-surface smoothing. Given that our calibrated diffusivity (Fig. 13) is on the low end of previously observed values (Hurst et al., 2013; Richardson, 2015), infilling of topographic lows by dust is unlikely to be the most significant contributor to nondiffusive modification of fan surfaces, despite its overall importance to pedogenesis.

The deep incision into the surfaces of some fans necessitates that at some scales, erosion by overland flow is capable of impacting fan morphology (Fig. 1). While we attempted to mask out these regions in our delineation of fan units, our ability to perform that masking was limited by the scale at which we mapped fan units. Field observations on intermediate and older surfaces reveal fine sediment impounded behind roots (Fig. 14B) and small, tributary channels extending up from the base of fans (Figs. 1, 2, and 14C). These observations suggest that sediment transport and erosion by overland flow occur and extend to finer scales than we were able exclude in our mapping. Blair and McPherson (2009) recognized reworking by overland flow as the most frequent secondary process modifying fan surfaces, and even in hyperarid environments, overland flow is thought to reshape fans by transporting some coarse sediment from ridges farther down the fan (de Haas et al., 2014).

The relief between the thalweg of incised channels and the surface of the fan tends to decrease up fan (Figs. 1 and 2). Therefore, one consequence of surface modification by overland flow on morphologic ages could be a decrease in morphologic age down fan, as lower parts of the fan are more likely to experience greater incision, producing roughness on the fan surface that would yield a younger morphologic age in the framework of Equation 12 and perhaps contributing to the persistence of high-frequency topography on old fan surfaces (Fig. 9). This would tend to draw average morphologic ages down, producing a relatively small calibrated diffusivity like we observed, and it would tend to create a downstream decrease in morphologic age (Fig. 12). To a similar end, given that distal portions of alluvial fans tend to be smoother (Fig. 6), any overland flow erosion at positions farther down the fan will produce a greater relative impact on surface roughness. In line with this expectation, all of our oldest mapped fan units displayed significant negative correlations between down-fan distance and morphologic age (Table 1). While this phenomenon could alternatively arise from a difference between the decay in roughness (Fig. 6) of the reference fan and the initial conditions of these older fans, our field observations and the down-fan variability of morphologic ages support the inference that our morphologic ages have likely been impacted, and artificially depressed, by overland-flow–driven sediment transport.

Expanding access to detailed digital topographic data makes the idea of morphologic dating, i.e., establishing chronologies based on comparison between observed topography and a model for landscape development, increasingly appealing. Morphologic age determination has long been applied to fault scarps (e.g., Andrews and Hanks, 1985; Avouac, 1993; Hanks, 2000) and more recently has been extended to landslides (LaHusen et al., 2016; Booth et al., 2017). Here, we developed an approach for the morphologic dating of alluvial fans, although the method is extendable to other complicated landforms, and we provided a means of assessing some of the assumptions inherent in morphologic age determination.

As fluvial and debris-flow networks migrate across and sometimes incise into alluvial-fan surfaces, regions of fans become stranded and tend to become smoother with time (Frankel and Dolan, 2007). Linear, slope-dependent transport (e.g., Culling, 1960) can accomplish this smoothing, and although it ignores some important contributions to fan surface evolution, we utilized the predictions of this transport rule to develop a calibration between surface roughness and age on a suite of alluvial fans from the San Luis Valley, Colorado. Our morphologic ages largely honor the constraints imposed by crosscutting relationships (Fig. 10). Limited geochronologic constraints on fan surfaces further indicated a value for the diffusive transport coefficient, k_bf = 8 × 10^–4 (Fig. 13), that is in the relatively broad range of previously calibrated values (e.g., Hurst et al., 2013; Richardson, 2015).

While these broad-scale observations lend credence to our approach, some of the details of the actual observations point to complications. First, we demonstrate that the oldest fans studied here (ca. 100 ka) are at or beyond the limit of the resolving power of this technique (Fig. 8). Second, high-frequency, short-wavelength topography is more persistent on old fan surfaces than diffusion alone would predict (Fig. 9). Finally, morphologic ages show consistent variability with distance from the range front on the oldest fan surfaces, where a single value is expected (Fig 12). Down-fan variability highlights problems with either our choice of initial conditions or our description of fan evolution, casting some doubt on the morphologic ages assigned to these fans, but the ability to recognize such complications points to a strength of this method. One of the strengths of morphologic age dating is the ability to establish objective, quantitative timing constraints from analysis of topography alone. We believe that by formalizing the expectations of topographic changes as a function of morphologic age (e.g., Figs. 9 and 12), we have provided some internal consistency checks for the assumptions inherent in morphologic dating of alluvial fans or other initially rough landforms. This will aid future efforts that require confidence in established chronologies, such as those that wish to use surface morphologies to directly inform hazard assessments or to guide more detailed studies of hazard-prone areas, or those wishing to analyze remotely sensed data from other bodies in the solar system.

We thank G.E. Hilley for introduction to the solution to the diffusion equation used here, and J.P. Perkins for discussions that led to the recognition that this solution could be used for morphologic dating. Insightful reviews from Corina Cerovski-Darriau, Jonathan Harvey, and an anonymous reviewer and suggestions from Editor Shanaka de Silva greatly improved this manuscript. Support for S. Nicovich was partially provided by the U.S. Geological Survey (USGS) EDMAP program and Montana Space Grant Consortium, and support for R.M. Sare was provided by the Stanford-USGS Fellowship and National Earthquake Hazards Reduction Program grant G17AP00010. Primary support for this work came from the FEDMAP component of the USGS National Cooperative Geologic Mapping Program. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.