www.geosci-model-dev.net/3/293/2010/ © Author(s) 2010. This work is distributed under the Creative Commons Attribution 3.0 License.

Abstract. Geomorphic process modeling allows us to evaluate different methods for estimating moraine ages from cosmogenic exposure dates, and may provide a means to identify the processes responsible for the excess scatter among exposure dates on individual moraines. Cosmogenic exposure dating is an elegant method for estimating the ages of moraines, but individual exposure dates are sometimes biased by geomorphic processes. Because exposure dates may be either "too young" or "too old," there are a variety of methods for estimating the ages of moraines from exposure dates. In this paper, we present Monte Carlo-based models of moraine degradation and inheritance of cosmogenic nuclides, and we use the models to examine the effectiveness of these methods. The models estimate the statistical distributions of exposure dates that we would expect to obtain from single moraines, given reasonable geomorphic assumptions. The model of moraine degradation is based on prior examples, but the inheritance model is novel. The statistical distributions of exposure dates from the moraine degradation model are skewed toward young values; in contrast, the statistical distributions of exposure dates from the inheritance model are skewed toward old values. Sensitivity analysis shows that this difference is robust for reasonable parameter choices. Thus, the skewness can help indicate whether a particular data set has problems with inheritance or moraine degradation. Given representative distributions from these two models, we can determine which methods of estimating moraine ages are most successful in recovering the correct age for test cases where this value is known. The mean is a poor estimator of moraine age for data sets drawn from skewed parent distributions, and excluding outliers before calculating the mean does not improve this mismatch. The extreme estimators (youngest date and oldest date) perform well under specific circumstances, but fail in other cases. We suggest a simple estimator that uses the skewnesses of individual data sets to determine whether the youngest date, mean, or oldest date will provide the best estimate of moraine age. Although this method is perhaps the most globally robust of the estimators we tested, it sometimes fails spectacularly. The failure of simple methods to provide accurate estimates of moraine age points toward a need for more sophisticated statistical treatments.


Introduction
Cosmogenic exposure dating is an important technique for learning about glacier size changes during the last ∼10 5 yr of geologic time (Gosse and Phillips, 2001).Glaciers and ice sheets grow and shrink in response to climate change (Dyurgerov and Meier, 2000;Oerlemans, 2005;Jansen et al., 2007).Therefore, reconstructions of past glacier sizes over time yield information on past climates and rates of sea level rise.As glaciers advance and retreat, they mark their former margins with ridges of debris, called moraines (Gibbons et al., 1984).In cosmogenic exposure dating, field geomorphologists collect samples from boulders on the crests of moraines, and the concentrations of certain rare chemical species (cosmogenic nuclides) are measured in the samples.These cosmogenic nuclides are produced at predictable rates in surface materials by cosmic rays (Lal, 1991;Gosse and Phillips, 2001).Under ideal conditions, the ages of the moraines can be calculated directly from the nuclide concentrations (e.g., Gosse et al., 1995a).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Unfortunately, geomorphic processes bias cosmogenic exposure dates (see review in Ivy-Ochs et al., 2007).If the boulders contain some preexisting concentration of cosmogenic nuclides when they are deposited on the moraine, then the exposure dates will tend to overestimate the moraine's age.Most other processes tend to reduce the apparent exposure times of the boulders.For example, cover by snow or sediment reduces the flux of cosmic rays through the upper surfaces of the boulders.The exposure dates from these shielded boulders will underestimate the true age of the moraine on which they rest.Similarly, erosion of boulders removes the most nuclide-rich part of the rocks (Lal, 1991); therefore, eroded boulders also yield exposure dates that underestimate the age of their host moraine.
The effects of these processes on the distributions of exposure dates from moraines are not known a priori, and this uncertainty is reflected in the variety of procedures for estimating the ages of moraines that are described in the literature.Many workers prefer to use some measure of the central tendency of a data set; such estimators include the arithmetic average, the mean weighted by the inverse variance, and the mode (e.g., Kaplan et al., 2005;Licciardi et al., 2004;Kelly et al., 2008).Other investigators prefer extreme estimators, including both the youngest and the oldest dates (e.g., Benson et al., 2005;Briner et al., 2005).For data sets with large ranges, the choice of estimator has a profound effect on the estimated ages of the moraines (for example, compare Chevalier et al., 2005, with Brown et al., 2005).The choice of estimator is typically informed by geomorphic observations.However, without knowledge of the underlying parent distribution from which the dates are drawn, we cannot evaluate the effectiveness of these different procedures.
We might evaluate the effects of geomorphic processes on cosmogenic exposure dating by performing a positive control experiment.In such an experiment, we would identify a moraine whose age was known independently, perhaps from bracketing radiocarbon dates.We would then collect many samples from this moraine for cosmogenic exposure dating, and compare a histogram of the exposure dates to the independently known age of the moraine.The distribution of the exposure dates about the true age of the moraine would tell us the effects of geomorphology on the exposure dates from that moraine, other factors being equal.
Unfortunately, such a positive control experiment is impractical.To achieve robust results, we would need many samples from one moraine.The exact number of samples required is poorly defined, but it seems likely that 50 samples are insufficient (see Murphy, 1964, his Fig. 6).Because cosmogenic exposure dates are expensive, the necessary number of samples is probably not achievable.In addition, the geomorphic processes that affect exposure dating are likely to be highly variable between field sites.Thus, we would need to repeat the experiment on a large sample of moraines, multiplying the cost many times.Moreover, there are few sites where the ages of moraines are known independently, and these sites are already included in the nuclide production rate calibration database (Balco et al., 2008).Last, there are potential confounding effects.The difference between the independently determined age of a moraine and any individual exposure date is influenced by errors in estimating both the age of the moraine and the local production rates of cosmogenic nuclides.These errors interfere with our ability to separate out the effects of geomorphology on exposure dating.Thus, a positive control experiment to isolate the effects of geomorphic processes on exposure dating is prohibitively expensive, probably cannot be done for a representative sample of moraines, and is subject to strong confounding effects from uncertainties in moraine age estimates and nuclide production rates.
Monte Carlo-based numerical models offer a means of assessing the effects of geomorphic processes on cosmogenic exposure dating that avoids the disadvantages of positive control experiments.Although these models can never replace field observations, they provide a test bed for understanding existing exposure dates.Such models can generate thousands of synthetic exposure dates in a few minutes on desktop computers.Thus, these models do not have the large costs associated with collecting a representative number of samples from individual moraines.In these models, the user prescribes the age of the moraine and the nuclide production rate.Therefore, there are no confounding effects in the model experiments from errors in estimating these values.
In this paper, we use Monte Carlo models of two geomorphic processes that introduce biases into exposure dating to evaluate different methods of estimating moraine age from cosmogenic exposure dates.These processes are moraine degradation and inheritance, which we describe below.Our models are based on earlier work (e.g., Zreda et al., 1994;Hallet and Putkonen, 1994;Putkonen and Swanson, 2003;Benson et al., 2005; see also Muzikar, 2009).We expand on these groundbreaking studies in several ways.First, we provide explicit descriptions of the mathematical formulations of the models, pointing out the simplifying assumptions that are inherent in these formulations.We test the models' sensitivity to changes in their input parameters.Last, we provide code for these models that is written in MATLAB, an easily understood, high-level programming language.The model codes and documentation, as well as representative output from the models, are contained in the supplement (http://www.geosci-model-dev.net/3/293/2010/gmd-3-293-2010-supplement.zip).
Our results suggest that the skewnesses of individual data sets may provide a basis for deciding whether to take the mean, oldest date, or youngest date as the best estimate of moraine age.Preliminary tests indicate that this method is more robust in the presence of geomorphic effects than any of the commonly applied methods.
However, the skewness-based method for estimating moraine age sometimes yields misleading results.Thus, more sophisticated statistical methods are needed.Our models can be compared directly to data sets from individual moraines.This inverse modeling procedure (Applegate, 2009) yields explicit estimates of moraine age, as well as other parameters of geomorphic interest.

