Accurate estimates of past topography are required to reliably reconstruct past ice sheets to infer palaeoclimate. For this reason, understanding erosion rates across East Greenland is crucial to constrain landscape evolution driven by tectonics and climate-dependent erosion rates. Here we analyse published apatite fission-track (AFT) data to constrain the spatial pattern of AFT bedrock ages across the landscape. We compare these bedrock ages with published detrital distribution to highlight ambiguity in the pattern of erosion. In contrast to earlier work, we regress a simple model of exhumation pace through the bedrock ages such that age can vary as a function of both elevation and position. The resulting iso-age surfaces allow us to determine potential source areas for detrital AFT age distributions. We find that old ages observed in detrital distributions are just as likely to be sourced from low-elevation locations that are far from the coast as from high-elevation locations close to the coast. Additional data from lower temperature systems are thus required to make firm conclusions on landscape evolution in the region and distinguish between the two landscape-forming scenarios.

High-elevation low-relief plateau surfaces separated by deeply incised fjords characterize many high-latitude continental margins, including the Antarctic Peninsula, Greenland and Norway. The origin of this characteristic landscape continues to be debated (e.g. Nielsen et al. 2009; Steer et al. 2012; Japsen et al. 2018; Pedersen et al. 2021), with important implications for understanding long-term glaciation and sea-level. In particular, the dimensions of past glaciations provide a record of palaeoclimate but also depend on underlying topography (Kaplan et al. 2009; Anderson et al. 2012; Sternai et al. 2013; Clague et al. 2020). Do the plateau surfaces represent Mesozoic peneplanation and subsequent Cenozoic uplift, or are they formed as a result of glaciation? The former would suggest that the plateau surfaces are eroding very slowly and can be used to map Cenozoic rock uplift, whereas the latter would suggest that they can neither be used to estimate surface uplift nor to place constraints on valley incision. In Greenland, studies using cosmogenic nuclides have found that erosion of these high-elevation plateau surfaces has been slow (Strunk et al. 2017; Andersen et al. 2020; Skov et al. 2020), with cosmogenic nuclide abundances preserved and inherited from previous interglacials. However, landscape-evolution modelling efforts have indicated that such observations are in accordance with abundant long-term erosion on the plateau surfaces (Egholm et al. 2017). Indeed, cold-based non-erosive ice at these high-elevation plateau surfaces is an expected consequence of fjord formation that will localize glacial erosion in valley troughs and result in erosion-driven isostatic uplift of the plateau surfaces.

Detrital thermochronometry has the unique potential to reveal where sand is coming from and how this pattern of erosion has changed through time. A recent study by Olivetti et al. (2022) presented impressive detrital apatite fission-track (AFT) data from offshore South and East Greenland (Ocean Drilling Project sites 918 and 987) to provide unique constraints on late Cenozoic erosion across East Greenland. The principle behind this method is that AFT ages record cooling during exhumation through the upper crust and bedrock ages tend to increase with elevation (Wagner et al. 1979). If the AFT age distribution can be mapped across a bedrock area, detrital AFT ages from glacial and fluvial sediments can be used as a provenance tool to track patterns of erosion (Hurford and Carter 1991; Stock et al. 2006; Vermeesch 2007; Fox et al. 2015). Detrital AFT ages thus provide an indication of locations being eroded as a function of elevation, with relatively old ages at high elevation (Wagner et al. 1979). It is this potential that makes the AFT method ideal to determine whether or not the plateau surfaces have been eroding and contributing to detrital material in the late Cenozoic (Garver et al. 1999). Olivetti et al. (2022) highlighted that relatively old ages are found in detrital material since 8 myr ago, suggesting prolonged erosion at high elevations across East Greenland (Figure 1). They also argued that even older age populations found in stratigraphic sections closer to the present suggest that progressively higher elevations are contributing detrital material.

