Earth and Planetary Science Letters Paired-cosmogenic nuclide paleoaltimetry

Goal – The reconstruction of past topographies remains challenging and only a few methods allow accurate determination of past surface elevations. We propose here a new technique for deriving paleoelevations, in which multiple cosmogenic nuclides are measured in the same geological sample exposed at the Earth’s surface. This method relies on the altitude dependence of the cosmogenic nuclides’ production rates combined with the radioactive decays of nuclides with different half-lives. Theory – The position of the two cosmogenic nuclide exposure curves ( 26 Al/ 10 Be vs 10 Be or 10 Be/ 21 Ne vs 10 Be) depends on the altitude of exposure. If the studied surfaces have been exposed for suﬃciently long durations ( > 500 ka), or have been affected by low erosion rates ( < 1 mMa − 1 ), measurement of two cosmogenic nuclides with different half-lives thus allow accurate elevations to be determined with a reasonable uncertainty ( < 1000 m at 1 σ ).


Introduction
Cosmogenic nuclide techniques are powerful tools for addressing a wide variety of scientific problems in Earth sciences (e.g. Dunai, 2010;Granger et al., 2013). The majority of these studies have been based on the surface exposure dating of geomorphic surfaces (e.g. Gosse et al., 1995)  averaged erosion rates (e.g. Granger et al., 1996). For these quite straightforward applications, it is generally sufficient to measure only one nuclide in the suitable rock samples (e.g. Gosse and Phillips, 2001;Lal, 1991). Measuring two cosmogenic nuclides with different half-lives in the same rock allows more complex exposure histories to be addressed, and notably, the determination of burial ages, i.e. the time since a geological sample has been shielded from cosmic rays (e.g. Granger and Muzikar, 2001). This application has led to many breakthrough publications since it is one of the most accurate and efficient methods available for dating Hominid remnants in continental deposits (e.g. Granger et al., 2015;Lebatard et al., 2014).
Nevertheless, it is crucial to keep exploring new applications, in particular to tackle unsolved scientific questions in Earth sciences. Developing a method for determining accurate paleoelevations over geological timescales remains a challenge. Reconstructing past elevations of a mountain range may provide important clues regarding the geodynamic processes involved (e.g. Husson and Sempere, 2003;Molnar et al., 1993). Relief changes are also thought to have a major impact on atmospheric circulation and global paleoclimates (e.g. Molnar and England, 1990). Existing paleoaltimetry methods frequently suffer from several limitations in that many of them often requiring loosely constrained assumptions (Clark et al., 2007). This has led to contrasting results in some regions and leaves open the debate on the amplitude of the impact of tectonics on global climatic changes (e.g. Boos and Kuang, 2010;Licht et al., 2014).
Because the Earth's atmosphere significantly attenuates the cosmic-ray flux, the production rate of a cosmogenic nuclide is sensitive to elevation (Lal, 1991). This property has been explored as a means for retrieving paleoelevation, either from continuously exposed surfaces (Brook et al., 1995;Evenstar et al., 2015) or from buried paleo-surfaces that have been exposed for a known duration in the past (Blard et al., 2005). However, these exploratory works had to address several important limitations, notably because a method based on the measurement of only one nuclide corresponds to a mathematically underdetermined system: this approach thus requires the use of surfaces that have not been affected by erosion, and for which the duration of paleo-exposure can be independently and accurately measured with high precision (Blard et al., 2006). Such a combination of favorable conditions is rarely fulfilled in nature, making the method difficult to apply.
Here we present a novel approach for using and interpreting datasets of multiple-cosmogenic nuclides measured in the same rock samples. We show that, when the exposure time is long enough (typically > 500 ka), the isotopic ratio of two cosmogenic nuclides with contrasting half-lives, along with their respective concentrations, allows calculation of precise and accurate elevations. Under the most favorable conditions, the method enables altitudes to be derived with a precision (<500 m) that is useful for most geological applications. The strength and main advantage of this new technique is that its accuracy is only slightly sensitive (bias from this assumption is less than 1000 m) to the amount of erosion that has affected the exposed paleosurface, a parameter that is generally unknown. We tested the method on an existing dataset of samples collected in the arid part of the Western Andes, using bedrock and detrital samples in which pairs of cosmogenic nuclides ( 10 Be-21 Ne and 26 Al-10 Be) had been measured in previous studies (Kober et al., 2007;Nishiizumi et al., 2005;Placzek et al., 2010;Ritter et al., 2018a). The calculated altitudes are in good agreement with the present-day elevation of the samples, suggesting both that the method is adequate and that the altitude change of the Central Andes has been limited over the last hundreds of thousands of years. We also performed Monte Carlo simulations to explore the detection limits, the uncertainties and the range of applicability of the method over the geological timescale, for both the 26 Al-10 Be and the 10 Be-21 Ne pairs. This modeling shows that the precision and accuracy of this paleoaltimetry method both increase significantly with exposure time or with the low erosion rates of the pre-burial paleo-surface. Surfaces that have been exposed for longer than 500 ka years (or affected by erosion rates lower than 1 mMa −1 ) should be favored, as these yield accurate elevations with a precision better than 1000 m.
In this paper we also make a case that the elevation at the time of initial exposure to cosmic rays should be carefully taken into account when computing burial ages.