Numerical models
We describe models of two geomorphic processes that influence cosmogenic exposure dates from moraine boulders.These processes are moraine degradation and inheritance.In this section, we describe how our models treat these two processes, and we present preliminary results from these models.
These models are deliberately simplified.As an alternative approach, we might build a comprehensive model that would incorporate all of the processes that influence exposure dates on moraines.However, we wish to invert these models against observations, to allow direct estimation of moraine ages from collections of cosmogenic exposure dates (Applegate, 2009).In a model inversion, the maximum number of model parameters that can be estimated from a data set is typically smaller than the number of observations.Our models have three to five parameters each, and most collections of cosmogenic exposure dates from moraines contain about five observations (Putkonen and Swanson, 2003).Therefore, our models are already at the complexity limit imposed by the sizes of most available data sets.

The moraine degradation model
In moraine degradation, slope processes remove material from the crests of moraines and redeposit this material at the bases of the moraine slopes.The theoretical basis for understanding the redistribution of sediment on moraine slopes comes from observations made on fault scarps, wave-cut bluffs, and other landforms composed of unconsolidated sediment.These landforms become less steeply inclined and more rounded over time, suggesting that hillslope evolution can be modeled as a diffusive process (Nash, 1986;Hanks, 2000;Pelletier et al., 2006;Pelletier, 2008).That is, material moves downhill at a rate that is proportional to the local gradient.This observation implies that a sharp-crested moraine will lose height over its lifetime, as material moves from the moraine's crest to the toe of its slope (Anderson and Humphrey, 1989;Hallet and Putkonen, 1994;O'Neal, 2006;Putkonen et al., 2007;Pelletier, 2008).
Moraine degradation imparts a bias to cosmogenic exposure dates because it exposes boulders at the moraine crest that have been buried in sediment for some part of the moraine's history.Moraines typically contain large rocks distributed throughout a fine-grained matrix (Dreimanis, 1988;Benn and Evans, 1998).Because slope processes preferentially move fine-grained material, the boulders become concentrated on the crest of the moraine.Some of these boulders have been partly shielded from cosmic rays by the overlying sediment; they therefore contain smaller concentrations of cosmogenic nuclides than the boulders that have rested on the moraine crests since deposition of the moraine.The exhumed boulders yield cosmogenic exposure dates that underestimate the age of the moraine.
The model framework that we describe here builds on earlier studies.The use of slope evolution models to study moraines was first considered by Anderson and Humphrey (1989); Zreda et al. (1994) developed a model for the production of nuclides in boulders buried in an eroding surface.The first model of cosmogenic nuclide production on a diffusively evolving moraine was presented by Hallet and Putkonen (1994).This model was later developed further by Putkonen and Swanson (2003).Our model is closest to that of Putkonen and Swanson (2003).
To model the effects of slope processes on the height of moraines over time, we assume that moraines have an initial cross-section that is triangular, with an initial height h 0 and an initial slope S 0 , which is the (dimensionless) tangent of the slope in degrees.This profile evolves over time according to the one-dimensional diffusion equation, (1) (Hanks, 2000), where z(x, t) is the height of the moraine as a function of horizontal distance from the moraine crest x and time t; k is the topographic diffusivity (m 2 /yr).This rule assumes that k is constant over t and x (Pelletier et al., 2006;cf. Hallet and Putkonen, 1994;Roering et al., 2001).Solving this differential equation with our "sawtooth" initial moraine profile yields where Pelletier, 2008, his Eq. 2.45).Equation (2) agrees well with a Crank-Nicolson solution to Eq. ( 4) of Hallet and Putkonen (1994) if their β=0; compare Eq. ( 4) of Hallet and Putkonen (1994) to Eqs. (9.56) and (9.67) of Fletcher (1991).This analytical solution can be evaluated very quickly.As time goes on, the moraine's profile changes most at the crest and at the toe of the slope, becoming generally more rounded.As material is transported from the crest to the toe of the slope, the moraine becomes less tall.The moraine loses height rapidly at first, then more slowly.In (a), only one-half of the moraine's profile is shown; the modeled moraine is symmetrical about the y-axis.Note that the moraine loses more than 10 m of its initial height over 20 ka.Compare (a) to Fig. 1 of Hallet and Putkonen (1994); compare (b) to Fig. 1 of Putkonen and Swanson (2003).In this figure, the moraine's initial height is 50 m, its initial slope is 34 • , and its topographic diffusivity is 10 −2 m 2 /yr.
Setting x=0 in Eq. ( 2) yields an expression for the height of the moraine's crest as a function of time, (3) Figure 1 shows solutions to Eqs.
(2) and (3) for selected parameter values.The left panel (Fig. 1a) shows the moraine half-profile for elapsed time values of 5 ka, 10 ka, and 20 ka.
The moraine starts with a triangular profile, but becomes more rounded and less tall over time.The right panel (Fig. 1b) shows the height of the moraine as a function of time.The rate of crest lowering is rapid at first, then slows.
In both panels, the initial moraine height is 50 m, the initial moraine slope is 34 • (Putkonen and Swanson, 2003), and the topographic diffusivity is 10 −2 m 2 /yr (Hanks, 2000;Putkonen et al., 2007).These values seem reasonable for the large, last-glacial moraines of the western United States.Given Eq. (3), we can calculate the nuclide concentration in a boulder buried to some specified depth d 0 below the moraine's surface at the time of deposition.For purposes of calculating nuclide production rates, the depth of a boulder d(t) is given by Note that d 0 and d here refer to the depth of the top of the boulder, which is the point that will be sampled for cosmogenic nuclide measurements.in quartzite, following Granger and Muzikar (2001).The total production rate of beryllium-10 is roughly exponential as a function of depth; production is greatest at the surface, and falls off below the surface with an e-folding length of a few tens of centimeters (Lal, 1991).Most production near the surface is caused by high-energy protons and neutrons, which produce beryllium-10 by splitting atoms of oxygen and silicon in quartz (Gosse and Phillips, 2001).At greater depths, most production is due to muons, which do not interact with target atoms in the rock as easily as high-energy protons and neutrons.Compare this figure to Fig. 2a of Gosse and Phillips (2001).This figure assumes surface beryllium-10 production rates corresponding to sea level and high latitude and a rock density of 2.65 g/cm 3 .
Values of d 0 that exceed h 0 − h f are not meaningful, because these boulders will still be buried in the moraine at the time of sampling.By h f , we mean the final height of the moraine, achieved when t reaches the moraine's age.In addition, field geomorphologists typically do not sample boulders that stand less than some minimum height h b above the moraine crest (∼1 m; e.g., Gosse et al., 1995b).Thus, all the boulders that are sampled have values of d 0 that satisfy the criterion The production rate of most cosmogenic nuclides declines exponentially as a function of depth below material surfaces (Lal, 1991;see Zreda et al., 1994, for an important exception).That is, where P 0 is the production rate of the nuclide at the surface (atoms/g rock/yr), and is the attenuation length of cosmic rays in the material (∼160 g/cm 2 , divided by the material's density).We use the Lal/Stone production rates from the CRONUS online calculator (Balco et al., 2008) to estimate P 0 .
Equation ( 4) is a good approximation only at shallow depths, where nucleon production dominates; at greater values of d(t), muon production becomes important (Gosse and Phillips, 2001).To account for muon production, we use the parameterization of Granger and Muzikar (2001, their Eqs. 1-3).This scheme represents production at a given depth as the sum of four exponential terms, each with its own P 0 and .That is, We scale these terms relative to their values at sea level and high latitude, again using the CRONUS online calculator (Balco et al., 2008).This expression is a parameterization; Heisinger et al. (2002a, b) present alternative expressions that resolve the underlying physics.We use the relationship presented in Eq. ( 5) because it can be evaluated very quickly as a vector calculation in MATLAB.The speed of evaluation is important because this calculation must be performed millions of times for each model run.
Figure 2a shows the production rate of the cosmogenic nuclide beryllium-10 as a function of depth.Nucleon production dwarfs muon production at the surface, but muon production becomes increasingly important at greater depths (Fig. 2b).
Given Eqs. ( 3) and ( 5), we can calculate the final concentration of cosmogenic nuclides in a moraine boulder.For an unstable nuclide accumulating in a rock surface, the change in concentration C over an infinitesimal time is given by (Lal, 1991).λ is the decay constant of the appropriate nuclide (yr −1 ; Gosse and Phillips, 2001;Balco et al., 2008).After the surface has been exposed to cosmic rays for a time t f , the final nuclide concentration C f is given by That is, the final concentration is the production integrated over time, less the amount lost to nuclear decay.For simple forms of P (d(t)), analytical expressions for C f can be written (e.g.Lal, 1991, his Eq. 6).In contrast to this simple case, we model production in a large number of boulders that have different depth-through-time trajectories.The model production rate in a given boulder is a piecewise function of time, because the production rate stops changing when the boulder breaks the surface of the moraine (that is, when d becomes 0; Fig. 3).Therefore, we break the lifetime of the moraine into n time steps, each having a duration t.We then evaluate the change in concentration in all the modeled boulders during each of these time steps.For the model runs shown in this paper, we used values of t ranging from 25 yr to 100 yr.Note that the initial concentration C 0 is taken to be zero here; we treat inheritance in the next section.If boulders are uniformly distributed throughout the till, then some boulders will be at the surface when the moraine is deposited, whereas other boulders will be present in the till at greater depths.As time goes forward, the moraine loses height (Fig. 1), and the boulders approach the surface.At the same time, cosmogenic nuclides are produced in the boulders (Fig. 2).For buried boulders, production rates increase slowly as the surface lowers, then become constant after the boulders are exposed at the surface.Note that the majority of the cosmogenic nuclides in each boulder are produced after the boulder reaches the surface, even for the most deeply buried boulder.In (b), the dots indicate the time when each boulder reaches the surface.As in Fig. 1, the moraine's initial height is 50 m, its initial slope is 34 • , and its topographic diffusivity is 10 −2 m 2 /yr.The final heights of all the boulders are 1 m.
Figure 3a shows the depths of four boulders within the moraine as a function of time, assuming the same model parameters as in Fig. 1.At the beginning of the simulation, one boulder is at the surface (d 0 =0 m), another boulder is buried to a depth of 9 m (d 0 =9 m), and the other two boulders are evenly spaced between these depths.As the moraine loses height over time, the boulders approach the surface and are eventually exposed at the surface.Compare this figure to Fig. 1b.
Figure 3b shows the concentrations of beryllium-10 as a function of time in each of the boulders whose depth www.geosci-model-dev.net/3/293/2010/Geosci.Model Dev., 3, 293-307, 2010 trajectories are shown in Fig. 3a.Again, the model parameters used to generate this figure are the same as those in Fig. 1.The concentrations in the boulders increase slowly while the boulders are still buried in the moraine; after they reach the surface, the concentration increases roughly in proportion to surface residence time, less a small amount for nuclear decay (Eq.6).Note that the bulk of the final nuclide concentration in each boulder is acquired only after the boulder reaches the surface, even for the boulder that is buried most deeply in the moraine at the beginning of the simulation.This figure assumes that the beryllium-10 concentrations in all the boulders are zero when the simulation begins.
Although we do not consider boulder erosion in this paper, the model treats erosion by the progressive removal of thin shells of material from boulder surfaces after they are exhumed from the till.In contrast to Hallet and Putkonen (1994), we do not allow boulders to shrink below the observed boulder height h b (see Zreda et al., 1994).Instead, we determine the amount of time that each boulder will be exposed to surface weathering from Eq. ( 3), then specify initial sizes for the boulders that will result in the boulders having the observed height.
This model assumes that exhumed boulders do not topple or rotate as the crest of the moraine deflates.It also neglects the effects of cryoturbation (Lal and Chen, 2005).Toppling or rotation of boulders on a degrading moraine would produce a larger range of exposure dates than degradation alone, because these processes effectively reduce the measured nuclide concentrations in sampled boulders (Ivy-Ochs et al., 2007;Schaefer et al., 2008).Conversely, cryoturbation might bring boulders to the moraine surface sooner than would be predicted by diffusive removal of the moraine crest, thereby reducing the range of exposure dates from the moraine.In this paper, we assume that these processes are of secondary importance compared to the ones we do treat.Some moraines have geomorphic characteristics that are inconsistent with the assumptions used in constructing the moraine degradation model.For example, it would be inappropriate to apply our model of moraine degradation to the large Pinedale terminal moraines near Pinedale, Wyoming (Richmond, 1973;Gosse et al., 1995a), particularly in the Halls Lake (Mud Lake) drainage.These moraines have broad, flat crests, where the local slope is close to zero.Consequently, the downhill flux of material at the crests of these moraines should be small.We expect that these moraines have lost little material from their crests over time.Moreover, limited exposures in roadcuts at Fremont Lake show that there are few or no boulders in the subsurface till (E.Evenson, personal communication, 2008).These observations invalidate the assumption that the boulders are uniformly distributed throughout the outermost Pinedale-age moraine at Fremont Lake.