However, a few complexities obscure the relationships between bedrock AFT age, elevation and detrital thermochronometry. First, a single fission-track age for a bedrock sample comprises several tens of individual crystal ages (pooled sets of track counts). Pooling of grain data to calculate the pooled fission-track age is a means of capturing the stochastic (Poisson) process of fission of U atoms per unit volume. A Poisson distribution for a fission-track age will have a standard deviation of single grain ages (counts), thus a single grain age is a single point estimate within the true age distribution. This means that crystals derived from a specific bedrock location will contribute a wide range of single crystal ages to a detrital distribution. Ultimately, this blurs the mapping of detrital age to a specific location (Fox et al. 2015). However, in the case of Greenland, the pooled ages appear to be very consistent, as highlighted by the reproducible age–elevation relationships (Olivetti et al. 2022). Second, pooled bedrock AFT ages may not vary as a simple function of elevation. If bedrock age varies only as a function of elevation, surfaces of constant age would be flat and horizontal. These surfaces of constant age are termed iso-age surfaces (Vernon et al. 2008; McPhillips and Brandon 2010). Spatial variations in bedrock AFT age for a specific elevation may result from spatial variations in exhumation rates as the rocks passed through the AFT closure temperature in the Mesozoic. Spatial variations might also be the result of regional tilting and warping of previously near-horizontal iso-age surfaces. As noted by Olivetti et al. (2022), bedrock AFT ages in East Greenland become progressively older from the coast for the same elevations, highlighting that ages do not vary simply as a function of elevation. For our purposes, we focus on this second complexity and test whether a simplified method to predict the bedrock age map limits the conclusions of Olivetti et al. (2022). To do this, however, a new interpolation procedure is required.

To understand the source of sediment and to infer spatial patterns in erosion rate from a detrital dataset, the inferred bedrock age distribution needs to be as accurate as possible. Our aim here is therefore to resolve bedrock age distributions and determine potential source areas for the detrital AFT ages, considering that the distribution of bedrock ages across the East Greenland landscape varies as a function of elevation as well as longitude and latitude. Our approach is based on modelling the slope of the age–elevation relationship in space, under the assumption that the slope of this age–elevation relationship varies smoothly. With this approach, we show that bedrock AFT ages that would yield old detrital AFT ages are found both across plateau surfaces and in steep parts of the current topography where glacial erosion might be most efficient.

In the approach we develop here, we regress a spatially variable model of exhumation pace through elevation–age data using a linear inverse approach. First, we highlight how a single age can be converted to a simple summation equation before we link multiple ages in space. We then apply our approach to available AFT bedrock ages across East Greenland.

A plot of thermochronometric age on the x-axis v. elevation on the y-axis is commonly referred to as an age–elevation relationship (AER). A linear model regressed through these data has a slope with units of kilometres per million years, often interpreted as an exhumation rate, and an intercept with units of kilometres, often interpreted as a closure elevation. The closure depth is the difference between the elevation of the sample and the elevation of this closure elevation, which might be below sea-level. In contrast, it is convenient to consider this inverse problem as an elevation–age relationship (EAR) because this ensures that the errors in age are associated with the dependent elevation variable. The slope of this line is exhumation pace (Ma km–1) and the intercept has units of million years. There are additional benefits of using an EAR over an AER as discussed by Fox and Carter (2020).

Formally we can write that a single age, τ, is the integral of exhumation pace, which varies as function of elevation, p(z), between the elevation of the closure isotherm, h0, and the elevation of the sample, h:
τ=h0hp(z)dz.
(1)
This expression can be discretized into blocks of constant elevation (Δz) and written as a summation equation for a single age τi:
τi=j=1MpjΔz+pM+1zR
(2)
where zR is a remainder term so that the distance between the closure elevation and the sample height (commonly referred to as the closure depth) is h – h0 = MΔz + zR and j and M are integers used in this summation. A similar approach has previously been developed, termed GLIDE (Fox et al. 2014a). However, in contrast to GLIDE, we do not calculate the closure depth using a thermal model. Instead, we approach the problem as a geometric problem with the closure depth simply representing a point in elevation. The disadvantage of this is that multiple thermochronometric systems cannot be interpreted simultaneously. An advantage of this is that it reduces the complexity of the model. For example, it is clear that the bedrock ages are the result of a complex thermal history with multiple reheating events (Green et al. 2018; Bernard et al. 2016). By simply regressing a geometric model, these complex thermal histories can be ignored, and our aim is simply to reproduce the bedrock ages. Furthermore, as we are simply interested in the expected ages at different elevations, we prefer this geometric approach and do not interpret exhumation pace in terms of exhumation rate. We therefore fix h0 at a constant elevation that does not vary in space or time. Because we do not interpret the inferred EAR between this elevation and the lowest elevation ages, this choice does not influence our ability to reproduce the data and predict bedrock ages across the landscape.

A discrete expression, similar to equation (2), can be written for each age in a dataset and these can be combined as a matrix–vector product. If the samples are all from the same location, this solution can be solved with additional constraints to ensure the problem is well posed. Here the total number of elevation blocks required is determined by the highest elevation sample or Mmax.