General theory: determining altitudes from paired-cosmogenic nuclides
All of the calculations (i.e. simulations and data interpretations) presented here were performed using the parameters summarized in Table 1. Scaling factors were computed using the CREp calculator (https://crep .otelo .univ-lorraine .fr /#/) (Martin et al., 2017).
In a companion paper (Blard et al., 2019), we provide a Mat-lab© code, (Paleoaltitude.m), which allows paleoaltitudes to be computed from any dataset of cosmogenic nuclides (either the 26 Al-10 Be or the 10 Be-21 Ne pair), under both assumptions of continuous exposure and erosion (see Blard et al., 2019 for the mathematical description of the used algorithm). This code also plots data and simple exposure curves (ratio of two cosmogenic nuclides vs concentration of one cosmogenic nuclide) for variable altitudes in the same diagram. Fig. 1 shows the ratios of two radioactive cosmogenic nuclide (for the pairs 26 Al-10 Be and 10 Be-21 Ne) plotted against 10 Be concentration at different altitudes. These curves are defined by the equations determining the production of cosmogenic nuclides at the surface for different exposure durations with either (i) zero erosion (higher thick curve) or (ii) steady erosion (lower thin curve) ( Fig. 1) (Lal, 1991). In such a diagram, all samples that are exposed at constant depth and elevation define a region that has the shape of "banana". These plots were initially proposed by D. Lal, K. Nishiizumi and J. Klein (https://cosmognosis .wordpress .com / 2017 /03 /30 /upside -down -bananas/). If the exposure duration of a surface is shorter than 500 ka (for the 26 Al-10 Be and 10 Be-21 Ne pairs, Fig. 1), radioactive decay is negligible, and the isotope concentration ratios equal the production rates ratio of the two nuclides, meaning that these simple exposure curves overlap each other on the left part of the diagram (Fig. 1). In that case, it is only possible to determine a minimum elevation (Blard et al., 2019). In contrast, when the exposure time is longer (>500 ka, Fig. 1), the curves computed for different altitudes are no longer superimposed. The positions of these two-nuclides curves are thus altitude dependent ( Fig. 1): for similar exposure durations, samples standing at different elevations have similar 26 Al/ 10 Be (or 10 Be/ 21 Ne) ratios, but different cosmogenic nuclides concentrations. Although this property has been mentioned in few papers (e.g. Fig. 4 in Stone, 2000), this is often ignored in the literature. The positions of these simple exposure curves correspond to a two-equationstwo-unknowns problem, the two unknowns being the duration of exposure (or erosion) and the elevation. Therefore, measuring two nuclides in a steadily eroding, or non-eroding surface, can in theory allow the sample's exposure altitude to be determined (see section 2.2). The fact that a set of observations from surface samples should lie in this simple exposure region has already been Table 1 Parameters used to perform calculations (numerical simulation and data inversion) in this article, from (Balco et al., 2008;Chmeleff et al., 2010;Granger, 2006;Kelly et al., 2015;Kober et al., 2011;Lal, 1991;Lifton, 2016;Martin et al., 2017Martin et al., , 2015Nishiizumi et al., 1989Nishiizumi et al., , 2007Stone, 2000 Lal, 1991;Nishiizumi et al., 1989;Stone, 2000;Martin et al., 2017;Muscheler et al., 2005 Spatial and time-dependent scaling --Stone time dependent -Standard atmosphere -Atmospheric 10 Be VDM Fig. 1. Theoretical plots of the ratios of cosmogenic nuclides vs 10 Be concentrations for the (A) 26 Al-10 Be and (B) 10 Be-21 Ne isotopic systems. The ratios are plotted for continuous exposure (thick curves) and steady state denudation (thin curves) of the sampled surfaces, defining the so called "simple exposure curves". These curves were computed at 60 • latitude using the production parameters defined in Table 1. used to refine production rates (e.g. Kober et al., 2011), but this is the first time it is used to determine paleoelevations.

Physical and mathematical systematics of the paired-cosmogenic nuclides altimeter
The concentration of cosmogenic nuclides at or near the Earth's surface as a function of time, t, and depth below the surface x is described by the following equation (Lal, 1991): where N(x, t) is the nuclide concentration (at g −1 ), N(x, 0) the initial nuclide concentration, x (cm) is the depth below the surface. This depth x can be defined for any sort of covering material: rock, soil, ice, liquid water. P (at g −1 yr −1 ) is the cosmogenic nuclide production rate that is a function of x, and the spatial position of the sample, i.e. its altitude, latitude and longitude. λ is the radioactive decay constant (yr −1 ), μ = ρ/Λ is the cosmic ray atten-uation constant (cm −1 ), Λ being the attenuation length (g cm −2 ), ρ (g cm −3 ) the density of the covering material, ε (cm yr −1 ) the erosion rate of the surface.
For surface samples (x = 0) and under the assumption of no previous exposure or inheritance (N(x, 0) = 0), equation (1) can be simplified into: where P SLHL is the cosmogenic nuclide production rate normalized to Sea Level High Latitude (SLHL) (at g −1 yr −1 ) and f is the spatial scaling factor, which depends on the position of the sample (altitude, latitude and longitude). In this theoretical case, we assume that the altitude change is small over the exposure time t The scaling factor f accounts for the spatial variations in the cosmogenic nuclide production rates at the Earth's surface, due to the strong influence of the Earth's environment on cosmic par-ticles (e.g. Lal, 1991): with increasing elevation the atmospheric depth through which cosmic rays must travel before reaching the surface is smaller, yielding a higher cosmic ray flux. The impact of altitude is significant, since f increases exponentially with elevation: for each 1000 m of elevation gain, the production rate increases by a factor of two (Lal, 1991). To a lesser extent, f is also controlled by the intensity and the orientation of the Earth's magnetic field (Dunai, 2001;Lal, 1991), implying that latitude also affects the value of f : at the Equator, production rates are about two times lower than they are at the poles (Lal, 1991).
f is therefore a function of the sample altitude and latitude and may also in some cases include past variations of the magnetic field (e.g. Balco et al., 2008). If the latitude and the magnetic field variations at a sample site can be independently constrained, solving equation (2) for f allows the elevation of a sample to be constrained.
Equation (2) has a total of three unknowns: the scaling factor f , the surface erosion rate ε and the exposure time t. Thus, two extreme cases may be considered: i) No erosion or ii) steady-state erosion (t → +∞) (Lal, 1991): It is possible to reduce the degrees of freedom of equation (3a) and (3b) by combining two isotopic systems with different half-lives and, under favorable circumstances, solve for f . These conditions are satisfied for long exposure times (>500 ka) or low erosion rates (<1 mMa −1 ) and correspond to the case of nonsuperimposed simple exposure curves.
A companion methodology paper (Blard et al., 2019) provides a complete description of the different conditions that must be considered when solving for f .
i) The case of a continuously exposed, non-eroded surface In the zero-erosion case, equation (3a) can be solved for time by considering the combination of two nuclides with production rates P SLHL,1 and P SLHL,2 , respective decay constants λ 1 , λ 2 and concentrations N 1 , N 2 : In this case, there is no analytical solution, and equation (4) must be numerically solved to determine f . It is important to consider the analytical uncertainties attached to N 1 and N 2 to assess the ability to determine f from equation (4) with its appropriate uncertainties (detailed methods in Blard et al., 2019).

