Matt D Tomkins
email: matthew.tomkins@manchester.ac.uk

Jason M Dortch

Philip D Hughes

Jonny J Huck

Raimon Pallàs

Ángel Rodés

James L Allard

Andrew G Stimson

Didier Bourlès

Vincent Rinterknecht

Vincent Jomelli

Laura Rodríguez-Rodríguez

Ramon Copons

Iestyn D Barr

Christopher M Darvill

Thomas Bishop

J.-P Avouac

Moraine crest or slope: An analysis of the effects of boulder position on cosmogenic exposure age

Keywords:

Terrestrial cosmogenic nuclide dating of ice-marginal moraines can provide unique insights into Quaternary glacial history. However, pre-and post-depositional exposure histories of moraine boulders can introduce geologic uncertainty to numerical landform ages. To avoid geologic outliers, boulders are typically selected based on their depositional context and individual characteristics but while these criteria have good qualitative reasoning, many have not been tested quantitatively. Of these, boulder location is critical, as boulders located on moraine crests are prioritised, while those on moraine slopes are typically rejected. This study provides the first quantitative assessment of the relative utility of moraine crest and moraine slope sampling using new and published 10 Be and 36 Cl ages (n = 19) and Schmidt hammer sampling (SH; n = 635 moraine boulders, ∼19,050 SH R-values) in the northern and southern Pyrenees. These data show that for many of the studied moraines, the spatial distribution of "good" boulders is effectively random, with no consistent clustering on moraine crests, ice-proximal or -distal slopes. In turn, and in contrast to prior work, there is no clear penalty to either moraine crest or moraine slope sampling. Instead, we argue that landform stability exerts a greater influence on exposure age distributions than the characteristics of individual boulders. For the studied landforms, postdepositional stability is strongly influenced by sedimentology, with prolonged degradation of matrix-rich unconsolidated moraines while boulder-rich, matrix-poor moraines stabilised rapidly after deposition. While this pattern is unlikely to hold true in all settings, these data indicate that differences between landforms can be more significant than differences at the intra-landform scale. As ad hoc assessment of landform stability is extremely challenging based on geomorphological evidence alone, preliminary SH sampling, as utilised here, is a useful method to assess the temporal distribution of boulder exposure ages and to prioritise individual boulders for subsequent analysis.

Introduction