However, if the samples are from different locations, spatial variability is expected. To account for this, we discretize geographical space into a grid of pixels (NxNy) with uniform exhumation paces in each pixel, and at each position the exhumation pace varies with elevation. Now the total number of unknown exhumation pace parameters is given by NxNyMmax. Solving the resulting system of equations requires smoothness constraints. We minimize curvature in space described with the weighting matrix Ws, and curvature in elevation described with the weighting matrix Wt. These matrices calculate the curvature: across the central pixel and the surrounding four pixels for Ws; and the points above and below for Wt. The relative weights, α and λ, applied to these curvatures are determined using a trial-and-error approach to produce a model that fits the bedrock data. We note that this is potentially an ill-posed inverse problem and multiple combinations of model parameters will fit the data equally well. For this reason, we seek the smoothest model that fits the data and is reasonable. Therefore, we solve
(GWsWt)(p)=(τ00).
(3)
Here the top rows attempt to predict the ages by the spatial variability in exhumation pace. The middle rows (Wsp = 0, where 0 is a vector of zeros of length NxNyMmax) attempt to minimize spatial curvature, and the lower rows (Wtp = 0, where 0 is a vector of zeros of length NxNyMmax) minimize curvature of local EARs.

We solve equation (3) for the unknown exhumation pace model using weighting values of α = 500 and λ = 10. We use the bedrock data compiled by Olivetti et al. (2022) in this analysis to be consistent with earlier work. In this way, our goal is simply to accurately reproduce the observed bedrock ages and map bedrock ages across the landscape. We do note, however, that areas of our solution will be less reliable than other parts owing to the density of surface samples. In particular, in the western limits of our solution, the model is less well constrained and this can be inspected by simply looking at the distance from the bedrock samples. Additional data collected in more locations would improve our interpolation.

Solving equation (3) for the unknown exhumation pace vector results in a model that matches the observed EARs locally and varies smoothly between these locations (Fig. 2). This is in contrast to the maps produced by Olivetti et al. (2022), which cannot honour the observation that multiple ages may be found at the same, or very similar, geographical location (latitude and longitude) owing to the fact that age varies as a function of elevation. The most illustrative example of this would result from ages obtained from samples from a truly vertical borehole. Of course, this is an extreme example. In reality, bedrock ages may vary rapidly over short spatial distances, such as with a near-vertical profile, and this may obscure any long-wavelength trends in the ages. Ignoring this observation and simply mapping age as a function of space, as carried out by Olivetti et al. (2022), means that ages have previously been poorly reproduced by the interpolation procedure. In contrast, the AERs produced by Olivetti et al. (2022) reproduce the local ages perfectly, but do not respect the long-wavelength variability in age. Thus with previous methods it is unclear how representative an AER is of an area and how to combine AER constraints in the spaces between the bedrock sample locations. Our approach allows us to account for long-wavelength variation in ages owing to underlying geodynamic processes and short-wavelength variation owing to elevation changes. Furthermore, we are able to reproduce the bedrock ages and predict ages at all surface locations. An ability to predict ages at all surface locations is a requirement because the detrital ages integrate bedrock data across space.

The exhumation pace maps are relatively stable through elevation except where changes in the slope of the AER have previously been detected. We show three representative maps here (Fig. 3). We stress that these maps should not be interpreted in terms of an exhumation rate history because the exhumation pace parameters do not constrain exhumation rate over the same time intervals. Furthermore, we have not accounted for the topographic effects on the elevation of the closure isotherm or changes in the advection of heat owing to exhumation. At high elevations, the pace maps show relatively high pace and slow rate across most of the area. At lower elevations, the pace maps show lower values towards the south.

Finally, we can predict age as a function of elevation and geographical position. To do this we solve equation (2) for each elevation point in a digital elevtion model using the recovered exhumation pace from our analysis of the bedrock data. Figure 4 shows topography (a), topographic slope (b) and predicted age (c). We see that old ages are found at high elevations and further from the coast, in agreement with earlier work. However, here we are able to interpolate and extrapolate between datapoints because of the simplicity of the model.

Olivetti et al. (2022) argued that the old age populations they found within offshore detrital samples indicate that the high-elevation plateau surfaces are producing detrital material as far back as 8 Ma. This is contentious because it raises the question of whether these are formed in situ or can be used to map uplift of a low-relief surface in space. We can highlight that the old ages might simply be coming from high-slope parts of the landscape simply by hiding parts of Figure 4 that have ages that would not contribute to this old population. We define a mask as ages between 250 and 300 Ma. This corresponds to the ‘old’ age population identified by Olivetti et al. (2022). Figure 5 shows that the old population could correspond to both high and low elevations (Fig. 5a), as well as low- and high-slope areas (Fig. 5c). This is consistent with the idea that sediment is produced where glaciers are sliding fastest (Herman et al. 2015). Indeed, to determine where the sediment is coming from, additional data are required.