ii) Case of steady state erosion
In the case of steady erosion, equation (3b) leads to the following rather straightforward analytical solution, for a pair of radioactive cosmogenic nuclides: In theory, both equations (4) and (5) can be solved for f . If latitude and the past geomagnetic field variations can be independently constrained with reasonable confidence, f can in turn be used to determine the sample's altitude of exposure.
In practice, nuclide concentrations are only known through laboratory measurements and are thus affected by analytical uncertainties (or biases) that may lead to cases that are not solvable. In the next section 2.3, we present numerical simulations to compute the minimum duration of (paleo)-exposure required to solve the paleoelevation (i.e. the detection limit of the system), as well as the impact of exposure duration on the uncertainty in the computed elevation.

The detection limit of the paleoaltimeter and the impact of the exposure duration on the uncertainty
The mathematical conditions required to solve equations (4) and (5) are described in detail in the companion publication in MethodsX (Blard et al., 2019). Because of the analytical uncertainties attached to measurements of 10 Be, 21 Ne and 26 Al concentrations, it is possible that cases may be encountered where the system cannot be solved for elevation, or can only be used to determine minimal elevation with a large uncertainty, notably in the case of short exposure (<500 ka) or erosion rate To evaluate the influences of the exposure duration and the erosion rate on the detection limit, accuracy and the overall uncertainty of this paleoaltimetry method, we performed a numerical Monte Carlo simulation, with analytical uncertainties of 5% on each cosmogenic nuclide, for the pairs 26 Al-10 Be and 10 Be-21 Ne (Fig. 2).
The modeling shows that both 26 Al-10 Be and 10 Be-21 Ne have similar characteristics: longer exposure duration (or lower erosion rates) improve both the accuracy and the precision of the computed elevation. Exposure durations longer than 1 Ma or erosion rates lower than 1 mMa −1 guarantee the determination of accurate elevations, with absolute 1σ uncertainties that are smaller than 500 m. For exposure durations shorter than 500 ka, and erosion rates higher than 1 mMa −1 , the computed elevations are lower than the real altitudes of exposure. This systematic uncertainty results from the non-linearity of equations (4) and (5). In other words, for low nuclide concentrations, the two-nuclides curves are superimposed and the method thus only allows determination of a minimum elevation ( Fig. 1 and 2 and mathematical description in Blard et al., 2019). This bias increases when the exposure duration is shorter -or the erosion rates is higher -and may reach several hundreds of meters (Fig. 2). Using nuclide pairs with a cosmogenic nuclide having a shorter half-life (e.g. 36 Cl/ 10 Be) would improve the method sensitivity for higher erosion rates and shorter exposure times.
In summary, this Monte Carlo simulation shows that the accuracy and the precision of the method are much better for long exposure durations ( 100 ka) or low erosion rates ( 1 m Ma −1 ). However, shorter exposure durations may also provide useful information in that knowing the minimum elevation at which a terrain was exposed could also be a useful outcome in certain geological contexts.