Ice-marginal moraines are classic features of glaciated mountain ranges and are prominent terrestrial records of glacial history [START_REF] Hallet | Surface dating of dynamic landforms: young boulders on aging moraines[END_REF]. By constraining the timing of moraine deposition, it is possible to reconstruct the growth and decay of glaciers and ice sheets through the Quaternary and the palaeoclimatic drivers of glacial cycles. Recent developments in terrestrial cosmogenic nuclide (TCN) dating have transformed our understanding of Quaternary glaciations by permitting direct analysis of the fragmentary glacial stratigraphic record (Zreda and Phillips, https://doi.org/10.1016/j.epsl.2021.117092 1995). Despite this progress, TCN dating can be complicated by geologic processes which result in pre-or post-depositional exposure of rock surfaces and which account for apparent TCN ages that pre-or post-date the assumed age of the landform [START_REF] Applegate | Modeling the statistical distributions of cosmogenic exposure dates from moraines[END_REF]. Of these, post-depositional erosion, exhumation and shielding have been shown to profoundly influence TCN age distributions [START_REF] Briner | Cosmogenic exposure dating of late Pleistocene moraine stabilization in Alaska[END_REF][START_REF] Zech | Evidence for long-lasting landform surface instability on hummocky moraines in the Pamir Mountains (Tajikistan) from 10 Be surface exposure dating[END_REF][START_REF] Heyman | Too young or too old: evaluating cosmogenic exposure dating based on an analysis of compiled boulder exposure ages[END_REF][START_REF] Stübner | Middle and Late Pleistocene glaciations in the southwestern Pamir and their effects on topography[END_REF][START_REF] Chevalier | Deciphering old moraine age distributions in SE Tibet showing bimodal climatic signal for glaciations: marine isotope stages 2 and 6[END_REF].

To avoid geologic outliers, researchers select samples based on the depositional context and characteristics of individual surfaces. Previous studies have advocated sampling:

• boulders on moraine crests or on flat, stable surfaces [START_REF] Gosse | Precise cosmogenic 10 Be measurements in western North America: support for a global Younger Dryas cooling event[END_REF],

• the tallest boulders, to minimise the likelihood of postdepositional shielding [START_REF] Heyman | Boulder height -exposure age relationships from a global glacial 10 Be compilation[END_REF],

• the largest boulders or boulders embedded in the moraine matrix [START_REF] Ivy-Ochs | Cosmogenic nuclides and the dating of lateglacial and early Holocene glacier variations: the Alpine perspective[END_REF], to minimise the likelihood of post-depositional instability,

• well-rounded boulders which preserve evidence of glacial transport [START_REF] Darvill | Geomorphology and weathering characteristics of erratic boulder trains on Tierra del Fuego, southernmost South America: implications for dating of glacial deposits[END_REF], to minimise the likelihood of pre-depositional exposure. However, while these criteria have good qualitative reasoning, many have not been tested quantitatively. In turn, further work is required to test existing criteria for sample selection and to develop quantitative methods which minimise the effects of geologic processes [START_REF] Dortch | Timing and climatic drivers for glaciation across semi-arid western Himalayan-Tibetan orogen[END_REF][START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF]. These developments have the potential to significantly improve the robustness of TCN datasets and the chronological utility of the moraine record [START_REF] Applegate | Improved moraine age interpretations through explicit matching of geomorphic process models to cosmogenic nuclide measurements from single landforms[END_REF]. Within this context, this paper focuses on a fundamental component of TCN sample selection; the effect of moraine crest sampling on boulder exposure age. Of the above criteria, boulder location is critical, as boulders on moraine crests are prioritised, while those on moraine slopes are typically rejected, irrespective of their individual characteristics.

While this approach is qualitatively sound, early numerical models of moraine evolution predicted the greatest groundlowering at moraine crests [START_REF] Hallet | Surface dating of dynamic landforms: young boulders on aging moraines[END_REF][START_REF] Putkonen | Accuracy of cosmogenic ages for moraines[END_REF] with a period of maximum instability as glaciers retreat and as oversteepened ice-proximal slopes erode and stabilise [START_REF] Porter | 36 Cl dating of the classic Pleistocene glacial record in the northeastern Cascade Range, Washington[END_REF]. However, moraines continue to degrade through time as a function of moraine height and sedimentology [START_REF] Putkonen | Accuracy of cosmogenic ages for moraines[END_REF][START_REF] Putkonen | Landscape evolution degrades the geologic signature of past glaciations[END_REF][START_REF] Schaller | Quantifying glacial moraine age, denudation, and soil mixing with cosmogenic nuclide depth profiles[END_REF], as diffusive processes remove finegrained material from moraine crests and deposit material at the base of moraine slopes [START_REF] Applegate | Modeling the statistical distributions of cosmogenic exposure dates from moraines[END_REF]. Over time, these processes drive exhumation of boulders which have been shielded from cosmogenic exposure. In turn, the age distribution of moraine crest boulders may primarily reflect an initial stabilisation phase (∼1 ka; [START_REF] Briner | Cosmogenic exposure dating of late Pleistocene moraine stabilization in Alaska[END_REF][START_REF] Dortch | Beryllium-10 surface exposure dating of glacial successions in the Central Alaska Range[END_REF], modified by the ongoing process of moraine degradation, rather than the timing of initial moraine deposition. In contrast, slope diffusion models and lichenometric methods predict relative stability on moraine slopes [START_REF] Hallet | Surface dating of dynamic landforms: young boulders on aging moraines[END_REF][START_REF] Putkonen | Degradation of unconsolidated quaternary landforms in the western North America[END_REF], but these are rarely sampled for TCN, in part due to the perceived risk that boulders may rotate, shift or roll throughout the lifetime of the moraine. This dichotomy between model predictions and sampling procedures raises a fundamental and currently unanswered question: should moraine crests or moraine slopes be prioritised in TCN sample selection?

To address this uncertainty, we utilise 19 new and published 10 Be and 36 Cl TCN ages and 635 Schmidt hammer calibratedexposure ages (SH; 19,050 SH R-values) from ice-marginal moraines in the northern and southern Pyrenees. Weathering-based analyses are utilised here to enable intensive sampling of boulders across the moraine surface, with results verified against independent TCN ages. In total, these data provide the first quantitative assessment of the relative utility of moraine crest and moraine slope sampling.

Methods

Moraine selection

Six moraines of varying age, geomorphology and sedimentology were selected in the Pyrenees (Fig. 1); a mountain range which was extensively glaciated during Pleistocene cold stages (see Fig. 1F; [START_REF] Calvet | Recent advances in research on quaternary glaciations in the Pyrenees[END_REF]. Moraines were selected to encompass the primary deglaciation phases of the Pyrenees since the global Last Glacial Maximum (LGM) and all feature large populations of quartz-rich granitic moraine boulders, sourced from Axial zone granite outcrops in the Arànser, Gave de Pau and Noguera Rigaborçana glaciated valleys (Fig. 1E), and which are suitable for 10 Be dating. While this focused approach does not comprise all moraine types or depositional settings, these sites do encompass a range of moraine types commonly found in cirque and valley landsystems and which are often priority targets for TCN dating (i.e. ≤ LGM).

Selected moraines include both left (north) and right (south) latero-frontal moraines in the Arànser catchment, Cerdanya (Fig. 2A). These moraines are matrix-rich (matrix-supported), steep-sided (30-40 • ), heavily forested (Mountain pine: Pinus uncinata), and record the maximum ice extent of the Arànser glacier during the Würmian cold stage (110 -11.7 ka;[START_REF] Calvet | Recent advances in research on quaternary glaciations in the Pyrenees[END_REF]. The right latero-frontal moraine has previously been dated using 36 Cl (n = 2; [START_REF] Palacios | Maximum extent of Late Pleistocene glaciers and last deglaciation of La Cerdanya mountains[END_REF]. To supplement these data, a further 10 boulders were selected for 10 Be analysis (Table 1). Methods used for sample preparation, 10 Be measurement and exposure age calculation are detailed in the Supplementary Information.

On the north side of the Pyrenees, matrix-rich lateral moraines were selected in the Gave de Pau catchment (Fig. 2B). At least two neighbouring (∼60 m) but distinct lateral moraine ridges have been identified (Soum d'Ech moraines; Fig. 1D), with the outer moraine previously dated using 10 Be (n = 4; [START_REF] Rodés | La última deglaciación en los pirineos: de superficies de exposición mediante 10 Be, y modelado numérico de paleoglaciares[END_REF]. As at Arànser, these moraines likely correspond to the Würmian MIE but their distinctive morphologies (multiple nested ridges vs. a single large moraine) likely reflects a topographic control on moraine deposition (open topography vs. confined valley; [START_REF] Barr | A review of topographic controls on moraine distribution[END_REF][START_REF] Palacios | Maximum extent of Late Pleistocene glaciers and last deglaciation of La Cerdanya mountains[END_REF].

On the south side of the Pyrenees, and in the Val de Molières catchment of the Noguera Rigaborçana, sampled sites include the boulder-rich (clast-supported), matrix-poor Outer Pleta Naua terminal moraine (Fig. 3B,1B), previously assigned to Greenland Stadial 1 based on 10 Be (n = 3; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF], and the Tallada cirque moraine (Fig. 3A,1A), which consists of a single sharpcrested, arcuate terminal moraine with two minor ice-proximal ridges. Both the Outer Pleta Naua and Tallada moraines are primarily composed of wedged angular boulders, with little or no sediment matrix. Although undated, the Tallada moraine is assumed to be late-Holocene in age, as evidenced by minimal boulder weathering [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF], while its elevation (∼2400 m), topographic setting (enclosed cirque; ∼0.16 km 2 ) and aspect (NNE) likely contribute to the inter-annual preservation of a small snowfield (∼0.03 km 2 ). These factors may have enabled glacier growth or re-advance during more recent climatic periods (e.g. the Little Ice Age).

Sampling approach

To investigate the depositional and post-depositional histories of these moraines, glacial boulders were selected to cover the entire moraine surface, including the moraine crest (C), the inner (11.7 -110 ka;[START_REF] Calvet | Recent advances in research on quaternary glaciations in the Pyrenees[END_REF]. Locations and sample names for TCN dated boulders are shown (white circles; [START_REF] Rodés | La última deglaciación en los pirineos: de superficies de exposición mediante 10 Be, y modelado numérico de paleoglaciares[END_REF][START_REF] Palacios | Maximum extent of Late Pleistocene glaciers and last deglaciation of La Cerdanya mountains[END_REF]. In (A), the locations of the proximal Fornell (F) and Setut (S) moraines are highlighted. These moraines are stratigraphically distinct from the sampled Arànser moraines but are currently undated. The margins of Arànser glacier can be traced further up valley but sampling was focused on the illustrated moraine area (light purple shading) in which the moraine margins are easily delineated (≤ 2 km from glacier terminus).

ice-proximal slope (IS) and the outer ice-distal slope (OS), while the number of selected boulders varied as a function of moraine size (n = 60 -275). In turn, boulder selection was primarily motivated by spatial location and the construction of a dense matrix of sampling points, rather than individual boulder characteristics. Each boulder was sampled using an N-type Schmidt hammer (SH) to assess the relative degree of weathering following the sampling approach of Tomkins et al. (2018a). All boulders were of sufficient size (> 25 kg; [START_REF] Sumner | The effect of rock moisture on Schmidt hammer rebound: tests on rock samples from Marion Island and South Africa[END_REF] and sampled areas were free of surface discontinuities [START_REF] Williams | The effect of surface texture on the determination of the surface hardness of rock using the Schmidt hammer[END_REF] and lichen [START_REF] Matthews | Endolithic lichens, rapid biological weathering and Schmidt hammer R-values on recently exposed rock surfaces: Storbreen Glacier Foreland, Jotunheimen, Norway[END_REF]. Thirty R-values were recorded for each boulder by a single operator and no outliers were removed following [START_REF] Niedzielski | A minimum sample size required from Schmidt hammer measurements[END_REF]. Schmidt hammer functioning was assessed regularly using the manufacturer's test anvil, with instrument and age calibration performed following Tomkins et al. (2018a). In total, 635 moraine boulders were sampled across the selected sites and 19,050 SH R-values were generated. To complement these data, the dimensions, surface features and depositional context of each sampled boulder were recorded (see Supplementary Information).

Schmidt hammer R-values correspond to the degree of surface weathering, assuming minimal lithological variation between tested rock surfaces [START_REF] Mccarroll | Potential and limitations of the Schmidt hammer for relativeage dating: field tests on neoglacial moraines, Jotunheimen, Southern Norway[END_REF], and are inversely proportional to surface exposure age. The rate and style of weathering and returned recalculated 10 Be ages of 6.3 ± 0.9 ka (IPN01) and 16.0 ± 2.5 ka (IPN02). Given the stratigraphic position of this deposit and limiting ages from the Outer Pleta Naua and Molières moraines (MUL01 = 14.9 ± 2.6 ka, MUL03 = 14.9 ± 1.9 ka; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF]), it appears likely that IPN02 is affected by inheritance.

may also be modified by climate [START_REF] Riebe | Erosional and climatic effects on longterm chemical weathering rates in granitic landscapes spanning diverse climate regimes[END_REF][START_REF] Portenga | Understanding Earth's eroding surface with 10 Be[END_REF][START_REF] Marrero | Controls on subaerial erosion rates in Antarctica[END_REF]. At the intra-landform scale, lithologic-climatic variability is absent as all sampled boulders share a common source area and climatic regime. At the interlandform scale, variability in rock type is minimal, as all sampled boulders were coarse-to medium-grained granites and granodiorites. Finally, while regional climatic variability could account for variability in weathering rates across the studied sites, previous work has shown that rates of sub-aerial weathering of granite are consistent over large spatial scales for regions of similar climate (Tomkins et al., 2018b).

Calculating SH-calibrated exposure ages

As granitic lithologies have proved effective for calibratedrelative age dating, SH R-values are used here as a proxy for exposure age based on a 10 Be-SH calibration dataset developed by Tomkins et al. (2018b). This dataset comprises 52 10 Be ages, distributed between 4.2 ± 0.3 ka and 51.8 ± 4.5 ka (Fig. 4), obtained from granite and granodiorite glacial boulders and glaciallysculpted bedrock from across the central and eastern Pyrenees and their corresponding SH R-values (Tomkins et al., 2018b). This dataset has been updated to include two additional 10 Be dated surfaces from the Val de Molières (MUL01 and MUL03; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF]; see Supplementary Table 1).

To utilise these data, 10 Be ages were recalculated using the CRONUS Earth Web Calculator (Version 2.0; [START_REF] Marrero | Cosmogenic nuclide systematics and the CRONUScalc program[END_REF], available at: http://cronus .cosmogenicnuclides .rocks /2 .0/, accessed: 01/09/2020), relative to the production rate dataset in [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF] and the time-dependent Lm scaling scheme [START_REF] Lal | Cosmic ray labeling of erosion surfaces: in situ nuclide production rates and erosion models[END_REF][START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF], and assuming 0 mm ka -1 erosion. Recalculated 10 Be ages are minimum estimates, as no corrections were made for shielding by snow, sediment or vegetation, surface erosion, or isostatic adjustment. To ensure consistency, all 10 Be and 36 Cl TCN ages discussed in this paper have been recalculated using these input parameters. This includes published ages from [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF], [START_REF] Rodés | La última deglaciación en los pirineos: de superficies de exposición mediante 10 Be, y modelado numérico de paleoglaciares[END_REF][START_REF] Palacios | Maximum extent of Late Pleistocene glaciers and last deglaciation of La Cerdanya mountains[END_REF], in addition to the 10 new 10 Be dated samples from the Arànser catchment. Full sample details used for exposure age calculation are provided in the Supplementary Information.

In turn, a 10 Be-SH calibration curve was constructed using logarithmic orthogonal distance regression (ODR, [START_REF] Boggs | Statistical Analysis of Measurement Error Models and Applications[END_REF] which minimises orthogonal residuals to account for measurement uncertainties in both the independent and dependent variables. We utilise Monte Carlo simulations to explicitly incorporate measurement errors; an approach which is preferable to a weighted ODR which requires unnecessary assumptions regarding weighting constants and is biased by TCN age-uncertainty collinearity [START_REF] Ivy-Ochs | Cosmogenic nuclides and the dating of lateglacial and early Holocene glacier variations: the Alpine perspective[END_REF][START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF]. Our analytical procedure, which returns prediction estimates (1σ ) of ± 2.0 -2.3 ka, is described fully in the Supplementary Information and has been implemented in SHED-Earth (http://shed .earth), an online calculator developed to enable wider and more consistent application of our approach (Tomkins et al., 2018a). To assess the accuracy of the 10 Be-SH calibration curve, 15 10 Be and 36 Cl ages from the studied moraines were located and re-sampled with the SH.

Based on this calibration curve, mean R-values from the 635 sampled boulders were converted into "SH-calibrated exposure ages" through interpolation. While uncertainty estimates for individual SH-calibrated exposure ages are larger than typical uncertainties associated with individual TCN exposure ages, landform age estimates can be of comparable precision to established techniques when derived from large SH datasets (e.g. n boulders ≥ 30; Tomkins et al., 2018b) and when appropriate statistical approaches for outlier identification and error propagation are employed [START_REF] Applegate | Improved moraine age interpretations through explicit matching of geomorphic process models to cosmogenic nuclide measurements from single landforms[END_REF][START_REF] Dortch | Timing and climatic drivers for glaciation across semi-arid western Himalayan-Tibetan orogen[END_REF][START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF].

Calculating landform ages

To determine the timing of moraine deposition at each site, we analysed the distribution of SH-calibrated exposure ages using the Probabilistic Cosmogenic Age Analysis Tool (P-CAAT Version 1.0; [START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF]. This method builds on the earlier work of [START_REF] Dortch | Timing and climatic drivers for glaciation across semi-arid western Himalayan-Tibetan orogen[END_REF] and utilises non-linear curve fitting and a Monte Carlo style approach to isolate component Gaussian distributions to account for positive (prior exposure) and negative skew (incomplete exposure) of age datasets. The results of this analysis are presented in Fig. 5 and Table 2. To assess the validity of these landform ages, we compared these data to the distribution of 10 Be and 36 Cl ages from the studied landforms (Table 1; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF][START_REF] Rodés | La última deglaciación en los pirineos: de superficies de exposición mediante 10 Be, y modelado numérico de paleoglaciares[END_REF][START_REF] Palacios | Maximum extent of Late Pleistocene glaciers and last deglaciation of La Cerdanya mountains[END_REF].

Based on landform age analysis, individual boulders were sorted into "good" and "bad" groups, which are defined by the 2σ (95%) age boundaries of the calculated landform age. Boulders which returned SH-calibrated exposure ages within 2σ of the landform age were classed as "good", while those younger or older than the landform age (> 2σ ) were classed as "bad". Selection of a broad 2σ threshold is appropriate given the measurement uncertainties associated with SH sampling, in addition to the systematic and geologic uncertainties inherited from TCN dating. Logistic analysis is used to distinguish boulders which correspond to the timing of moraine deposition or initial stabilisation ("good") from those which are likely compromised by pre-or post-depositional exposure ("bad").

Spatial analysis

The spatial distribution of "good" and "bad" boulders was analysed using global and local Moran's I spatial autocorrelation and based on a row-standardised distance-band weights matrix, where the distance band threshold is the minimum distance required to ensure that each boulder has at least two neighbours (Table 3). The Python implementation is available on GitHub: https://github .com / matt -tomkins /moraine -crest -or-slope. At the global level, Moran's I was used to assess whether the overall clustering of the data The summed probability density estimate (PDE) and lower probability component Gaussians are denoted by black and grey distributions, respectively. For each moraine, we include the bandwidth estimator used and its associated numeric bandwidth, the P-CAAT model fit (R 2 ), the total number of SH-calibrated exposure ages (n) and in brackets, the number of ages which are enclosed by the selected component Gaussian distribution at 2σ . Based on this approach, selected component Gaussians are interpreted to reflect the timing of moraine deposition or initial stabilisation. In contrast, younger component Gaussians may reflect post-depositional processes (e.g. moraine degradation, boulder exhumation or instability) while older component Gaussians likely incorporate pre-depositional processes (e.g. reworking of glacial deposits). was significantly different from a random distribution. For datasets that are non-random (p < 0.05), local Moran's I was used to identify the location of statistically significant boulder clusters (Fig. 6).

Current sampling approaches are based on the qualitatively-sound but quantitatively-untested assumptions that (i) the distribution of "good" boulders is non-random, and that (ii) "good" clusters are more likely on moraine crests. These assumptions can be explicitly tested for the studied moraines using global and local Moran's I respectively.

Sensitivity analysis

The above analyses provide important information on the relative occurrence and spatial clustering of "good" and "bad" boulders for moraines of varying age and morphology. However, this logistic classification is ultimately dependent on the calculated landform age, which will vary depending on the choice of numeric bandwidth estimator and the size and clustering of the input dataset [START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF]. To evaluate the reproducibility of our results, sensitivity testing was performed to evaluate the number of samples required to reproduce the estimated landform age based on 1σ and 2σ thresholds. The full analytical approach is described in the Supplementary Information and the results are presented in Fig. 7B.

Results

SH-calibrated exposure ages

There is a strong correlation between recalculated 10 Be ages and their corresponding SH R-values (Fig. 4; n = 54). Moreover, of the 15 10 Be and 36 Cl dated boulders re-sampled with the SH, the majority closely match the existing calibration dataset (n = 13). These observations indicate that when lithological variation is minimised, the relative degree of rock surface weathering can be used as a proxy for surface exposure age.

Exceptions to this correlation are samples ECH03 (17.2 ± 3.5 ka) and ECH04 (16.8 ± 3.3 ka) from the Soum d'Ech moraines [START_REF] Rodés | La última deglaciación en los pirineos: de superficies de exposición mediante 10 Be, y modelado numérico de paleoglaciares[END_REF] which are significantly more weathered (∼38 R) than their corresponding 10 Be ages would predict (∼47 R). This difference likely reflects sub-surface weathering prior to boulder exhumation. However, the scale of this influence is unlikely to be universal given the close correspondence between sample ECH01 (19.7 ± 3.6 ka) and the existing calibration dataset (see Fig. 4).

While sub-surface weathering of boulders under thin soil cover (∼25 cm) can occur [START_REF] Darmody | Chemical weathering and boulder mantles, Kärkevagge, Swedish Lapland[END_REF], boulders are often protected from weathering by sediment burial, as evidenced by the emergence of unweathered boulders from glacial tills and alluvium [START_REF] Ehlmann | Quantitative morphologic analysis of boulder shape and surface texture to infer environmental history: a case study of rock breakdown at the Ephrata Fan, Channeled Scabland[END_REF]. In turn, as SH-calibrated exposure ages from the Soum d'Ech moraines may well incorporate the effects of both sub-aerial and sub-surface weathering, and could also be influenced by weathering rate variability (e.g. differences between the Atlantic-(wet) and Mediterranean-influenced (dry) Pyrenees), it is possible that the estimated depositional age is an overestimate.

Landform ages

Landform ages derived from SH-calibrated exposure ages and associated P-CAAT model parameters are reported in Table 2 [START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF]. Based on this approach, latero-frontal moraines in the Arànser catchment were deposited at 23.3 ± 1.1 ka (left) and 22.3 ± 0.9 ka (right). As these estimates are consistent within measurement uncertainties, and given the comparable morpho-stratigraphy of these deposits (Fig. 2A), we consider moraine deposition to be contemporaneous. No independent dating evidence is available for the left lateral moraine, but 12 TCN ages are now available for the right lateral moraine ( 36 Cl, n = 2; 10 Be, n = 10). Using P-CAAT and selecting the oldest component Gaussian distribution that contains ≥ 3 ages to represent the age of the landform (see Fig. 3 in [START_REF] Dortch | Timing and climatic drivers for glaciation across semi-arid western Himalayan-Tibetan orogen[END_REF], these data return a landform age of 21.5 ± 2.2 ka (Mean bandwidth estimator; Numeric bandwidth = 0.8108, R 2 = 0.9997, p < 0.01), while the oldest sample is 22.4 ± 1.8 ka (SAL-10). Both estimates are consistent within measurement uncertainties with the SH-derived landform ages.

In the Gave de Pau catchment, SH-calibrated exposure ages from the proximal Soum d'Ech lateral moraines return landform ages of 26.2 ± 2.5 ka (outer, n = 61) and 26.1 ± 1.7 ka (inner, n = 39). While these moraines are morpho-stratigraphically distinct, they cannot be statistically distinguished. It is possible that moraine deposition occurred within the resolution of our sampling approach, or that differences in moraine age have been masked by moraine stabilisation, degradation or sub-surface boulder weathering. As the temporal distribution of SH-calibrated exposure ages is near identical (Table 2), we assign these deposits a landform age of 27.3 ± 1.8 ka based on P-CAAT (n = 100; STD / IQR bandwidth estimator; Numeric bandwidth = 0.9877, R 2 = 0.9989, p < 0.01), and perform subsequent analyses on the combined dataset for computational ease. While this estimate is significantly older than the corresponding 10 Be ages (16.8 ± 3.3 ka, 17.2 ± 3.5 ka, 19.7 ± 3.6 ka; [START_REF] Rodés | La última deglaciación en los pirineos: de superficies de exposición mediante 10 Be, y modelado numérico de paleoglaciares[END_REF], it is consistent with limiting 14 C ages obtained from a proximal palaeolake sediment sequence at Lac de Lourdes [START_REF] Reille | The late Pleistocene and Holocene in the Lourdes Basin, Western Pyrénées, France: new pollen analytical and chronological data[END_REF]. While the oldest radiometric 14 C ages from this over-deepened glacial basin are now considered suspect due to contamination from mineral carbon [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF] likely that the ECH samples are representative of final moraine stabilisation, rather than initial deposition. This interpretation is supported by sensitivity analysis (see Section 3.5), as the number of TCN ages (n = 3) is below the threshold required to consistently reproduce the landform age derived from the full dataset at 2σ (n = 11) and 1σ (n = 41).

In the Val de Molières catchment of the Noguera Rigaborçana, recalculated 10 Be ages on the Outer Pleta Naua moraine range from 12.6 ± 1.5 ka to 13.2 ± 1.6 ka (n = 3; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF]. These estimates are consistent with the SH-calibrated exposure ages, which range from 11.8 ± 2.0 ka to 13.1 ± 2.0 ka (n = 60). As the SHcalibrated exposure ages conform to a normal distribution (Fig. 5B; Shapiro-Wilk test, W = 0.96, p = 0.07), are well-clustered (IQR = 0.6 ka), and return an excellent P-CAAT model fit with a single com- In turn, the Outer Pleta Naua moraine was likely deposited at 12.5 ± 0.4 ka. Applying the same analytical approach (x ± t) to the corresponding 10 Be ages produces 12.9 ± 1.0 ka, which is statistically indistinguishable. Moreover, these estimates are stratigraphically consistent with independent landform ages in the Val de Molières catchment [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF], with maximum and minimum limiting ages for moraine deposition provided by samples from the Molières (MUL01 = 14.9 ± 2.6 ka, MUL03 = 14.9 ± 1.9 ka) and Inner Pleta Naua moraines respectively (Fig. 3B; IPN01 = 6.3 ± 0.9 ka; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF].

Finally, the Tallada cirque moraine returned a landform age of 3.2 ± 0.7 ka. While this SH-derived estimate cannot be independently verified, the limited weathering of the moraine boulders (SH R ≥ 60), in combination with the topographic setting of the Tallada cirque, appears consistent with a late-Holocene origin.

Temporal distribution

Estimated landform ages are generally consistent with independent TCN ages (n = 19) but the age distribution of SH-calibrated exposure ages varies significantly between the sampled moraines (Fig. 5). For the Arànser and Soum d'Ech moraines, the distribution of SH-calibrated exposure ages is strongly negatively skewed (Table 2), in line with exhumation models [START_REF] Applegate | Improved moraine age interpretations through explicit matching of geomorphic process models to cosmogenic nuclide measurements from single landforms[END_REF], while Tallada is normally distributed with a slight positive skew (Shapiro-Wilk test, W = 0.98, p = 0.56); a trend which may reflect prior exposure or reworking of glacial material [START_REF] Applegate | Modeling the statistical distributions of cosmogenic exposure dates from moraines[END_REF].

In light of these trends, the proportion of "good" and "bad" boulders, as defined by the 2σ age boundaries of the corresponding landform age, varies between the sampled moraines. The proportion of "good" boulders is highest on the Outer Pleta Naua moraine (100%) and lowest on the Arànser left (56%) and Arànser right moraines (49%). For moraines corresponding to the LGM, most "bad" boulders are younger than the assumed age of deglaciation (Table 2), while the Holocene Tallada moraine contains a small but significant component of boulders which are older than the assumed age of deglaciation (14%). Logistic analysis indicates that boulder characteristics (e.g. boulder height) did not have a consistent statistically significant effect on the distribution of "good" and "bad" boulders across the sampled moraines (see Supplementary Information).

Spatial distribution

Summary statistics for spatial analysis are presented in Table 3. This approach reveals marked inter-landform variation, with statistically significant spatial clustering absent from the Tallada, Outer Pleta Naua and Arànser right moraines (simulated p > 0.05). In turn, the spatial distribution of "good" and "bad" boulders for these moraines is effectively random.

One exception to this rule is the Arànser left moraine where statistically significant clustering is evident (simulated p < 0.05) and where clusters identified using local Moran's I have plausible geomorphological explanations (Fig. 6A). Clusters of "young" boulders occur:

(i) at the moraine terminus, (ii) where the moraine crest has been cross-cut and incised by a minor stream and, (iii) where boulders have accumulated at the base of the moraine slope.

Additional clusters are also evident on the ice-proximal slope (Fig. 6A). Clusters (i) and (ii) are likely fluvial in origin, with the former explained by incision of the terminal deposits, which may have led to degradation of the lateral flanks and exhumation of moraine boulders. This pattern of post-depositional degradation matches the spatial clustering of 36 Cl ages on a comparable LGM moraine deposited in the nearby Duran valley (see Fig. 11 in Pala-cios et al., 2015). The second cluster may be partially explained by meltwater erosion, given the proximity of the incised area to the former terminus of the Setut glacier (Fig. 2A). The origins of the remaining "young" clusters are less clear, but these ultimately reflect instability of the ice-proximal slope, although it is not yet clear whether this was driven by autogenic moraine stabilisation or external factors (e.g. subsequent glacial advance, fluvial erosion). Clusters of "good" boulders were also identified on the Arànser left moraine but these are distributed across moraine crests and ice-proximal and -distal slopes and follow no clear spatial pattern. Finally, local Moran's I identified both "young" and "good" clusters on the outer Soum d'Ech moraine (Fig. 6B) but there is no clear geomorphological evidence which explains their distribution.

The proportion of "good" boulders varies markedly between the studied moraines, but this overall trend is relatively consistent across boulder groups (C, IS, OS) at the intra-landform scale (Fig. 7A). While there are clear differences between boulder groups at the Arànser left and right moraines, there are no consistent trends at the inter-landform scale and no single boulder group performs optimally across all landforms.

Sensitivity results

Based on the sensitivity analysis described in Section 2.6, there are clear differences in the number of SH samples required to reproduce the landform ages obtained from the full datasets (n = 60 -275; Fig. 7B). The Outer Pleta Naua landform age requires only three samples at both 1σ and 2σ . Landform ages for both the Arànser left and right moraines can be reproduced with relatively few samples at both 1σ (n ≤ 26) and 2σ (n ≤ 16), while both the Soum d'Ech and Tallada moraines require ≥ 40 samples to reproduce the landform age at 1σ . These trends are largely explained by the degree of overlap between component Gaussian distributions (see Fig. 5). Both the Tallada (Fig. 5A) and Soum d'Ech moraines (Fig. 5E) feature lower probability component Gaussians, centred on 4.7 ± 0.9 ka and 24.4 ± 1.7 ka respectively, which overlap with the highest probability component Gaussian. In contrast, there is minimal overlap between component Gaussians for the Arànser left moraine (Fig. 5C), despite the high degree of dataset skew and the large number of "bad" boulders (44%). The Arànser right moraine is intermediate in character (Fig. 5D), with clear unidirectional skew but a greater degree of overlap between the highest probability Gaussian (22.3 ± 0.9 ka) and younger lower probability component Gaussians (17.6 ± 2.9 ka; 20.9 ± 0.9 ka). This distribution explains the larger number of samples required at both 1σ and 2σ relative to the Arànser left moraine. Ultimately, as the degree of overlap between component Gaussians increases, more samples are required to isolate the highest probability component Gaussian and eliminate PDE skew. Despite this, all landform ages could be reproduced with relatively few samples at both 1σ (n ≤ 40) and 2σ (n ≤ 26). While these values exceed typical sample size recommendations for TCN dating [START_REF] Putkonen | Accuracy of cosmogenic ages for moraines[END_REF], they are based upon strict thresholds (≥ 95% of simulated landform ages within 1σ or 2σ of the full dataset landform age) and should be utilised by researchers when pre-screening a larger population of boulders prior to targeted TCN sampling.

Discussion

Efforts to minimise sampling bias of moraine TCN datasets may significantly improve the utility of moraine chronologies in determining glacial history and the climatic drivers of glacial cycles. However, while careful geomorphological assessment of individual boulders is necessary to isolate those influenced by pre-or post-depositional processes, many criteria for TCN sample selection have not been tested quantitatively. Of these, boulder location is traditionally thought to be critical, as moraine crest boulders are prioritised due to perceived stability (e.g. [START_REF] Gosse | Precise cosmogenic 10 Be measurements in western North America: support for a global Younger Dryas cooling event[END_REF][START_REF] Hallet | Surface dating of dynamic landforms: young boulders on aging moraines[END_REF], while those deposited on ice-proximal or -distal slopes are typically rejected. This study is the first to quantitatively assess this approach.

Based on 10 Be (n = 10) and Schmidt hammer sampling (n = 635) of ice-marginal moraines in the Pyrenees, it is clear that the spatial distribution of SH-calibrated exposure ages is both complex and site-specific. For many moraines, the distribution of "good" and "bad" boulders is effectively random (p > 0.05), as assessed using global Moran's I (Table 3), while in others, clusters of "good" and "bad" boulders have clear geomorphological explanations. More fundamentally, the likelihood of selecting a "good" boulder is comparable for moraine crests, ice-proximal and -distal slopes (Fig. 7A). Although statistically significant spatial clustering is evident for the Arànser left and Soum d'Ech moraines (p < 0.05; Fig. 6), the distribution of "good" boulder clusters is complex, with clusters distributed across moraine crests and moraine slopes.

While there are no consistent spatial patterns at the interlandform scale, the temporal distribution of SH-calibrated exposure ages varies markedly between the studied landforms, with a number of important observations. First, moraine sedimentology appears to place a key control on post-depositional stability [START_REF] Zreda | Cosmogenic 36 Cl accumulation in unstable landforms: 2. Simulations and measurements on eroding moraines[END_REF][START_REF] Putkonen | Degradation of unconsolidated quaternary landforms in the western North America[END_REF], as age distributions for matrix-rich moraines (e.g. Arànser, Soum d'Ech) are strongly negatively skewed (Fig. 6), with many boulders younger than the assigned age of the landform. The frequency of "young" boulders for these moraines (Table 2) likely reflects the influence of diffusive slope processes [START_REF] Applegate | Modeling the statistical distributions of cosmogenic exposure dates from moraines[END_REF], as the transfer of sediment to the base of moraine slopes drives exhumation of entrained boulders [START_REF] Porter | 36 Cl dating of the classic Pleistocene glacial record in the northeastern Cascade Range, Washington[END_REF] and erosion of moraine crests [START_REF] Schaller | Quantifying glacial moraine age, denudation, and soil mixing with cosmogenic nuclide depth profiles[END_REF] and leads to increasingly subdued moraine topography [START_REF] Putkonen | Degradation of unconsolidated quaternary landforms in the western North America[END_REF]. The clearest signal of moraine degradation is evident at the Arànser left (IQR = 7.9 ka; Skew = -1.02) and Arànser right moraines (IQR = 6.9 ka; Skew = -1.13) and this trend may be partially explained by forest growth and boulder toppling [START_REF] Ivy-Ochs | Cosmogenic nuclides and the dating of lateglacial and early Holocene glacier variations: the Alpine perspective[END_REF], as well as the effects of fluvial incision (Fig. 6), while historic land use may also play a role [START_REF] Pallàs | Small, isolated glacial catchments as priority targets for cosmogenic surface exposure dating of Pleistocene climate fluctuations, southeastern Pyrenees[END_REF].

In contrast, the boulder-rich, matrix-poor Outer Pleta Naua moraine stabilised rapidly after glacial retreat, as evidenced by the distribution and clustering of both its SH-calibrated exposure ages (IQR = 0.6 ka; Shapiro Wilk W = 0.96, p = 0.07) and the corresponding 10 Be dataset [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF]. The sedimentology of the Outer Pleta Naua moraine is likely a function of catchment size and glacier area, and the short transport distance from the bedrock source area (Fig. 3; ≤ 300 m). In the absence of a supporting sediment matrix, boulder-rich moraines stabilise quickly and appear less susceptible to subsequent erosion [START_REF] Ivy-Ochs | Cosmogenic nuclides and the dating of lateglacial and early Holocene glacier variations: the Alpine perspective[END_REF][START_REF] Pallàs | Small, isolated glacial catchments as priority targets for cosmogenic surface exposure dating of Pleistocene climate fluctuations, southeastern Pyrenees[END_REF]. Finally, for moraines deposited by niche cirque glaciers, reworking of glacial, periglacial or rockfall material appears more significant than post-depositional modification, in line with previous studies [START_REF] Heyman | Too young or too old: evaluating cosmogenic exposure dating based on an analysis of compiled boulder exposure ages[END_REF]. In these environments, the age of the oldest boulder may overestimate the "true" age of the moraine [START_REF] Putkonen | Accuracy of cosmogenic ages for moraines[END_REF][START_REF] Briner | Cosmogenic exposure dating of late Pleistocene moraine stabilization in Alaska[END_REF].

Implications for TCN sampling of moraines

The results described above have implications for future sampling approaches. First, while "good" boulders are not more likely on moraine crests, we find there is no clear penalty to moraine crest sampling, as initial differences between moraine crests and ice-proximal and -distal slopes appear to be masked by continued moraine degradation. Thus, in the absence of detailed geomorphological assessments of individual landforms, restricting sampling to moraine crests is a viable strategy to minimise the likelihood of boulder instability, assuming there are sufficient numbers of boulders to select from. This finding is unlikely to hold true for recently deposited (< 1 ka) unconsolidated landforms [START_REF] Putkonen | Degradation of unconsolidated quaternary landforms in the western North America[END_REF], whose over-steepened ice-proximal slopes have yet to stabilise [START_REF] Briner | Cosmogenic exposure dating of late Pleistocene moraine stabilization in Alaska[END_REF][START_REF] Dortch | Beryllium-10 surface exposure dating of glacial successions in the Central Alaska Range[END_REF].

Second, our results show that sampling boulders on iceproximal and -distal slopes can be as effective as sampling moraine crests (Fig. 7A). While boulder density is typically highest at moraine crests [START_REF] Putkonen | Landscape evolution degrades the geologic signature of past glaciations[END_REF], there is no guarantee that these boulders are the best options for TCN dating. Moreover, if sample selection criteria are rigorously applied, the number of suitable boulders available for dating could fall below a critical level. Without robust statistical identification of outliers, this could lead to unclear results given the ubiquity of post-depositional modification of moraines [START_REF] Zech | Evidence for long-lasting landform surface instability on hummocky moraines in the Pamir Mountains (Tajikistan) from 10 Be surface exposure dating[END_REF][START_REF] Heyman | Too young or too old: evaluating cosmogenic exposure dating based on an analysis of compiled boulder exposure ages[END_REF]. One strategy which is rarely utilised is to select boulders for TCN dating from ice-proximal and -distal slopes, but evidence from the studied moraines indicates that this is a viable strategy, as the proportion of "good" boulders is comparable to moraine crests (Fig. 7A). For many moraines, the spatial distribution of "good" boulders is random, while statistically significant clusters of "good" boulders are distributed across moraine crests and moraine slopes (Fig. 6). These observations indicate that redefining selection criteria to include the entire population of moraine boulders would have no clear negative effect and could prove beneficial for moraines where ideal boulders are rare or are distributed away from moraine crests.

Third, our data indicate that landform characteristics have a clear impact on the temporal distribution of SH-calibrated exposure ages (Fig. 5; [START_REF] Putkonen | Degradation of unconsolidated quaternary landforms in the western North America[END_REF][START_REF] Ivy-Ochs | Cosmogenic nuclides and the dating of lateglacial and early Holocene glacier variations: the Alpine perspective[END_REF][START_REF] Pallàs | Small, isolated glacial catchments as priority targets for cosmogenic surface exposure dating of Pleistocene climate fluctuations, southeastern Pyrenees[END_REF]. Within this context, we suggest that landform stability should be prioritised, as differences between landforms appear far greater than differences between boulder groups on an individual landform (C vs. IS vs. OS). Differences are evident as a function of moraine sedimentology [START_REF] Zreda | Cosmogenic 36 Cl accumulation in unstable landforms: 2. Simulations and measurements on eroding moraines[END_REF], with rapid stabilisation of matrix-poor, boulder-rich moraines (e.g. Outer Pleta Naua; [START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF][START_REF] Pallàs | Small, isolated glacial catchments as priority targets for cosmogenic surface exposure dating of Pleistocene climate fluctuations, southeastern Pyrenees[END_REF][START_REF] Ivy-Ochs | Cosmogenic nuclides and the dating of lateglacial and early Holocene glacier variations: the Alpine perspective[END_REF] but prolonged degradation of unconsolidated landforms (e.g. Arànser; [START_REF] Putkonen | Degradation of unconsolidated quaternary landforms in the western North America[END_REF][START_REF] Dortch | Beryllium-10 surface exposure dating of glacial successions in the Central Alaska Range[END_REF]. Although moraine sedimentology has explanatory power for the studied moraines, the observed trends are unlikely to hold true in all settings due to climatic and topographic controls on moraine stability [START_REF] Barr | A review of topographic controls on moraine distribution[END_REF]. Moreover, restricting sampling to matrix-poor landforms could have unintended adverse effects, as moraines may incorporate supraglacial rock avalanche debris and may primarily preserve a non-climatic signal. Alternatively, sampling unconsolidated landforms does not guarantee poor clustering (e.g. χ 2 > 1), particularly in regions where moraine denudation is limited by climate [START_REF] Zech | Evidence for long-lasting landform surface instability on hummocky moraines in the Pamir Mountains (Tajikistan) from 10 Be surface exposure dating[END_REF][START_REF] Morgan | Degradation of glacial deposits quantified with cosmogenic nuclides, Quartermain Mountains, Antarctica[END_REF][START_REF] Balter | A 14.5-million-year record of East Antarctic Ice Sheet fluctuations from the central Transantarctic Mountains, constrained with cosmogenic 3 He, 10 Be, 21 Ne, and 26 Al[END_REF] or where topographic factors promote moraine stability [START_REF] Barr | A review of topographic controls on moraine distribution[END_REF]. Finally, restricting sampling to landforms with specific characteristics is often not viable, as key glacial chronological markers may be represented by only a small number of landforms.

Within this context, we suggest that landform selection is critical, and care should be taken to select methods which are appropriate for its assumed age and stability and to collect a sufficient number of samples to enable robust outlier identification [START_REF] Putkonen | Accuracy of cosmogenic ages for moraines[END_REF]. However, it is often challenging to assess landform stability based on geomorphological evidence alone. Our approach, in light of strong regional evidence for an inverse correlation between SH R-values and exposure ages for granitic surfaces (Tomkins et al., 2018a(Tomkins et al., , 2018b)), indicates that preliminary SH sampling could be a useful method to assess landform stability, to identify boulders affected by post-depositional processes, and to prioritise individual boulders for analysis based on R-value clustering [START_REF] Tylmann | Erratics selection for cosmogenic nuclide exposure dating -an optimization approach[END_REF].

Based on the sensitivity approach described in Section 2.6, the number of SH samples required scales with the complexity of the underlying distribution (Fig. 7B), from those which are approximately normal to those which feature overlapping component Gaussian distributions (Fig. 5) or multi-directional skew (i.e. preand post-depositional skew). However, given that it is not possible to ascertain the underlying distribution a priori, a relatively large sample size is ultimately required. For most landforms, sampling a minimum of ∼30 boulders would be a reasonable approach to estimate a depositional age within 2σ (n ≥ 23), but more would be required (n ≥ 40) to improve precision to 1σ for complex datasets or if Schmidt hammer R-values were being used as a basis for cosmogenic nuclide sample selection [START_REF] Tylmann | Erratics selection for cosmogenic nuclide exposure dating -an optimization approach[END_REF]. Collecting a minimum of 30 -40 samples is necessary to ensure a full understanding of the underlying age distribution, even for complex datasets. Based on this preliminary sampling, statistical approaches could be used to isolate component Gaussian distributions [START_REF] Dortch | Timing and climatic drivers for glaciation across semi-arid western Himalayan-Tibetan orogen[END_REF][START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF] and to identify individual boulders which are consistent with the age of the landform and to reject those which are "young" or "old" [START_REF] Heyman | Too young or too old: evaluating cosmogenic exposure dating based on an analysis of compiled boulder exposure ages[END_REF].

Finally, it is important to note that the effectiveness of this approach may vary as a function of lithology and climate (Mc-Carroll, 1989), while the underlying measurements are sensitive to factors which have only a minor effect on cosmogenic nuclide concentrations (e.g. surface discontinuities, [START_REF] Williams | The effect of surface texture on the determination of the surface hardness of rock using the Schmidt hammer[END_REF]lichen coverage, Matthews and Owen, 2008). However, when these limitations are accounted for, Schmidt hammer R-values can be used as a proxy for surface exposure age (Fig. 4). Given the ubiquity of geologic scatter (e.g. exhumation, erosion, shielding), incorporating time-and cost-efficient preliminary SH sampling as an additional tool for TCN sample selection could ultimately improve the chronological utility of the moraine record and enable a deeper understanding of the climatic drivers of glacial cycles.

Conclusions

Based on 10 Be and Schmidt hammer sampling of ice-marginal moraines in the Pyrenees, this study provided the first quantitative analysis of the relative utility of moraine crest and moraine slope sampling for terrestrial cosmogenic nuclide dating. Using spatial analysis of SH-calibrated exposure ages, we show that there is no clear penalty to moraine slope sampling. However, contrary to current sampling approaches, which typically prioritise moraine crest boulders due to perceived stability, we show that the proportion of "good" boulders is comparable between moraine crests and ice-proximal and -distal slopes, while for many moraines, the spatial distribution of "good" boulders is effectively random. Crucially, however, differences between landforms appear more significant than differences at the intra-landform scale; a result which indicates that the stability of the landform can have a far greater impact on the distribution of boulder exposure ages than the characteristics and depositional context of individual boulders. In this study, moraine sedimentology likely accounts for the observed differences between landforms, with rapid stabilisation of matrix-poor, boulder-rich moraines and prolonged degradation for unconsolidated landforms. Although these trends are unlikely to be universally applicable given climatic and topographic controls on moraine stability, our data indicate that preliminary SH sampling is a valuable tool to assess landform stability and to prioritise individual boulders for further analysis.

Fig. 1 .

 1 Fig. 1. Site photographs of the (A) Tallada, (B) Outer Pleta Naua, (C) Arànser and (D) Soum d'Ech moraines (denoted by red arrows). (E-F) Topographic maps of the Pyrenees (ASTER GDEM V3, WGS 84 UTM 31N), showing the locations of the studied catchments and selected moraines and the distribution of Axial Zone granites within those catchments. The latter was derived from a 1:400,000 geological map produced by the IGME (Spain) and the BRGM (France). Also shown are the locations of major summits ( Aneto, C arlit, Estats, Plana de Lles, M onte P erdido, P osets, V ignemale) and the extent of glaciers during the global Last Glacial Maximum (LGM; Calvet et al., 2011). (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Fig. 2 .

 2 Fig. 2. Geomorphological maps for the (A) Arànser and (B) Soum d'Ech moraines (WGS 84 UTM 31N). These moraines likely correspond to the maximum ice extent (MIE) during the Würmian glacial stage (11.7 -110 ka; Calvet et al., 2011). Locations and sample names for TCN dated boulders are shown (white circles; Rodés, 2008; Palacioset al., 2015). In (A), the locations of the proximal Fornell (F) and Setut (S) moraines are highlighted. These moraines are stratigraphically distinct from the sampled Arànser moraines but are currently undated. The margins of Arànser glacier can be traced further up valley but sampling was focused on the illustrated moraine area (light purple shading) in which the moraine margins are easily delineated (≤ 2 km from glacier terminus).

Fig. 3 .

 3 Fig. 3. Geomorphological maps for the (A) Tallada and (B) Outer Pleta Naua moraines in the Val de Molières catchment of the Noguera Rigaborçana (WGS 84 UTM 31N).Locations and sample names for TCN dated boulders are shown (white circles;[START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF]. The Inner Pleta Naua moraine was also investigated by[START_REF] Pallàs | Late Pleistocene and Holocene glaciation in the Pyrenees: a critical review and new evidence from 10 Be exposure ages, South-central Pyrenees[END_REF]

Fig. 4 .

 4 Fig. 4. (A) Location of exposure age calibration sites (blue points) in the Bassies (B, n = 6), Carlit (C, n = 3), Noguera Rigaborçana (N, n = 4), Maladeta (M a , n = 9), Malniu (M n , n = 21), Molières (M o , n = 2), Orri (O, n = 3) and Querol catchments (Q, n = 6). Underlying topography is ASTER GDEM V3 (WGS 84 UTM 31N). Also shown are the locations of sampled moraines (orange points; see Fig. 1F) and the maximum ice extent (MIE) during the global Last Glacial Maximum (LGM; Calvet et al., 2011). (B) Monte Carlo-derived orthogonal distance regression (ODR) between 54 10 Be exposure ages (blue points ± external age uncertainty) and their corresponding SH R-values (mean of 30 R-values ± Standard Error of the Mean; Tomkins et al., 2018b), plus 1σ (blue dashed lines) and 2σ prediction limits (grey dashed lines). Independent TCN samples ( 10 Be, 36 Cl) from the studied moraines (n = 15) are shown as orange points. Inherited outliers from the original calibration dataset (n = 2; Tomkins et al., 2018b) are not shown for clarity. (C) Example of a 10 Be dated boulder from the Arànser right moraine (SAL-10).

Fig. 5 .

 5 Fig. 5. Gaussian decomposition of SH-calibrated boulder exposure ages for the Tallada (A), Outer Pleta Naua (B), Arànser (C-D) and Soum d'Ech moraines (E). Following P-CAAT guidelines[START_REF] Dortch | Timing and climatic drivers for glaciation across semi-arid western Himalayan-Tibetan orogen[END_REF][START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF], we selected the highest probability component Gaussian (red shading) to represent the age of the landform as all are ≤ LGM.

Fig. 6 .

 6 Fig. 6. Results of local Moran's I spatial autocorrelation for the Arànser left (A) and Soum d'Ech moraines (B).Points denote the location of sampled boulders, with neighbouring boulders linked by grey lines. Neighbours are calculated based on a fixed distance, defined as the minimum distance required to ensure that each boulder has at least two neighbours, and were analysed using inverse distance weighting (IDW). Points are coloured based on the results of local Moran's I , with regions of no statistically significant spatial clustering shown as white, while clusters of "good" (HH) and "bad" boulders (LL) and their contributing neighbours are shown in blue and red, respectively. Outlier points (HL and LH) are not shown for clarity. A histogram illustrating the distribution of calibrated boulder exposure ages is included for each moraine, coloured by the "good" (blue) and "bad" components (grey).ponent Gaussian (R 2 = 1, p < 0.01), we use the arithmetic mean (x) to represent the age of the landform and estimate the total uncertainty (t) following[START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF] as follows:

Fig. 7 .

 7 Fig. 7. The likelihood of sampling a "good" boulder (%; within 2σ of the landform age) for each of the studied moraines (A), subset by boulder position (inner ice-proximal slope, moraine crest, outer ice-distal slope). Sensitivity results are shown for each moraine (B), illustrating the number of samples required to reproduce the associated landform age within 1σ and 2σ thresholds.

Table 1

 1 Summary data for terrestrial cosmogenic nuclide exposure ages from the sampled moraines a .

	Moraine	Name	Isotope	Latitude	Longitude	Elevation	Age	Internal	External	SH R ± SEM b
				( • )	( • )	(m)	(ka)	± (ka)	± (ka)
	Outer Pleta Naua c	OPN01	10 Be	42.6365	0.7399	2217	13.2	1.3	1.6	-
		OPN02	10 Be	42.6365	0.7406	2197	13.0	1.7	2.0	51.68 ± 0.5
		OPN03	10 Be	42.6365	0.7409	2195	12.6	1.2	1.5	-
	Arànser (Right) d	SAL-01	10 Be	42.4283	1.6300	2000	17.6	0.6	1.5	47.57 ± 0.83
		SAL-02	10 Be	42.4273	1.6321	1983	19.2	0.6	1.5	45.07 ± 0.84
		SAL-03	10 Be	42.4270	1.6326	1975	21.1	0.6	1.7	44.07 ± 0.82
		SAL-04	10 Be	42.4254	1.6358	1933	18.0	0.6	1.5	47.57 ± 0.84
		SAL-05	10 Be	42.4240	1.6389	1912	17.0	0.9	1.6	48.9 ± 0.77
		SAL-06	10 Be	42.4237	1.6395	1908	19.2	0.6	1.5	44.53 ± 0.74
		SAL-07	10 Be	42.4229	1.6415	1896	16.7	0.5	1.4	47.43 ± 0.96
		SAL-08	10 Be	42.4223	1.6447	1863	17.1	0.6	1.4	47.7 ± 0.9
		SAL-09	10 Be	42.4215	1.6481	1820	20.7	0.9	1.7	44.77 ± 0.8
		SAL-10	10 Be	42.4213	1.6489	1808	22.4	0.7	1.8	43.03 ± 0.95
		PIR-11-13	36 Cl	42.4213	1.6495	1809	18.2	1.6	2.1	-
		PIR-11-14	36 Cl	42.4209	1.6499	1805	17.3	1.7	2.2	47.6 ± 0.83
	Soum d'Ech e	ECH01	10 Be	43.0863	-0.0870	776	19.7	3.2	3.6	42.43 ± 0.98
		ECH02	10 Be	43.0858	-0.0880	778	59.0	43.2 f	43.0	-
		ECH03	10 Be	43.0862	-0.0873	779	17.2	3.3	3.5	38.86 ± 1.11
		ECH04	10 Be	43.0865	-0.0867	781	16.8	3.0	3.3	38.77 ± 1.05

a Full sample information used for exposure age calculation is provided in the Supplementary Information or is available on GitHub: https://github .com /matt -tomkins /moraine -crest -or-slope. b Mean of 30 SH R-values ± the Standard Error of the Mean. c OPN samples from Pallàs et al. (2006). d PIR samples from Palacios et al. (2015). e ECH samples from Rodés (2008). f Measurement error, see Rodés (2008).

Table 2

 2 Age statistics for the sampled moraines. Method used for kernel density estimation (see[START_REF] Dortch | Probabilistic Cosmogenic Age Analysis Tool (P-CAAT), a tool for the ages[END_REF] and its associated numeric bandwidth. c All model p values < 0.01. d Reported uncertainty (±) is the 1σ bounds (68%) of the highest probability component Gaussian, unless stated otherwise. e Interquartile range. f Shapiro-Wilk test for normality p values. g Based on the landform age ± 2σ . h Arithmetic mean of 60 samples ± total uncertainty. i Calculation based on a reduced dataset of 274 samples. Sample ARL-192 (1.97 ± 2.06 ka) is more than three standard deviations from the mean of the remaining samples and was removed for program stability.

	Moraine	Group	Method a	Bandwidth b	Model fit c	Age	IQR e	Skew	Normality f	Young	Good	Old
						(ka) d				(%) g	(%) g	(%) g
	Tallada	-	STD / IQR	0.3731	0.9985	3.2 ± 0.7	1.2 ka	0.34	0.44	6	80	14
	Outer Pleta Naua	-	Mean	2.016	1	12.5 ± 0.4 h	0.6 ka	-0.24	0.07	0	100	0
	Arànser	Left	MAD	0.7003	0.9978	23.3 ± 1.1 i	7.9 ka	-1.02	< 0.01	44	56	0
		Right	MAD	0.6796	0.9991	22.3 ± 0.9	6.9 ka	-1.13	< 0.01	51	49	0
	Soum d'Ech	Outer	STD / IQR	1.0734	0.998	26.2 ± 2.5	3.5 ka	-1.49	< 0.01	-	-	-
		Inner	STD / IQR	1.1661	0.9996	26.1 ± 1.7	3.5 ka	-1.05	< 0.01	-	-	-
		Combined	STD / IQR	0.9877	0.9989	27.3 ± 1.8	3.6 ka	-1.49	< 0.01	24	76	0
	a,b

Table 3

 3 Spatial statistics for the sampled moraines. Inner ice-proximal slope (IS), moraine crest (C) and outer ice-distal slope (OS). b Defined as the minimum distance required to ensure that each boulder has at least two neighbours. c p values > 0.05 support no statistically significant spatial clustering. p values ≤ 0.05 are consistent with a non-random distribution and spatial clustering of the input data. d Spatial autocorrelation was not possible for the Outer Pleta Naua moraine as all boulders were classed as "good" based on the 2σ threshold.

			Number of samples		Global Morans I			"Good" boulder (%)
	Moraine	Type	Total	IS a	C a	OS a	Distance threshold (m) b	I	Simulated p value c	IS a	C a	OS a
	Tallada	Terminal	70	16	29	25	21.6	0.0980	0.0719	80	79	81
	Outer Pleta Naua	Terminal	60	20	20	20	23.9	NA d	NA d	100	100	100
	Arànser (Left)	Latero-frontal	275	199	51	25	59.5	0.0915	0.0064	53	57	76
	Arànser (Right)	Latero-frontal	130	57	33	40	66.3	0.0651	0.1194	63	36	40
	Soum d'Ech	Laterals	100	37	50	13	51.1	0.1519	0.0106	76	72	81

a

Acknowledgements

MT was funded by a University of Manchester President's Doctoral Scholar Award. The fieldwork for this paper was supported by a British Society for Geomorphology Postgraduate Research grant awarded to MT. 10 Be analysis was performed at the ASTER AMS national facility (CEREGE), which is supported by the CNRS/INSU, by the ANR, through the "projets thématiques d'excellence" programme for the "Equipements d'excellence" ASTER-CEREGE action, and by the Institut de Recherche pour le Développement (IRD). MT, JD, PH and RP designed the project, with refinements from IB and CD. MT, JA and AS conducted Schmidt hammer sampling. 10 Be sampling and analysis was undertaken by RP, AR, DB, VR, VJ, LRR and RC. MT and JH conducted spatial analysis. MT, JH and TB improved error propagation. All authors contributed to the manuscript. We would like to thank Professor W. Amidon for his constructive review. We thank Dr. G. Evatt and Dr. A. Smedley at the University of Manchester for their comments, D. Tomkins and Dr. B. van Bodegraven for fieldwork support, and the contributors to the open-source libraries PySAL and GeoPandas.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary material

Supplementary material related to this article can be found online at https://doi .org /10 .1016 /j .epsl .2021.117092.