It is also possible to estimate the probability of erosion based on the distribution of detrital ages and the bedrock maps. To do this, we fit a kernel density estimate through the detrital ages that is proportional to probability (Vermeesch 2012). Here we use the youngest detrital sample from Site 987, which has a depositional age of 0.7 ± 0.1 Ma (Fig. 6a). We then map the age probabilities onto the landscape to produce a simple map of probability of erosion (Fig. 6b). The most striking feature of this map is the high probability of erosion from very low elevations that are below sea-level. These are highlighted by the red colours at low elevations. This is challenging to explain in terms of modern erosional processes because we would expect the base of deep glacial fjords to be depositional environments. However, given that this sample is from sediment deposited at 0.7 Ma it makes sense that the high probability of erosion in the base of the fjords could be the result of glacial erosion and deepening when glaciers occupied the fjords (Sternai et al. 2013; Fox et al. 2014b; Pedersen et al. 2014). Interestingly, high erosion probability values are predicted within the overdeepened fjord even during these late stages of glaciation, whereas it might be expected that the focus of erosion would migrate headwards through time (Shuster et al. 2011). We also see high probabilities of erosion across higher elevation parts of the catchment that are potentially incising where old bedrock ages are found. This is more consistent with the hypothesis that valley erosion would migrate headwards through time. Across the high-elevation low-relief surfaces, we mostly see low probability of erosion.

Bedrock ages with different spatial patterns would be required to leverage the detrital dataset to determine whether old ages are coming from low-relief, high-elevation surfaces or incising valleys. In this way, the production function (where sediment is coming from) is the same but the bedrock age function is different. The apatite (U–Th)/He system with a lower closure temperature would be ideal because the ages will be controlled by topographic evolution as opposed to earlier tectonic events. Furthermore, radiation damage influences the closure temperature of this system so that damaged crystals have higher closure temperatures than less damaged crystals. With sufficient data, the apatite (U–Th)/He system combined with radiation damage effectively provides a suite of slightly different thermochronometric systems. Therefore, a specific detrital age fingerprints different locations depending on the amount of radiation damage. Physically based models accounting for heat flow, topographic evolution and damage would be required for this sort of analysis (Clinger et al. 2020, 2022). Furthermore, accounting for these complexities, along with collecting additional data, would allow additional processes to be investigated. These include whether sediment is mixed during glacial transportation (Bernard et al. 2020), how lithology influences glacial erosion (Bernard et al. 2021) or how overdeepened basins trap sediment and nutrients (Cook and Swift 2012; Delaney and Adhikari 2020; Swift et al. 2021).

By interpreting bedrock data with simple tools to map age as a function of space, the bedrock age distribution can be interpolated and extrapolated. We note that aspects of our model, particularly towards the west, may be poorly resolved owing to large distances between bedrock ages. Ultimately, the only way to improve this is to collect additional bedrock data. We have used a geometrical model because we simply want to reproduce the bedrock age distribution, but models that solve heat transport could also be used. In either case, earlier simplified models relating age to elevation may yield conclusions that are dependent on a model that is too unrealistic.

Our ability to estimate the source locations of sediment offshore of inaccessible locations has important implications for understanding landscape evolution. If the low-relief surfaces are not eroding and have not been eroding over the last 8 myr, as suggested by our analysis, these surfaces can be used to estimate the magnitudes of valley incision and patterns of geodynamic processes. This also implies that glaciers in Greenland have cut deep fjords into a relatively low-relief landscape that existed prior to glaciation. Ultimately, such a reconstructed low-relief landscape could be used to simulate glacial conditions at the time of glacial inception with improved fidelity.

We have shown that the current combination of bedrock AFT data and detrital AFT data cannot determine whether the low-relief, high-elevation surfaces of Greenland have been eroding for the last 8 myr. By interpolating a simple model through available data, we are able to predict bedrock data and extrapolate the exhumation pace between data points. Resulting predicted bedrock AFT ages highlight that old ages are found in low-elevation and relatively steep parts of the landscape, in addition to the high-elevation plateau surfaces that have previously been put forward as their source. At these locations, ice is expected to be thick and fast flowing, potentially producing high rates of glacial erosion. Propagating uncertainties from the measured ages, interpolation model and detrital distributions would simply make conclusions drawn from the data less discriminating and this is not required here.

We would like to thank V. Olivetti for sharing the detrital dataset and providing helpful comments on the paper.

MF: conceptualization (equal), formal analysis (equal), funding acquisition (lead), methodology (equal), writing – original draft (equal); VP: investigation (equal), writing – review & editing (equal)

This work was funded by the Natural Environment Research Council (NE/N015479/1).

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.

All data generated or analysed during this study are included in this published article.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)