Bias due to poor knowledge of the degree of surface preservation: is the surface eroded or not?
It is difficult, often impossible, to tell a posteriori whether a paleosurface has undergone continuous exposure or steady erosion. The unknown state of the paleosurface may induce a systematic uncertainty in the reconstructed elevation of that surface (independent of analytical uncertainties that are random), since the positions of simple exposure curves differ slightly for a continuous exposed surface and a steadily eroded surface ( Fig. 1 and 3A). In other words, if a preserved surface is wrongly assumed to have and erosion rate (C and D) on the accuracy and the precision of the paleoaltimeter, for the 26 Al-10 Be and the 10 Be-21 Ne pairs. These simulations were performed using the parameters presented in Table 1, assuming analytical uncertainties of 5%, at a latitude of 60 • , using the Lal/Stone (Stone, 2000) time independent scheme, at 0, 2000 and 4000 m. For each time step, 500 random draws (5000 draws for each erosion rate step) were made for nuclides 1 and 2, assuming that (N1, σ 1) and (N2, σ 2) follow normal distributions, σ 1 and σ 2 being the analytical uncertainties (5% here). Next, the randomly picked cosmogenic nuclides pairs were solved for elevation using equations (4) and (5), in the case of continuous exposure and steady-state erosion, respectively. Central curves represent the medians (50%) and the means, error envelopes are bounded by the 0.16 and 0.84 confidence intervals (1σ ). been eroded, this will lead to overestimation of the actual elevation of exposure (Fig. 3A). Conversely, assuming continuous exposure when the surface has actually undergone erosion will lead to underestimation of elevation.
We estimated this bias by modeling the difference in the reconstructed altitudes for a pair of surface nuclide concentrations by considering both the continuous exposure and the steady erosion scenarios (Fig. 3B). This bias results from the difference: (4) and (5)). Fig. 3B shows how this systematic uncertainty in the reconstructed paleoaltitude is directly linked to the unknown state of the paleosurface for both the 10 Be-26 Al and 10 Be-21 Ne pairs.
The modeling highlights that for exposure durations shorter than 100 ka, (or erosion rates higher than 1 m Ma −1 ), this bias may reach 900 m for a surface exposed at sea level and up to 1900 m for exposures at 6000 m (Fig. 3B). However, for surfaces that have been exposed long enough ( 100 ka) (or that have been slowly eroded 1 m Ma −1 ) this systematic uncertainty in the reconstructed paleoaltitude decreases significantly, until it is null for surfaces with radioactively-saturated 10 Be, 26 Al and 21 Ne concentrations.
Longer exposures (or lower erosion rates) thus significantly improve both the accuracy and the precision of the paleoaltimetric method.