The inheritance model
Boulders that are deposited on a moraine with nonzero concentrations of cosmogenic nuclides are said to have inheritance.The inherited nuclides were produced in each boulder during one or several periods of "pre-exposure" (Ivy-Ochs et al., 2007).That is, the boulders were incompletely shielded from cosmic rays before being deposited on the moraine.These boulders contain larger concentrations of cosmogenic nuclides than boulders that were completely shielded from cosmic rays at all times before being incorporated into the moraine.Exposure dates from boulders with inherited nuclides tend to overestimate the age of the moraine.
There are at least two potential sources of pre-exposed boulders in glaciated landscapes (Ivy-Ochs et al., 2007).First, boulders may topple onto the glacier surface from cirque headwalls or adjacent, oversteepened valley walls (Seong et al., 2009).These boulders then ride the glacier's surface to the terminus, where they fall onto the moraine.Second, glaciers may re-entrain boulders deposited in the valley bottom during an earlier advance, or pluck boulders from bedrock outcrops at the glacier bed.These boulders are then transported subglacially to the glacier terminus, where they may be emplaced at the moraine surface by thrusting (e.g., Krüger, 1996) or other ice-marginal processes.
The mathematical descriptions of these two situations are nearly identical.In both cases, the concentration measured in each boulder is the sum of the inherited component acquired during pre-exposure, and the post-depositional component that reflects the exposure history of the boulder after moraine construction.
The model that we describe here is based on an earlier model presented by Benson et al. (2005), which treated inheritance in boulders derived from cirque headwalls.Our model uses a mathematical formulation that is similar to the one used by Benson et al. (2005), but treats a larger set of geomorphic situations.In addition, our model of inheritance is similar to the model of nuclide concentrations in sediment over time used in cosmogenic burial dating (Granger et al., 2001;Granger and Muzikar, 2001).Following this pioneering work, we assume that the sampled clasts had two distinct periods of residence in the landscape, and that the rate of change of nuclide concentrations in the clasts was different during these two periods.
For simplicity, we begin by describing the model treatment of inheritance in reworked boulders.We then point out a slight change in the model formulation that allows it to treat inheritance in boulders derived from cirque headwalls and valley walls.
For reworked boulders, the inherited concentration in each boulder depends on the time between deposition of the boulder by the retreating ice and entrainment of the boulder by the readvancing glacier t pre , and on how deeply the boulder was buried during this time d pre .Both these parameters are unknown for any individual boulder, but it is reasonable to say that they must range from zero to some maximum.0 ≤ t pre ≤ max(t pre ), and 0 ≤ d pre ≤ max(d pre ).
The maximum time max(t pre ) represents the time between the beginning of the penultimate glacial retreat and the time of moraine deposition; the maximum depth max(d pre ) is the maximum thickness of material eroded by the glacier during its readvance.
Note that d pre refers to the depth of the point on each boulder that is eventually sampled, not the top of the boulder, during the predepositional exposure time.Field geomorphologists typically sample the upper surfaces of boulders, because those surfaces receive the maximum flux of cosmic rays.However, glacial transport rotates boulders, and so the sample point is not necessarily the same as the apex of the boulder during the predepositional exposure time.Sampling of the sides of moraine boulders yields a range of nuclide concentrations (Schaefer et al., 2008), consistent with theoretical predictions of the distribution of nuclide production in solids (Masarik and Wieler, 2003;Lal and Chen, 2005).
For a boulder buried in a till sheet, Eq. ( 5) gives the production rate in the point that is eventually sampled.Given this production the inherited concentration C pre is (Lal, 1991, his Eq. 6), and the final concentration C f , achieved after the boulder has rested on the moraine for a time t, is (Lal, 1991, his Eq. 6; note that the depth is 0 for this application).Here, ε is the erosion rate of the boulders after they are delivered to the moraine (cm/yr; assumed negligible), and is the attenuation length of the nucleonic component of cosmogenic nuclide production (∼160 g/cm 2 , divided by the material's density; Lal, 1991;Gosse and Phillips, 2001).
Our model is readily adapted to treat inheritance in boulders derived from cirque headwalls and valley walls, as in Benson et al. (2005).From a nuclide production perspective, the angle of the overlying surface is the critical difference between a boulder buried in a till sheet and one that is still in a cirque headwall; for a till sheet, the overlying surface should be nearly horizontal, whereas cirque headwalls are quite steep.To model nuclide production as a function of depth below inclined surfaces, we use the parameterization of Dunne et al. (1999, their Eq. 18).This parameterization gives results within 3% of estimates from a more explicit model (Dunne et al., 1999), even for the steep slopes representative of cirque headwalls (∼30 • ; Benson et al., 2005).
The model implicitly accounts for the rotation of boulders during glacial transport.Because glacial transport mixes sediment and boulders, most previously exposed boulders will arrive on the moraine in a different orientation than they had during their predepositional exposure times.Thus, for a cube-shaped boulder, there is a 1-in-6 chance that the face that is eventually sampled is the one with the largest concentration of inherited cosmogenic nuclides (Benson et al., 2005).Because our model is formulated in terms of the depth of the sampled point on each boulder below the predepositional exposure surface, the inherited nuclide concentrations are insensitive to the boulders' orientations during the predepositional exposure time.This statement will be true as long as the density contrast between the boulders and the surrounding material is small.If there is a large difference in density between the boulders and the surrounding material during the predepositional exposure period, the production rate in the sampled points will differ, depending on the orientations of the boulders.
This inheritance model relies on many assumptions.First, we assume that there are no nuclides inherited from any periods of residence in the landscape preceding the last glacial cycle.Because many cosmogenic nuclides have halflives that are long compared to glacial cycles (Gosse and Phillips, 2001;Shackleton, 2000), this assumption requires that glaciers sweep out most of the easily eroded material from their valleys during each advance.Second, we assume that surface production rates were the same during the predepositional exposure time as they are in the boulders' observed positions.Because some boulders are undoubtedly coming from higher elevations than the present-day moraine crests, this assumption tends to underestimate surface production rates during the predepositional exposure time.Future versions of this model will need to incorporate information on the elevation distribution of glaciated basins (e.g., Bierman et al., 2005).For boulders that travel to the moraine atop glacial ice, some cosmogenic nuclide atoms are produced during the transport time (Seong et al., 2009), and our model neglects this production.Moreover, glaciers do erode boulders during subglacial transport, and this model does not include that process.We tolerate these problems for the sake of developing this preliminary model.

Monte Carlo simulation
To determine a statistical distribution of apparent exposure dates from our models, we use Monte Carlo methods (Hilborn and Mangel, 1997;Bevington and Robinson, 2003).In Monte Carlo simulation, the values of model parameters are chosen randomly from predefined probability distributions.The model is then run for these parameter values, and the output is saved.This process is repeated many times; depending on the speed of the model and the desired precision, www.geosci-model-dev.net/3/293/2010/Geosci.Model Dev., 3, 293-307, 2010 Monte Carlo model evaluations may include thousands to millions of individual model runs.The model output is then plotted as a histogram.In our models, there are several free parameters that will be different for each boulder on a moraine.We have no way of determining, for example, how deeply buried any individual boulder was at the time of moraine deposition.The moraine degradation model has only one highly variable parameter, the initial depth d 0 ; the inheritance model has two highly variable parameters, the predepositional exposure time t pre and the depth during the predepositional exposure time d pre .
Because all these free parameters range from zero to some maximum, we choose random values for these parameters from continuous uniform distributions.In a continuous uniform distribution, all real numbers that lie between the minimum and maximum ends of the distribution are equally probable (Hilborn and Mangel, 1997;Bevington and Robinson, 2003).For our models, the minimum ends of these distributions are always 0; the maximum ends are specified by max(d 0 ), max(t pre ), and max(d pre ).
For each draw of these randomly chosen parameter values, we calculate the final concentration C f (Eqs.7 and 9, above) and the apparent exposure time t app , according to (Lal, 1991, his Eq. 6).This expression reflects the "naïve" estimate (Wolkowinsky and Granger, 2004) of moraine age from a single boulder sample, ignoring boulder erosion and all other geomorphic processes.
Note that we differentiate between moraine-level parameters and boulder-level parameters.Moraine-level parameters in the degradation model include the moraine age, topographic diffusivity, initial height, and initial slope; in the inheritance model, the moraine-level parameters are the moraine age, the maximum predepositional exposure time, and the maximum predepositional burial depth.The boulderlevel parameters are the initial depth of boulders below the moraine surface in the degradation model, and the predepositional exposure time and burial depth in the inheritance model.In estimating the probability distribution of cosmogenic exposure dates from a single moraine, we vary the boulder-level parameters, but the moraine-level parameters remain constant.

Plotting non-normal distributions
Many common methods of plotting collections of exposure dates from moraines implicitly assume that the dates are drawn from a normal distribution.However, the statistical distributions of exposure dates produced by our models are clearly not normal.Therefore, we represent the statistical distributions of modeled exposure dates using histograms, cumulative density functions, and box plots (Chambers et al., 1983;Croarkin and Tobias, 2006).These plotting methods are robust, even for statistical distributions that vary considerably from the normal distribution.
Histograms are probably the most familiar method of representing distributed data, but the choice of bin size exerts a strong control on the shape of the histogram.In a histogram, the synthetic observations are sorted into bins.The heights of the bars on the histogram are proportional to the number of observations in each bin.
Unlike histograms, plots of cumulative density functions do not require arbitrary choices about how to group the data.On a plot of a cumulative density function, the y-axis represents the probability that any individual observation is equal to or less than a particular value on the x-axis (Press et al., 1992, their chapter 14;Hilborn and Mangel, 1997;Croarkin and Tobias, 2006).The x-axis therefore ranges from the minimum to the maximum of the observations; the y-axis ranges from 0 to 1.0.
Box plots provide a compact way of representing distributed data; placing several box plots next to one another allows quick comparison of distributions.In a box plot, the position and width of the box indicates where the middle 50% of the observations lie.That is, the box represents the interquartile range of the data (Chambers et al., 1983;Croarkin and Tobias, 2006).The line in the box is the median, or the value that separates the lower half of the observations from the upper half.In this paper, the ends of the whiskers indicate the positions of the largest and smallest observations.Often, box plots indicate outliers as dots or small crosses outside the whiskers (Chambers et al., 1983), but we do not follow this practice.

Model output for representative parameter values
The output from the moraine degradation model is shown in Fig. 4. Figure 4 assumes the same parameter values used in Figs. 1 and 3; the initial height of the moraine is 50 m, the initial slope of the moraine is 34 • (Putkonen and Swanson, 2003), the topographic diffusivity is 10 −2 m 2 /yr (Hanks et al., 2000;Putkonen et al., 2007), and the age of the moraine is 20 ka.In addition, we specify that the tops of all sampled boulders must be at least 1 m above the crest of the moraine at the time of sampling.
Figure 4a illustrates the relationship between the initial depth of a given boulder and the apparent exposure time yielded by that boulder.As expected, the more deeply buried samples yield younger apparent exposure times.
Figure 4b and c shows the statistical distribution of the exposure dates produced by the degradation model for these parameter values.The distribution is strongly skewed toward old values; that is, more of the probability mass falls to the young side of the distribution's peak than would be the case if the distribution were normal.The corresponding  (Chambers et al., 1983).As in Figures 2 and 4, the moraine's initial height is 50 m, its initial slope is 34°, and its topographic diffusivity is 10 -2 m 2 /yr.The final heights of all the boulders are 1 m.(Chambers et al., 1983).As in Figs. 1 and 3, the moraine's initial height is 50 m, its initial slope is 34 • , and its topographic diffusivity is 10 −2 m 2 /yr.The final heights of all the boulders are 1 m.cumulative density function rises slowly, then more rapidly as it approaches the true age of the moraine (20 ka).The box portion of the box plot, which represents the position of the bulk of the data, falls on the right-hand side of the plot.
The output from the inheritance model is shown in Fig. 5ac.These plots assume a moraine age of 20 ka, a maximum predepositional exposure time of 100 ka, a maximum depth during the predepositional exposure period of 2 m, an overburden density of 2.0 g/cm 3 , and a flat surface geometry during the predepositional exposure period.Again, the total number of synthetic observations in each of these plots is 10 5 .
Figure 5a shows contours of the apparent exposure time produced by the inheritance model as a function of the model's free parameters, predepositional exposure time and predepositional exposure depth.As expected, the samples that yield the greatest apparent exposure times are those that had the greatest length of time to acquire inherited nuclides and were near the surface during that time.That is, the samples that appear oldest have the longest predepositional exposure times and the smallest predepositional exposure depths.
Figure 5b and c shows the statistical distributions of exposure dates expected from the inheritance model for these parameter values.The distribution is skewed toward old values; it contains a mode close to the true age of the moraine (20 ka), and a long, heavy tail to the old side, as shown in the histogram (Fig. 5b).These features of the distribution are reflected in the cumulative density function (Fig. 5c), which rises rapidly, then levels off.The box portion of the box plot falls near the left end of the plot.

Sensitivity of modeled distributions to input parameter values
Some of the parameters used in our models are either highly uncertain, or else vary considerably between moraines.In this section, we show how the modeled distributions of exposure dates change as individual parameters vary.Figures 6 and 7 illustrate the sensitivity of the two models using box plots (Chambers et al., 1983).
In both models, the moraine age controls the position of the box plot along the time axis.In the inheritance model, the spread of the exposure dates is independent of moraine age; the distance between the ends of the whiskers is the same for all values of moraine age.In contrast, the moraine age does affect the spread of exposure dates yielded by the degradation model; that is, younger moraines show less spread than older moraines (Fig. 6; Putkonen and Swanson, 2003).The increase in spread among exposure dates with age for degrading moraines happens because older moraines have more time to lose material from their crests (Fig. 1), and this process exposes more boulders that have spent progressively less time exposed to the full surface flux of cosmic rays (Fig. 3).
In the degradation model, the spread of dates is most strongly controlled by the topographic diffusivity, although the initial slope and initial height of the moraine also have some influence on the scatter (Fig. 6).Small diffusivities cause the moraine's height to change only slightly over its lifetime, and so few new boulders are exhumed at the crest of the moraine.Very large diffusivities flatten the moraine in a few thousand years after its construction; the reduced spread in exposure dates produced by the model for a diffusivity of 1 m 2 /yr happens because such a high diffusivity exposes most of the buried boulders within a few thousand years after the deposition of the moraine.Such large diffusivities cause the moraine to disappear almost totally over 20 ka, so they are inconsistent with the observed persistence in the landscape of topographically distinct moraines (see Hanks, 2000;Putkonen et al., 2007).The modeled distributions of exposure dates from tall moraines are wider than distributions from less tall moraines of the same age (Putkonen and Swanson, 2003), although the width of the distribution stops increasing as the initial height of the moraine is made greater than ∼35 m.The range of modeled exposure dates increases monotonically with the initial slope of the moraine.
In the inheritance model, the maximum predepositional exposure time controls the width of the distribution, and the maximum predepositional exposure depth controls where the bulk of the dates falls between the extreme ends of the As in Fig. 5, the base values for the input parameters specify that the true age of the moraine is 20 ka, the maximum predepositional exposure time is 100 ka, and the maximum predepositional burial depth is 2.0 m. distribution (Fig. 7).A large value for the maximum predepositional exposure time causes a wide range of exposure dates; a small value produces a narrow range.Large values of the maximum predepositional exposure time concentrate most of the observations near the young end of the range, whereas smaller values place more of the observations into the tail of the distribution.
Increasing the surface slope has only a small effect on the distributions of exposure dates produced by the inheritance model (Fig. 7).There is little difference between the distributions of modeled exposure dates for boulders derived from flat surfaces and those for boulders derived from sloped surfaces with inclinations of 30 • or less, because the depth dependence of nuclide production changes only slightly over this range of slopes (Dunne, 1999).A 30 • slope is representative of cirque headwalls (Benson et al., 2005), which are a likely source for supraglacial boulders.The model sensitivity to surface slope is not extreme, even for larger slope values.

Applications
The sensitivity analysis presented in Sect. 3 shows that there are clear differences between the statistical distributions pro- Even large data sets (n=21) can provide a misleading estimate of the skewness of the parent distribution.In particular, randomly chosen data sets will often yield skewnesses that do not have the same sign as the underlying parent distribution.The parent distribution in the top panel is the same as that shown in Fig. 5b; the parent distribution in the bottom panel is shown in Fig. 4b.The parent distribution in the middle panel is a normal distribution with a mean of 20 ka and a standard deviation of 1 ka.
duced by the moraine degradation and inheritance models.Moreover, these differences are robust for the parameter combinations we have tested.The statistical distributions of exposure dates produced by the moraine degradation model are always skewed toward young values (Fig. 4b); conversely, the distributions of exposure dates produced by the inheritance model are always skewed toward old values (Fig. 5b).That is, the cumulative density functions from the degradation model are always concave-up (Fig. 4c), and the cumulative density functions from the inheritance model are always concave-down (Fig. 5c).On the box plots, the box occurs near the right-hand end of the distribution in the degradation model (Fig. 6), and near the left-hand side of the plot for the inheritance model (Fig. 7).This result is obvious in hindsight; however, this contrast between the effects of moraine degradation and inheritance on exposure dating of moraine boulders has never been shown before.So far as we know, histograms of modeled cosmogenic exposure dates from a degrading moraine have appeared in the literature only once (Zreda et al., 1994), and there is no corresponding figure for inheritance.
Fig. 9.The reliability of different interpretive methods in estimating moraine ages from collections of cosmogenic exposure dates.Each box plot (Sect.2.3) represents age estimates from 10 6 randomly selected data sets containing eight synthetic cosmogenic exposure dates each.The heavy, vertical black line in each panel represents the true age of the moraine, which is 20 ka in each case.In each panel, the methods listed on the y-axis are listed according to how close the median age estimate falls to the true moraine age; the method listed at the top is the best for the indicated parent distribution, and the method listed at the bottom is the worst.This ordering is insensitive to the number of samples in each data set for reasonable data set sizes (3 ≤ n ≤ 21).The parent distribution in the middle panel is a normal distribution with a mean of 20 ka and a standard deviation of 1 ka, corresponding to a case where all of the scatter between the exposure dates is due to measurement error.The parent distributions in the other two panels are those shown in Figs. 4 and 5.

Skewness as a guide to geomorphic process
Our results suggest a potential method for distinguishing moraine degradation from inheritance.If the skewness of a data set is strongly positive, it suggests that the boulders contain inherited nuclides.Similarly, if the skewness is strongly negative, the dates may be biased by moraine degradation.We test this hypothesis by calculating the skewnesses of synthetic data sets drawn at random from the distributions produced by our models (Figs. 4 and 5), as well as a representative normal distribution.This normal distribution represents the case where there is no geomorphic bias, and all the scatter among exposure dates is due to measurement error.It has a mean of 20 ka and a standard deviation of 1 ka, corresponding to a moraine deposited during the Last Glacial Maximum and a measurement error of 5%.
The results of this numerical experiment are shown in Fig. 8.In general, the skewnesses of reasonably-sized data sets (n <∼ 20) are a poor guide to the skewness of the underlying parent distribution; in many cases, the skewness of a particular data set will have the opposite sign from the parent distribution.However, this figure does place rough bounds on the range of skewnesses that the measurement error-only case can produce.Skewnesses that lie outside the boxes in the middle panel of Fig. 8 should be taken as evidence for either moraine degradation or inheritance, depending on whether the skewness is negative or positive.

The effectiveness of different methods for estimating the ages of moraines from cosmogenic exposure dates
Next, we evaluate the robustness of different methods of estimating moraine age in the presence of geomorphic bias.
As noted in the introduction, many methods for estimating moraine age from cosmogenic exposure dates appear in the literature.These methods include the mean, the mean after excluding outliers, the oldest date, and the youngest date.In this case, we define outliers as those observations that are more than twice the standard deviation away from the mean of the exposure dates in a data set.Other methods also exist; the mode and the weighted mean (Bevington and Robinson, 2003) are both common choices.
To this list, we add the min/mean/max technique, which was suggested by our modeling results.Given a collection of exposure dates from the same moraine, we estimate the moraine's age using the min/mean/max technique as follows.If the skewness is greater than 0.5, we infer that the dates are biased by inheritance, so we take the youngest date.If the skewness is less than −0.5, we assume moraine degradation, and take the oldest date.If the skewness is between −0.5 and 0.5, we take the mean.The cutoff of 0.5 corresponds to the ends of the boxes in the middle panel of Fig. 8.Other cutoff values could also be tried.
We evaluate these methods by using them to estimate moraine ages for synthetic data sets drawn from our modeled distributions, for which the true ages are known.For the purposes of this study, we restrict ourselves to data sets containing eight exposure dates.However, we find that our results are robust for reasonable data set sizes (n < 21).
Figure 9 shows the results of this experiment.As we might expect, the mean performs best in the measurement erroronly case, and the extreme estimators (oldest date, youngest date) do very well when applied to data sets produced by the appropriate models.However, these methods fail when applied inappropriately.
Discarding outliers before taking the mean of exposure dates is a common practice, yet Fig. 9 suggests that this method performs no better than simply taking the mean in cases where there is geomorphic bias.This result may be understood from Figs. 4b and 5b.Some dates fall far from the moraine's true age, yet the bulk of the exposure dates have much smaller degrees of bias.Thus, discarding the extreme outliers and taking the average of the rest still yields an answer that is too young or too old, depending on the underlying geomorphic process (Applegate et al., 2008).It may be useful to think of geomorphic bias as a dial that is set to a different value for each boulder on a moraine, rather than a switch that is either on or off.
Thus, we may prefer an estimator that is not sensitive to the underlying distribution.The min/mean/max technique appears to satisfy this requirement; the median age estimate using this method is within a few thousand years of the true moraine age for each of the parent distributions we tested (Fig. 9).This statement is not true of any of the other methods we examined.However, the min/mean/max technique sometimes yields age estimates that are wrong by thousands of years.

Discussion
In this paper, we described numerical models of two geomorphic processes, moraine degradation and inheritance, and we used the results of these models to examine the effectiveness of different procedures for estimating the ages of moraines.Given a large number of samples, the statistical distributions of exposure dates from degrading moraines should be skewed toward young values, whereas the statistical distributions of exposure dates from boulders containing inherited nuclides should be skewed toward old values.Thus, the skewness may provide a means for determining whether the scatter among exposure dates from a particular moraine is due to moraine degradation, inheritance, or measurement error only.The min/mean/max technique for estimating moraine age uses a cutoff value for the skewness to decide which simple estimator to use.This method appears to be more robust to geomorphic biases than any of the other estimators we tested.
Despite the potential usefulness of the min/mean/max technique for estimating moraine age, more sophisticated statistical methods are probably needed.In using the skewness to decide which simple estimator to use, we sometimes guess incorrectly (Figs. 8 and 9).For example, we might infer from a positive skewness that a data set is biased by inheritance, when the scatter among the dates is actually due to moraine degradation.Thus, we pick the oldest date, whereas the youngest date would be more appropriate.In such cases, the min/mean/max technique will err by thousands of years.
Therefore, we have developed a method for inverting our process models against data sets from individual moraines (Applegate, 2009).This approach uses the Kolmogorov-Smirnov test statistic to find the best fit between observed and modeled distributions of cosmogenic exposure dates.From this exercise, we learn the most likely values for the model parameters, including moraine age, topographic diffusivity, and the depth of glacial erosion.Thus, this inverse method provide a means to learn about geomorphic processes at individual field sites from cosmogenic exposure dates.
The inverse method does require a large number of dates from each moraine to be useful, so we plan to develop the min/mean/max technique further.In particular, the choice of 0.5 for the skewness cutoff is somewhat arbitrary.It should be possible to determine the skewness cutoff on a case-bycase basis, in a way that minimizes the chance of incorrect process assessments.
We developed our models of moraine degradation and inheritance separately, to emphasize the contrasts between the modeled distributions; however, moraine degradation and inheritance undoubtedly co-occur on some moraines (Laabs et al., 2009).In addition, we have neglected some important processes, such as boulder erosion.Thus, our statements about the accuracy of different methods for estimating moraine ages should be reevaluated as more complete models are developed.

30Figure 4 :
Figure 4: Distribution of cosmogenic exposure dates produced by the moraine degradation model for a representative case.Panel (a) shows the exposure dates yielded by boulders as a function of their initial burial depth in the moraine (compare Fig. 3b).Panel (b) shows a histogram of these apparent ages.Most of the exposure dates cluster around the true age of the moraine (20 ka), but there is a long, heavy tail to the left.That is, the distribution of exposure dates produced by the moraine degradation model is skewed toward young values.The total number of observations shown in this histogram is 10 5 .Panels (c) and (d) show the cumulative density function and box plot of the 10 5 observations shown in the histogram.Panel (d) is immediately below panel (c).Dashed lines in (c) and (d) show the relationship of the box plot to the cumulative density function; breaks in the box plot represent the quartiles of the distribution(Chambers et al., 1983).As in Figures2 and 4, the moraine's initial height is 50 m, its initial slope is 34°, and its topographic diffusivity is 10 -2 m 2 /yr.The final heights of all the boulders are 1 m.

Fig. 4 .
Fig. 4. Distribution of cosmogenic exposure dates produced by the moraine degradation model for a representative case.Panel (a) shows the exposure dates yielded by boulders as a function of their initial burial depth in the moraine (compare Fig. 3b).Panel (b) shows a histogram of these apparent ages.Most of the exposure dates cluster around the true age of the moraine (20 ka), but there is a long, heavy tail to the left.That is, the distribution of exposure dates produced by the moraine degradation model is skewed toward young values.The total number of observations shown in this histogram is 10 5 .Panels (c) and (d) show the cumulative density function and box plot of the 10 5 observations shown in the histogram.Panel (d) is immediately below panel (c).Dashed lines in (c) and (d) show the relationship of the box plot to the cumulative density function; breaks in the box plot represent the quartiles of the distribution(Chambers et al., 1983).As in Figs.1 and 3, the moraine's initial height is 50 m, its initial slope is 34 • , and its topographic diffusivity is 10 −2 m 2 /yr.The final heights of all the boulders are 1 m.

Fig. 5 .Figure 6 Fig. 6 .
Fig.5.Distribution of cosmogenic exposure dates produced by the inheritance model for a representative case.Panel (a) shows contours of the apparent ages yielded by boulders as a function of the length of time that they were exposed to cosmic rays and the depth to which they were buried during that time.Panel (b) shows a histogram of exposure dates produced by random sampling of 10 5 synthetic observations from the contour plot in (a).In contrast to the distribution produced by the moraine degradation model (Fig.4), the inheritance model produces distributions that are skewed toward old values.The bulk of the exposure dates fall near the true age of the moraine (20 ka), but there is a long, heavy tail to the right.Panels (c) and (d) show the cumulative density function and box plot of the 10 5 observations shown in the histogram.Panel (d) is immediately below panel (c).The true age of the moraine is 20 ka, the maximum predepositional exposure time is 100 ka, and the maximum predepositional burial depth is 2.0 m.