Why muons can be (almost) safely neglected here
In the equations described in section 2.2, P SLHL only refers to the production rate due to spallogenic high-energy neutrons. This means that the different attenuation lengths of neutrons and muons (Groom et al., 2001) are ignored here. This is justified because the favorable surfaces for paleoaltimetry are either noneroding, and thus the muogenic production at depth does not affect the surface nuclide concentration, or are eroding very slowly, meaning that the deep production is negligible compared to the total surface concentration since the muon-produced nuclides will have decayed before reaching the surface. Stable 21 Ne represents a possible exception however, since muogenic 21 Ne produced at depth does not decay on its way to the surface and may therefore make a greater contribution to the total surface concentration than 10 Be and 26 Al do. Neglecting muogenic 21 Ne in a slowly eroding surface may therefore cause underesti- ). This systematic uncertainty is plotted against the 10 Be concentration, a proxy of increasing exposure duration (or decreasing erosion rate). The wrong assumption that a preserved surface has been eroded will lead to overestimation of the actual elevation. Conversely, assuming continuous exposure when the surface has actually undergone erosion will lead to underestimation of elevation. Plots show the two nuclides pairs: 10 Be-26 Al (solid line) and 10 Be-21 Ne (dashed line). These uncertainties are plotted for a surface exposed at constant elevations of 0, 2000, 4000 and 6000 m.a.s.l. mation of the 10 Be/ 21 Ne ratio, and thus underestimation of the correct elevation. This potential bias must be kept in mind when using this pair of nuclides as a paleoaltimeter. Ideally, the muogenic production term should be added to the equation describing 21 Ne production.
2.6. Paleoaltimetry using fossil exposed surfaces 2.6.1. Approach Continuously exposed surfaces only offer a long-term integrated signal that includes the most recent elevation changes (Blard et al., 2006). The altimetry method presented in section 2.2 may also be used to reconstruct the altitude of an ancient surface that underwent continuous exposure or steady erosion in the past and that has since been buried for a known period of time. The main interest of paleosurfaces is to provide a snapshot of the paleoelevation history of a massif at a particular moment in its uplift, or subsidence, history. In such a case, the present-day measured concentrations of radioactive nuclides ( 10 Be, 26 Al) in a completely shielded paleosurface must be corrected for the radioactive decay that occurred since burial (Fig. 4).
A number of requirements need to be fulfilled to derive the most accurate and precise paleoaltitude possible from a paleo-exposed surface: 1) like non-buried surfaces, the paleosurfaces need to have been exposed for long enough during the past ( 100 ka), or to have been affected by low erosion rates ( 1 m Ma −1 ); 2) the surface must have been rapidly and deeply Fig. 4. A) Schematic diagram of an ideal paleosurface for cosmogenic paleoaltitude reconstructions, where a paleosurface is rapidly buried under a datable layer. B) Principle of paleoaltitude evolution of a surface: the surface is exposed during a period of different elevation and has been buried since. The nuclide concentration on this surface reflects the elevation prior to the altitude change. C) Simple exposure curves evolution of the paired 26 Al and 10 Be concentrations of a buried surface. D) Simple exposure curves evolution of the paired 21 Ne and 10 Be concentrations of a buried surface. buried (ideally below tens of meters, Fig. 6; Blard et al., 2006) so that there is negligible post-burial nuclide accumulation; and 3) the burial age should be independently estimated. Under these conditions, it is possible to correct the cosmogenic nuclide concentrations measured on a paleosurface, N measured , for radioactive decay that has occurred since burial, t burial , to reconstruct the initial nuclide concentration, N initial : Equation (6) can then be combined with equations (4) and (5) to reconstruct the paleoaltitude of the surface. In the case of stable nuclides, such as 21 Ne, the measured concentration remains constant during burial and a radioactive decay correction is not necessary.
Unlike previous approaches to cosmogenic nuclide paleoaltimetry (Blard et al., 2005;Libarkin et al., 2002), this paired-nuclides method is in theory applicable for both steadily eroding and continuously exposed surfaces. As the preservation state of a surface is sometimes difficult to assess (section 2.4), this is therefore a major advantage of the approach proposed in this article. The method allows investigation of a much broader range of geomorphic objects that have been exposed over a large time-span in the past. Volcanic environments are likely good candidates for such paleoaltitude reconstructions as volcanic flows provide rapid burial of the underlying surfaces and can be dated over timescales of several million years using absolute radiometric methods such as K-Ar or Ar-Ar.

Method range
The method range directly depends on the analytical detection limit of the analyzed cosmogenic nuclides, and thus, on the combined influences of the paleoexposure duration and the age of the burial of the fossil exposed surface. Considering the current analytical detection limits of 10 Be, 21 Ne and 26 Al (∼5 × 10 3 , ∼10 5 and ∼10 4 at g −1 , respectively), the oldest buried surface that may be confidently used to apply this paleoaltrimetry method is ∼8 Ma in the case of the 26 Al-10 Be pair and ∼12 Ma in the case of the 10 Be-21 Ne pair (Fig. 4).

Integration time and response time after rapid uplift or subsidence scenarios
Given that this paleoaltimetry method requires relatively long exposure episodes (>500 ka), it is probable that the elevation has changed during exposure, meaning that the assumption of a time independent scaling factor is not valid: df /dt = 0.
The integration time t int can be defined as the mean age of a given cosmogenic nuclide contained in a rock sample. t int is defined by the following equation, where N mes (at g −1 ) is the measured concentration of the nuclide N and P (at g −1 yr −1 ) its local production rate: t int is different from the actual exposure time t exp that can be derived from equation (3a): Two extreme cases can be distinguished: i) If exposure time t exp is much shorter than the half-life of the nuclide (t exp 1/λ), then t int = t exp /2, ii) If t exp is much larger than the half-life of the nuclide (t exp 1/λ), then t int = 1/λ.
In the practical case of the paired-nuclide altimetry, t int is determined by the cosmogenic nuclide having the shortest half-life.
This integration time constrains the response time of the method after an altitudinal change. We performed several numerical tests with different "staircase" uplift and subsidence scenarios. These simulations shown on Fig. 5 constrain the method reactivity. Although these extreme staircase scenarios are not encountered in natural settings, they represent benchmarks and provide limit values for the response time of the altimeter.
These numerical simulations (Fig. 5) show that the response time is dependent on the half-lives of the nuclides considered: in the uplift case, the recorded elevation is ±10% similar to the correct altitude after 1.5 Ma of exposure for the 26 Al-10 Be pair (Fig. 5A), while it requires ∼3 Ma for the 10 Be-21 Ne pair (Fig. 5B).
This modeling also shows that the response time of the registered altitude is shorter in the case of uplift than in the case of subsidence (Fig. 5). This is due to the exponential increase in production rate with elevation, which gives a larger weight to the nuclides produced at high elevation.
3. Testing the method against the existing 10 Be-21 Ne-26 Al dataset from Atacama

Data description
In order to evaluate the accuracy of this new paleoaltimetric method as well as its overall uncertainty in real conditions, we tested this approach on an existing dataset of samples, in which multiple cosmogenic nuclides have been measured ( 26 Al and 10 Be, or 21 Ne and 10 Be). For this, we selected datasets from a geological setting that fulfilled the following criteria: i) samples have presumably been exposed at the surface for a sufficiently long and continuous duration ( 100 ka); ii) rocks have remained at the same sampling elevation without any significant elevation change during the exposure history; and iii) the dataset covers a quite large altitudinal range, from sea-level to several thousand of meters. For this initial test, we chose the Atacama desert, a very arid region of the Western Tropical Andes that meets the required geological and geomorphological criteria: i) rainfall is extremely low (<10 mm w.e.yr −1 ), making long exposure (>1 Ma) possible (Dunai et al., 2005;Ritter et al., 2018b), ii) several lines of evidence suggest that no altitudinal change greater than 1000 m has occurred since 5 Ma (Garzione et al., 2008;Kar et al., 2016), and iii) the steep west-east topographic gradient gives access to an altitudinal range from 0 to 5000 m.
We consider here all published datasets that have reported measurements of cosmogenic nuclides pairs ( 26 Al-10 Be and 10 Be-21 Ne) from surface samples collected in the Atacama desert: (Kober et al., 2007;Nishiizumi et al., 2005;Placzek et al., 2010;Ritter et al., 2018aRitter et al., , 2018b. Table S1 presents the main characteristics of these samples (nature, lithology, sampling coordinates). The dataset includes a total of 68 samples for the 26 Al-10 Be pair and 43 samples for the 10 Be-21 Ne pair. All of the objects sampled have a quartz rich lithology and variable geomorphologic characteristics: bedrock and lake shoreline samples, and detrital objects of different sizes and sources: clasts, cobbles and boulders deposited by alluvial or gravitational processes.

Comparison between calculated and present-day elevations
As a first screening, we plotted the data in the two isotopes diagram, along with the simple exposure curves for different elevations (0, 1500, 3000 and 4500 m). In Fig. 6, the sampling elevations are represented by the color of the ellipses. To allow a proper comparison, all cosmogenic nuclides concentrations were scaled to the same latitude of 20 • using the time independent model of (Stone, 2000) and the standard atmosphere model (N.O.A.A., 1976).
The plots indicate a good first order agreement between the sampling altitudes, shown by the colors of the ellipses, and the elevations deduced from the cosmogenic nuclides, which are represented by the relative position of these ellipses compared to those of the exposure curves. In the graphs, no sample, except one, plots in the upper forbidden zone -the area above and to the right of the curves -strongly suggesting that none of the samples originated from a higher elevation and that the cosmogenic nuclide production rates are well constrained. On the other hand, a few samples plot below or on the left of their respective exposure curves, suggesting that these samples have been affected by periods of burial or recently exhumed.
Importantly, the majority of the samples are located on their corresponding elevation-curves, suggesting that the dataset was not affected by any significant elevation change since the initiation of exposure.
For a more precise comparison, we also inverted all of the data to compute paleoelevations using the altimetry method, following the mathematical approach described in (Blard et al., 2019). Ta- ble S2 displays the detailed results and Fig. 7 is a plot of these computed elevations vs sampling elevations. Since we do not have any a priori knowledge of the amount of erosion that affected these samples, the two extreme cases of continuous exposure and steady-state erosion were considered ( Fig. 7A and B). Sample exposures that did not satisfy the conditions to compute elevations are not plotted on Fig. 7. Among the 43 21 Ne-10 Be samples, elevations could not be calculated for only 2 samples. In the case of the 68 26 Al-10 Be samples, mean altitudes could not be calculated for 18 samples in the case of continuous exposure, and 10 in the case of steady-state erosion. This difference mainly results from the shorter integration times of the 10 Be-26 Al dataset (0.55 ± 0.28 Ma) compared to those of the 21 Ne-10 Be dataset (1.22 ± 0.60 Ma). Fig. 7 plot shows that there is a first order good agreement between the computed elevations and the sampling elevations. The two extreme scenarios of continuous exposure and steady-state erosion yield quite similar results, indicating that the method is robust, whatever our knowledge of the surface preservation state. Best-fit regression lines were computed, along with their parameter uncertainties at the 2-sigma confidence level.
In the case of the 26 Al-10 Be dataset, the best-fit regression curves lie slightly below the 1:1 line, both for the continuous exposure ( y = (0.78 ± 0.31) · x + (10 ± 600), R 2 = 0.50, p-value < 6 ×10 −6 ) and the steady-erosion case ( y = (0.70 ±0.39) · x +(670 ± 770), R 2 = 0.39, p-value < 10 −3 ). This slight divergence with the 1:1 line is due to a group of 6 samples that today stand between 2300 and 3800 m, but that yielded computed elevations 1000 to 2000 m below their present position. This bias may result from unrecognized recent exhumation of samples that spent most of their exposure histories more than 50 cm below the rock surface (see section 4.1 below). However, if these 6 points are excluded, the agreement between the sampling and the computed elevation is excellent, indicating that the method is robust and accurate, even when using a large dataset that was not specifically sampled for paleoaltimetric purposes.
The 21 Ne-10 Be dataset yields best fits that are statistically in quite good agreement with the 1:1 line: y = (1.12 ± 0.21) · x − (800 ± 610), R 2 = 0.83, p-value < 2 × 10 −13 in the case of continuous exposure, and y = (1.42 ± 0.21) · x − (760 ± 540), R 2 = 0.83, p-value < 6 × 10 −16 , in the case of steady-erosion. The samples that are not aligned on the 1:1 line are those from the Quillagua-Llamara Soledad paleoake shorelines (Ritter et al., 2018a), which are located at low elevations (∼1000 m) and yield underestimated computed elevations (ranging from −1500 to 500 m). Several possibilities might be proposed to explain the fact that many of the shorelines samples of Ritter et al. (2018a) yielded computed altitudes lower than their present-day elevations: i) A significant part of the exposure may have started and occurred below the lake surface, or below soil cover that has been removed during the upper Pleistocene; ii) We did not consider the muogenic production in the equation. As muon production is greater than spallation below a depth of several meters, it is plausible that neglecting this process has lowered the 10 Be/ 21 Ne ratios after the decay of the muogenic 10 Be produced at depth, leading to erroneously low elevations. Finally, it is interesting to note that the detrital samples with the longest exposure (∼20 Ma) and integration times (∼2 Ma) also yield computed elevations that are almost identical to (or slightly lower than) their present-day sampling elevations (Table S2). This shows that these samples, despite their detrital origin, have cosmogenic-nuclides inventories that do not result from exposures at higher elevations. Hence, the exposure times computed using the present-day elevations are probably correct and are not overestimated. Fig. 6. Plots of cosmogenic nuclide pairs (A -26 Al-10 Be and B -10 Be-21 Ne) measured in surface samples from the Atacama region (data from Kober et al., 2007;Nishiizumi et al., 2005;Placzek et al., 2010;Ritter et al., 2018aRitter et al., , 2018b. There are 68 samples for the 26 Al-10 Be pair and 43 samples for the 21 Ne-10 Be pair. This plot was realized using the production parameters given in Table 1. All cosmogenic nuclides concentrations have been scaled to the same latitude of 60 • using the time independent model of Stone (2000). Ellipses are plotted for 68% confidence intervals. The color of each ellipse represents the sampling elevation. Simple exposure curves are plotted for elevations of 0, 1500, 3000 and 4500 m. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Impact of the exposure depth: the equivalence between atmospheric and rock depth
The theory that describes how the altitude of exposure controls the position of the simple exposure curves (Section 2, Equations (3a) and (3b)) is also valid in the case of different depths of exposure below the rock surface (Gosse and Phillips, 2001). Considering that i) the attenuation lengths of fast neutrons in gas and in soils are similar (Gosse and Phillips, 2001), and ii) that soils have a density about 2000 times higher than the density of atmosphere at sea level (assuming a soil density of ∼2.4 g cm −3 , and an atmosphere density of ∼1.2 × 10 −3 g cm −3 at sea level), 1000 m of air is there-fore equivalent to a soil thickness of 50 cm. Given this equivalence, it is crucial that the thickness and the density of the overburden are well known if a rock is sampled below a paleo-exposed surface for paleoaltimetric investigation. If an unknown process has recently exhumed the rock, there is a risk that the actual altitude of exposure will be significantly underestimated using the paleoaltimetry method (Fig. 3). This important aspect should be kept in mind during sampling and any evidence that suggests the sporadic presence of loess, soils or any cover should be considered before sampling a paleosurface on a fresh outcrop. The possibility of landslides should likewise be carefully assessed.
A potential depth-meter. If the altitude of exposure is wellknown, the position of the simple exposure curve may also be turned into a sensitive and precise depthmeter (Gosse and Phillips, 2001). If exposure is sufficiently long ( 100 ka), the paleoaltime- Fig. 7. Sampling against elevations calculated with the paired cosmogenic nuclides method considering A) continuous exposure and B) steady-state erosion, for the whole Atacama dataset (data from Kober et al., 2007;Nishiizumi et al., 2005;Placzek et al., 2010;Ritter et al., 2018aRitter et al., , 2018b. Since no systematic difference is observed between bedrock, boulders and cobbles (Table S2), the different types of sample are not differentiated in this plot. The best-fit regression and the 1:1 lines are shown. Integration times are represented by the color-scale. Uncertainties of the best-fit parameters are two sigmas; uncertainty envelopes are shown in red. (Only those samples for which it was possible to calculate elevation, positive and negative error bars are shown.) try method has a 1σ uncertainty of ∼500 m of atmosphere for the 26 Al-10 Be and the 10 Be-21 Ne pairs (Section 2.3). Thus, in such cases, the method is potentially sensitive enough to measure the depth of exposures with a precision of ∼25 cm of rock/soil. (Hidy et al., 2018) used the 26 Al-10 Be couple to calculate a thickness of ∼80 cm of loess cover over a paleosol exposed during more than 1 Ma in Yukon, Canada. For exposures that occur below water (density of 1 g cm −3 for water), the uncertainty in the measured water depth would be slightly larger (∼50 cm). Note that the minimum required exposure time will be significantly reduced (and the precision improved) using a radioactive cosmogenic nuclide that has a relatively short half-life, such as 14 C (half-life of 14 C = 5730 years). Although presenting the limits of this specific aspect of the method is beyond the scope of this article, it is worth noting that the 14 C-10 Be pair has the best potential for depth determinations. It has notably been used to constrain the depth of partial snow shielding in the Gotthard Pass area in the Swiss Alps (Hippe et al., 2014). Several intriguing questions can thus be addressed with this depthmeter, such as measuring the mean depth of landslides or the thickness of paleocovers of any nature, for example soil, ash, loess, snow, ice or vegetation. Hydrological studies might even be conducted with this method, as it theoretically allows measurement of the (paleo)water depth of a (paleo)lake.

Impact of the elevation of exposure on the accuracy of burial ages
The burial age dating method has been widely used to place important geochronological constraints on several major problems in Earth sciences (e.g. Granger et al., 2015;Sartégou et al., 2018). One of the main strengths of this method is due to the fact that the preburial 26 Al/ 10 Be (or 10 Be/ 21 Ne) ratio can often be considered as "independent of latitude and altitude" (Granger and Muzikar, 2001). However, as illustrated by Fig. 8, in the case of long preburial exposure ages (>100 ka) or low erosion rates (<1 m Ma −1 ), altitude could affect the initial preburial 26 Al/ 10 Be (or 10 Be/ 21 Ne) ratio and, when it is known, it should therefore be taken into account in the burial age calculation (Equations (9) and (10)). In the Fig. 8. Example showing the impact of the altitude of exposure on the accuracy of burial ages, in the case of the 26 Al-10 Be system. In the given example, assuming an initial 26 Al/ 10 Be ratio similar to the production ratio leads to an overestimate the actual burial age. For long preburial exposure durations (>100 ka) or low erosion rates (<1 m Ma −1 ), it is important to take into account the elevation of exposure to compute the preburial 26 Al/ 10 Be (or 10 Be/ 21 Ne) ratio and, thus, to obtain accurate burial ages. specific case of the burial dating of cave sediments, an accurate calculation should use the elevation of the quartz-rich watershed, rather than the altitude of the cave. Fortunately, in the majority of cases, preburial ratios are in practice similar to the production ratios, implying that many burial ages computed with this approximation are not affected by this potential elevation-related bias.
For a material that has been exposed under steady-erosion conditions before burial, the equation that must be used to compute the burial age t burial is: If nuclides 1 and 2 are both radioactive, equation (9) has no analytical solution and must be numerically solved to determine t burial . For the particular case of a cosmogenic nuclides pair with only one radioactive isotope, such as 10 Be-21 Ne, equation (9) simplifies and can be analytically solved: A Matlab© code, (Burial.m), computing burial ages taking into account the altitude of the preburial exposure episode is available in a MethodsX companion paper (Blard et al., 2019).

Concluding remarks
• In this article, we propose a new means of determining paleoelevations using pairs of in situ cosmogenic nuclides that have different half-lives. Positions of cosmogenic nuclide simple exposure curves are altitude-dependent, a property that may be used to constrain the elevation of exposure. In practice, given current analytical capabilities, the best nuclide pairs for this purpose are 26 Al-10 Be and 10 Be-21 Ne in quartz, which have respective ranges of 6 Ma and 12 Ma, respectively.
• Both 26 Al-10 Be and 10 Be-21 Ne systems require minimum equivalent exposure durations that are longer than 100 ka (ideally > 500 ka, or erosion rates lower than 1 mMa −1 ) to yield accurate elevations, with a typical 1σ uncertainty of ∼500 m. • If the exposure duration is shorter than 500 ka (or if the erosion rate is less than <1 m Ma −1 ), the method allows minimum elevations of exposure to be calculated.
• A poor knowledge of the preservation state of the surface (i.e. the absence or presence of erosion) may induce a bias on the computed elevation of less than 1000 m at sea level.
• We applied the method to bedrock and detrital objects from the Atacama desert (Andes). Results show a good agreement between the computed altitudes and the present-day sampling elevations.
• Given that the majority of erosion rates on Earth are greater than 1 m Ma −1 , the method may not be applicable to obtain mean elevations in many regions. However, if erosion rates are greater than 1 mMa −1 , the method may be applied to detrital material to determine the minimum basin-averaged altitude -a constraint that could be useful in addressing several geodynamics-related problems.
• An intriguing direct application that can be derived from this paleoaltimetry method is depthmetry, i.e. the measurement of soil, snow, ice or water thickness during exposure.
• Since the positions of a two-nuclide curves is altitude dependent, then the pre-burial cosmogenic nuclide ratio of the buried material may be affected by the altitude of exposure. This must be considered before applying the burial dating method (Granger and Muzikar, 2001) in order to avoid any potential inaccuracy. This is particularly true in the case of watersheds that experience low erosion rates (<1 m Ma −1 ) and that span large elevation ranges (>1000 m).