Surface paleothermometry using low-temperature thermoluminescence of feldspar

Thermoluminescence (TL) of feldspar is investigated for its potential to extract temperature histories experienced by rocks exposed at Earth's surface. TL signals from feldspar observed in the laboratory arise from the release of trapped electrons from a continuous distribution of trapping energies that have a range of thermal stabilities. The distribution of trapping energies, or thermal stabilities, is such that the lifetime of trapped electrons at room temperature ranges from less than a year to several billion years. Shorter lifetimes are associated with low-temperature TL signals, or peaks, and longer lifetimes are associated with high temperature TL signals. Here we show that trapping energies associated with shorter lifetimes, or lower-temperature TL signals (i.e. between 200 and 250 ^∘C), are sensitive to temperature fluctuations occurring at Earth's surface over geological timescales. Furthermore, we show that it is possible to reconstruct past surface temperature histories in terrestrial settings by exploiting the continuous distribution of trapping energies. The potential of this method is first tested through theoretical experiments, in which a periodic temperature history is applied to a kinetic model that encapsulates the kinetic characteristics of TL thermometry. We then use a Bayesian approach to invert TL measurements into temperature histories of rocks, assuming that past temperature variations follow climate variations observed in the δ¹⁸O records. Finally, we test the approach on two samples collected at the Mer de Glace (Mont Blanc massif, European Alps) and find similar temperature histories for both samples. Our results show that the TL of feldspar may be used as a paleothermometer.

1 Introduction

Earth's climate fluctuates in a cyclic way, on timescales of years to millions of years, driven by Earth's orbital processes and rare aberrant shift and extreme climate transients during the last 10³ to 10⁵ years (e.g. Zachos et al., 2001). Reconstructions of past terrestrial climates often rely on the use of climate proxies that preserve the physical and/or chemical characteristics related to Earth's past climate. Examples of such proxies include ice cores, tree rings, sub-fossil pollen, boreholes, corals, lake sediments and carbonate speleothems (e.g. Jones and Mann, 2004 for a review). Although they have provided invaluable insights into past climates and their physical characteristics, very few of these proxies provide a direct measure of temperature variations in continental settings (e.g. glycerol dialkyl glycerol tetraether, GDGT; Tierney et al., 2012), and many of these methods often suffer from methodological limitations that limit reliable construction of terrestrial temperatures. For example, fossil pollen and plant macrofossils have provided key insights into past terrestrial climate at millennial timescales (e.g. Bartlein et al., 2011) but typically rely on numerous additional climate parameters, including precipitation, plant-available moisture, seasonality and the length of the growing season, which make the inference of rock temperature histories challenging. Similarly, GDGTs rely on the preservation of organic compounds. To address some of these issues, Tremblay et al. (2014a, b) recently introduced a new paleotemperature proxy based on the thermal stability and release of ³He and ²¹Ne noble gases in quartz. The system has a single energy level and is therefore only able to estimate a single equivalent diffusion temperature (EDT). As a result, reconstruction of complex rock temperature histories is challenging.

This study revives the use of thermoluminescence (TL) as a paleothermometer and takes advantage of recent progress made for TL thermochronometry (Biswas et al., 2018) to introduce a new approach for reconstructing temperature histories in terrestrial settings at millennial timescales. The idea is not new; it has already been tested through a series of feasibility studies (Li and Li, 2012; Ronca, 1964; Ronca and Zeller, 1965) and was applied to lunar samples (Durrani et al., 1977). Unfortunately, the method has remained underdeveloped since this early work. More recently, Guralnik and Sohbati (2019) used optically stimulated luminescence (OSL) to estimate paleotemperature. However, they only investigated a single thermal stability, and thus a single equivalent temperature, similar to the noble gas approach of Tremblay et al. (2014a, b).

The present study aims to establish the TL of feldspar as a paleothermometer to constrain temporal variations of Earth's surface temperature, within the last few tens of thousands of years. TL signals from feldspar arise from a series of traps that have different thermal stabilities, with lifetimes at room temperature ranging from less than a year to billions of years (Biswas et al., 2018). In the following, we first investigate the sensitivity of the TL signals to linear and periodic thermal histories. Then, we use a Bayesian approach to infer the thermal histories of rocks from TL measurements. Finally, we apply this method to two samples collected from the Mer de Glace glacier (Mont Blanc massif, European Alps) as proof of concept.

2 Theoretical background

In its simplest form, luminescence enables the measurement of the concentration of charge trapped in the impurity or defect centres of natural crystalline minerals (e.g. quartz or feldspar); the trapped charge population is proportional to the time elapsed since the material was either heated or exposed to sunlight. In the laboratory when the minerals are heated using a linear heating rate, trapped electrons are released from increasingly higher energy traps, some of which recombine radiatively with the holes, and the luminescence process generates a TL glow curve. Trapping of electrons happens in nature when samples are exposed to environmental radiation over geological time. Detrapping reflects the escape of electrons from traps, which can occur thermally or athermally. In nature, the TL system can reach a dynamic equilibrium between radiation-induced growth and decay via thermal and athermal pathways. For higher ambient temperatures, the thermal decay rate is greater, which results in a lower equilibrium level and vice versa (for more details see Aitken, 1985, and Ronca, 1964). Such a difference in equilibrium can be exploited to reconstruct the thermal history of rocks (Biswas et al., 2018; Brown et al., 2017; Guralnik and Sohbati, 2019; Herman and King, 2018; Herman et al., 2010; King et al., 2016a, b). Until now, this principle has only been exploited for inferring rock cooling due to exhumation. Here we use the same principle for inferring the paleotemperature of rocks exposed at Earth's surface. In this section, we start by outlining the theoretical model that describes the TL process of feldspar and how one can constrain all the model parameters. Then we assess how sensitive the model is to various surface temperature histories and ultimately show how surface temperatures in terrestrial settings can influence TL signals.

Figure 1 Inferred TL kinetic parameters from a TL glow curve of sample MBTP9. The method used to constrain these parameters in the laboratory is fully explained in Biswas et al. (2018) and described in Sect. S1.

2.1 TL kinetic model and TL thermometer

The kinetic model that describes the TL of feldspar was reviewed in Biswas et al. (2018). The model encapsulates the process of populating traps with electrons in response to surrounding radiations, and the processes of electron escape through thermal and athermal pathways. A key element of the method is that each kinetic parameter can be constrained in the laboratory for each sample. We briefly outline the kinetic model here and illustrate how the kinetic parameters may be constrained from laboratory measurements for one sample. The reader is referred to Sect. S1 in the Supplement and Biswas et al. (2018) for further information.

The model assumes general-order kinetics (Biswas et al., 2018; Guralnik et al., 2015a, b). The rate equation for the trapped charge population of a specific trapping centre, with single trap depth (E), frequency factor (s) and athermal fading parameter (ρ^′; Huntley, 2006; Tachiya and Mozumder, 1974) is as follows:

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(\stackrel{\mathrm{‾}}{n}\right(r^{\prime },t\left)\right)=& {\displaystyle \frac{\dot{D}}{D_{\mathrm{0}}}}(\mathrm{1}-\stackrel{\mathrm{‾}}{n}(r^{\prime },t){)}^a-s{e}^{-\frac{E}{kT}}(\stackrel{\mathrm{‾}}{n}(r^{\prime },t){)}^b\\ \text{(1)}& {\displaystyle }& {\displaystyle }-\stackrel{\mathrm{̃}}{s}{e}^{-{{\mathit{\rho }}^{\prime }}^{-\frac{\mathrm{1}}{\mathrm{3}}}{r}^{\prime }}\stackrel{\mathrm{‾}}{n}(r^{\prime },t),\end{array}$$

where $\stackrel{\mathrm{‾}}{n}$ is equal to n∕N (where n is the number of trapped electrons at time t and temperature T, and N is the total number of available traps), $\dot{D}$ is the dose rate due to ambient radioactivity (Gy ka⁻¹), D₀ is the onset of dose saturation (Gy), a and b are the kinetic orders of trapping and thermal detrapping respectively, E is the trap depth or activation energy (eV), s and $\stackrel{\mathrm{̃}}{s}$ are the thermal and athermal frequency factor respectively (s⁻¹), ρ^′ is the dimensionless athermal fading rate, and r^′ is a dimensionless distance that characterizes the probability of athermal escape (Huntley, 2006). Each of these parameters can then be constrained from laboratory experiments that are described in Biswas et al. (2018) and summarized in Table S1 in the Supplement for one sample.

Figure 2 Panels (a–c) are the three prescribed thermal histories: isothermal, cooling and warming respectively. Panels (d–f) are the corresponding dynamic equilibrium levels of trapped charge population ($\stackrel{\mathrm{‾}}{n}$) of 10 different thermometers in the range of 200–300 ^∘C. The temperatures of the thermometers (or TL temperature) are shown on the right.

To account for athermal loss, i.e. anomalous fading (Wintle, 1973), the total number of trapped electrons at any instant $\stackrel{\mathrm{‾}}{n}\left(t\right)$ is obtained by integrating over the whole range of dimensionless distances ($\mathrm{0}<r^{\prime }<\mathrm{2}$; Kars et al., 2008) over which electrons can athermally escape;

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{‾}}{n}\left(t\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}p\left(r^{\prime }\right)\stackrel{\mathrm{‾}}{n}(r^{\prime },t)\mathrm{d}{r}^{\prime },\end{array}$$

where p(r^′) is the probability of nearest recombination centre at a distance between r^′ and $r^{\prime }+\mathrm{d}{r}^{\prime }$ and expressed as $p\left(r^{\prime }\right)\mathrm{d}{r}^{\prime }=\mathrm{3}{r^{\prime }}^{\mathrm{2}}{e}^{{r^{\prime }}^{\mathrm{3}}}\mathrm{d}{r}^{\prime }$ (Huntley, 2006). This model was validated using rocks from the KTB borehole and applied to samples from Namche Barwa (Biswas et al., 2018), which gave results in agreement with other studies from the same area (King et al., 2016a).

The TL of feldspar arises from a continuous distribution of trapping energies (Biswas et al., 2018; Duller, 1997; Grün and Packman, 1994; Pagonis et al., 2014; Strickertsson, 1985) and can be assumed to be the sum of a large number of traps (Biswas et al., 2018; Pagonis et al., 2014); all follow the process described by Eq. (1). To constrain the kinetic parameters in Eq. (1), we measure full TL glow curves and see how the kinetic parameters are distributed along the TL glow curve (Biswas et al., 2018). For the modelling, we use the kinetic parameters for sample MBTP9 (Lehmann et al., 2020, for sample details). The experimental details are provided in Sects. 4 and S1. The distribution of kinetic parameters along the TL glow curve temperature is reported in Fig. 1. The results show that the kinetic parameters vary systemically with the glow curve temperature. Such data are then fitted using a spline function from which the kinetic parameters are then extracted for a specific TL temperature (Biswas et al., 2018), which we define herein as a specific “TL thermometer”.

2.2 Temperature response of TL thermometers

In this section, we use the model (Eq. 1) in a forward manner by prescribing a temperature history and predicting the trapped charge population through time ($\stackrel{\mathrm{‾}}{n}$). Both linear (isothermal, warming and cooling) and periodic thermal histories are used. The model was run for 1 Myr, with an initial condition of $\stackrel{\mathrm{‾}}{n}=\mathrm{0}$, which is long enough to ensure that $\stackrel{\mathrm{‾}}{n}$ reaches equilibrium. We investigate how $\stackrel{\mathrm{‾}}{n}$ changes for 10 different TL thermometers, in the temperature range of 200–300 ^∘C with 10 ^∘C intervals, each having an independent set of kinetic parameters (see Table S1).

2.2.1 Linear thermal history

The thermal response of the dynamic equilibrium level of the trapped charge population $(\stackrel{\mathrm{‾}}{n}$) of the 10 thermometers is tested for three different linear thermal histories: (1) isothermal holding at 20 ^∘C (Fig. 2a), (2) isothermal holding at 20 ^∘C for 900 kyr followed by linear cooling to 0 ^∘C during the last 100 kyr (Fig. 2b) and (3) isothermal holding at 0 ^∘C for 900 kyr followed by linear heating to 20 ^∘C during the last 100 kyr (Fig. 2c). The results are shown in Fig. 2d, e and f respectively. In all cases, $\stackrel{\mathrm{‾}}{n}$ is lowest for the lowest TL thermometers (or TL signals) and highest for the highest TL thermometers. For the isothermal scenario, $\stackrel{\mathrm{‾}}{n}$ remains constant over the entire recent time period (after reaching steady state) as there is no temperature change. For the cooling scenario, $\stackrel{\mathrm{‾}}{n}$ increases as temperature decreases, because of a decrease in thermal loss. For the warming scenario, $\stackrel{\mathrm{‾}}{n}$ decreases as temperature increases because of an increase in thermal loss. The increase and decrease in $\stackrel{\mathrm{‾}}{n}$ (Fig. 2d and e) with temperature are most pronounced for lower-temperature TL thermometers (<250 ^∘C) and negligible for higher-temperature TL thermometers (>250 ^∘C). It must be noted that if the amplitude of temperature change (here 20 ^∘C) were increased, the higher-TL thermometers would change more dramatically. However, this would be unrealistic and beyond temperature variations observed at Earth's surface at comparable timescales.

The previous section shows that the thermal sensitivities of the different TL thermometers are distinct. It is therefore expected that their temporal sensitivities are also different. To quantify the temporal response of each TL thermometer, we prescribe a simple step function in which the temperature is set to be equal to 0 ^∘C until a given time, t_change, and then increases to 10, 20 and 30 ^∘C, in three different cases. We then calculate the present-day trapped charge population (${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$) for each TL thermometer. t_change is varied from 100 kyr to the present.

The model predictions (${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$) are shown as a function of time of change (t_change) in Fig. 3. In addition, we calculate the t_change where ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ increases by 20 % compared to ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ predicted for a thermal history where t_change is equal to 100 kyr. We define this corresponding time as the memory time (t_memory) at which the TL thermometer can record a temperature change. As expected, the lower-TL thermometers record more recent temperature changes and the higher-TL thermometers record older temperature changes. For example, in the case of a 10 ^∘C temperature change (Fig. 3a), the 200–210 ^∘C thermometers record a change for ∼10 kyr, while the 240–250 ^∘C thermometers record a change for ∼50 kyr. Furthermore, the number of sensitive thermometers increases as we raise the final temperature. For a 10 ^∘C temperature change, only 5 TL thermometers (200–250 ^∘C) record the temperature change (Fig. 3a), 6 TL thermometers (200–250 ^∘C) record a 20 ^∘C temperature change (Fig. 3b), and 7 TL thermometers (200–250 ^∘C) record a 30 ^∘C temperature change (Fig. 3c).

Figure 3 The variation of present-day trapped charge population (${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$) of different thermometers with the times of change of temperature (t_change) of step functions (such as thermal history) with temperature changes of (a) 10 ^∘C, (b) 20 ^∘C and (c) 30 ^∘C. The asterisk symbols denote the memory time (t_memory) that a thermometer can record the temperature change history. Estimation of t_memory is illustrated in plot (a) for 200–210 ^∘C TL thermometer. The solid green line corresponds to ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ for t_change=100 ka, and the dashed green lines represent ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ for t_change=t_memory, which is 20 % higher than the solid green line.

Figure 4 Variation of trapped charge population ($\stackrel{\mathrm{‾}}{n}$) for different sinusoidal thermal histories. Panels (a–c) represent the prescribed thermal histories with the same mean amplitude (10 ^∘C) and three different mean temperatures (0, 15 and 30 ^∘C) but for different periods, 1, 10 and 100 ka respectively. The dashed lines are the isotherm (mean temperature of oscillation). Panels (d–f), (g–i) and (j–l) are the response of $\stackrel{\mathrm{‾}}{n}$ for the corresponding thermal field. The solid lines are for oscillating fields and dashed lines are for isothermal fields. The temperatures of the representative thermometers (or TL temperature) are shown inside in blue.

Figure 5 Present-day trapped charge populations (${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$) of different TL signals (or thermometers) with (a) mean temperature variation (amplitude and period are fixed), (b) amplitude variation (mean temperature and period are fixed) and (c) period variation (mean temperature and amplitudes are fixed).

2.2.2 Periodic thermal history

Earth's climate varies in cyclical way, at multiple timescales from years to decades, centuries, and millennia, influenced by periodic variations in the Earth's orbit, known as Milankovitch cycles, at 25.77, 41 and 100 kyr. Therefore, temperatures are dictated by periodic functions that include several harmonics comprising decadal and millennial periods. In order to assess how the trapped charge population is affected by a periodic temperature history, we prescribe the following function as a thermal history.

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }T\left(t\right)=T_{\text{mean}}+T_{\text{amp}}\times \text{sin}\left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{P}}t\right),\end{array}$$

where T_mean, T_amp and P are the mean temperature, the amplitude and the period of oscillation respectively.

Different combinations of arbitrary periodic thermal histories, with the same amplitude (10 ^∘C) but in three different periods (1, 10 and 100 kyr) and with three different mean temperatures (0, 15 and 30 ^∘C), were used (solid lines of Fig. 4a, b and c) to assess the effect of periodic temperature variations on $\stackrel{\mathrm{‾}}{n}$ (${\stackrel{\mathrm{‾}}{n}}_{\text{osc}}$) for each TL thermometer (solid lines of Fig. 4d–l). These predictions are compared to isothermal effects on $\stackrel{\mathrm{‾}}{n}$ (${\stackrel{\mathrm{‾}}{n}}_{\text{iso}}$) (dashed line in Fig. 4d–l) if the samples are kept at the mean temperature of the corresponding oscillation (dashed lines of Fig. 4a, b and c).

The results show that the $\stackrel{\mathrm{‾}}{n}$ always depletes more for the lower temperature TL thermometer (e.g. 200–210 ^∘C TL) than for the higher-temperature TL thermometer (e.g. 290–300 ^∘C TL), which results from a gradient in the thermal stabilities of lower- to higher-temperature TL thermometers. For every TL thermometer, $\stackrel{\mathrm{‾}}{n}$ decreases if the mean temperature, T_mean, increases (Fig. 4d–f, g–i and j–l). This is because the probability of detrapping increases with increasing temperature. Finally, the higher-temperature TL thermometers (near to 300 ^∘C) remain relatively insensitive to such periodic temperature forcing (T_mean up to 30 ^∘C); with increasing T_mean the higher-temperature TL thermometers become more responsive.

One can also describe how the system behaves by comparing the period of oscillation (P) and the lifetime (or resident time) of trapped electrons (τ) for a given temperature. For the 200 to 300 ^∘C TL thermometers, τ spans ∼10 kyr to 1 Gyr when the samples are at 0 ^∘C. At higher temperatures τ is reduced as the probability of electron escape increases, reducing the lifetime to between ∼0.1 kyr and 1 Myr when the samples are stored at 30 ^∘C (Biswas et al., 2018). For P≪τ (e.g. Fig. 4d, where P=1 kyr and τ spans ∼10 kyr to 1 Gyr for the 200 to 300 ^∘C TL thermometers for T_mean=0 ^∘C), the value of ${\stackrel{\mathrm{‾}}{n}}_{\text{osc}}$ exhibits small fluctuations but always remains lower than ${\stackrel{\mathrm{‾}}{n}}_{\text{iso}}$. This result implies that smaller periods (<1 kyr) and T_mean (<30 ^∘C) do not influence trapped charge equilibrium levels of 200 to 300 ^∘C TL thermometers in an oscillating fashion and cannot be differentiated from the trapped charge population resulting from an isothermal condition. We must mention here, if the amplitude of oscillation increases, the oscillating response to trapped charge equilibrium levels will be relatively prominent. However, the predicted values of ${\stackrel{\mathrm{‾}}{n}}_{\text{osc}}$ are lower than those predicted using constant temperature of T_mean, which is the mean temperature of the oscillation explored. Similarly, the present-day ${\stackrel{\mathrm{‾}}{n}}_{\text{osc}}$ remains indistinguishable from ${\stackrel{\mathrm{‾}}{n}}_{\text{iso}}$ when P≫τ (e.g. see the behaviour of the low-temperature TL thermometer shown in Fig. 4l). These two endmember scenarios are therefore not suitable for predicting the temporal variation of surface temperature. Interestingly, the response of ${\stackrel{\mathrm{‾}}{n}}_{\text{osc}}$ deviates from its temperature forcing when P∼τ (e.g. Fig. 4g–i and j–l). Under this condition, ${\stackrel{\mathrm{‾}}{n}}_{\text{osc}}$ is out of phase and asymmetric compared to the prescribed forcing, i.e. the thermal history. More importantly, the degree of deviation for different thermochronometers is different. Therefore, temperature variations can be reconstructed by targeting TL thermometers that have lifetimes of trapped electrons comparable to the period of surface temperature changes.

As discussed in the previous section, the response of trapped electron concentrations corresponding to a TL thermometer depends highly on the three characteristic parameters of the periodic forcing, i.e. T_mean, T_amp and P. We now test the sensitivity of the model to these three parameters. The present-day trapped charge population (${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$) is predicted for different arbitrary combinations of T_mean (0, 15 and 30 ^∘C), T_amp (5, 10 and 20 ^∘C) and P(1, 10 and 100 kyr). The results show that ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ is highly dependent on the mean temperature variation, and less dependent on the amplitude and the period (Fig. 5). Although the ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ is less sensitive to the amplitude and the period, the pattern of ${\stackrel{\mathrm{‾}}{n}}_{\text{present}}$ of different thermometers is unique. This ensures that complex thermal histories, that comprise multiple harmonics with periods of about tens of thousands of years but that are distinct from one another, can be reconstructed.

3 Inversion of TL data into realistic thermal histories

The objective of this section is to test whether a temperature history can be recovered by inverting TL data into a realistic thermal history. We start the exercise by predicting TL data using Eq. (1) for a specific thermal history (forward modelling) that we then invert using a Bayesian approach (Biswas et al., 2018; King et al., 2016a, b) (inverse modelling). This synthetic experiment is performed in two different cases. First, we describe the forward and inverse modelling in a general way and then report the two synthetic examples.

3.1 Forward modelling

Forward modelling is achieved by solving Eq. (1) and prescribing a thermal history, similarly to the previous sections. This approach enables us to predict the present-day trapped charge population for a specific TL thermometer using the kinetic parameters extracted for sample MBTP9. Using this approach, we generate a range of “observed values” (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) for a particular thermal history, which we then try to recover using an inversion method. We run the model for 1 Ma to ensure that $\stackrel{\mathrm{‾}}{n}$ reaches steady state assuming an initial condition of $\stackrel{\mathrm{‾}}{n}=\mathrm{0}$. For a specific thermal history, TL thermometers with lower thermal stability exhibit lower ${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$ than TL thermometers with higher thermal stabilities. It is worth noting that the predicted ${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$ values are mostly sensitive to two parameters, the trap depth (E) and athermal fading (ρ^′), such that these must be constrained carefully from laboratory experiments.

3.2 Inverse modelling

To invert TL data (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) into a thermal history, a Bayesian approach is used. We first generate a large number of random thermal histories (300 000). For each random path, the present-day TL signals are predicted by solving Eq. (1). Each predicted present-day TL value (${\stackrel{\mathrm{‾}}{n}}_{\text{predict}}$) is then compared with the observed TL (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) using the following misfit function (Wheelock et al., 2015) and likelihood:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\text{misfit}={\displaystyle \frac{\mathrm{1}}{l}}\sum _{i=\mathrm{1}}^l{\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\times {\displaystyle \frac{{\stackrel{\mathrm{‾}}{n}}_{\text{obs}}}{{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{n}}_{\text{obs}}}}}\times \log {\displaystyle \frac{{\stackrel{\mathrm{‾}}{n}}_{\text{predict}}}{{\stackrel{\mathrm{‾}}{n}}_{\text{obs}}}}\right]}^{\mathrm{2}},\text{(5)}& {\displaystyle }& {\displaystyle }\text{likelihood}=\text{exp}(-\text{misfit}),\end{array}$$

where l is the number of TL thermometers (here l=4) and ${\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{n}}_{\text{obs}}}$ is the uncertainty of corresponding ${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$. An arbitrary uncertainty of 20 % of ${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$ is assumed for synthetic test. For each random path, ${\stackrel{\mathrm{‾}}{n}}_{\text{predict}}$ for a specific thermometer is usually calculated using mean values of the specific set of kinetic parameters (Biswas et al., 2018; King et al., 2016a, b). However, the kinetic parameters have uncertainties, as shown in Fig. 1 and Table S1. To accommodate the measurement uncertainties in kinetic parameters, for each random thermal history, we randomly picked the kinetic parameters within its error range. Finally, ${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$ was also randomly picked within its error range, assuming that any value within error limit is equally probable (cf. Guralnik et al., 2015a). Since the randomization is applied to a large number of parameters, it is necessary to run the model for large number of thermal histories (300 000 iterations here).

The thermal histories that best fit the data are selected using a rejection algorithm that satisfies the criterion likelihood>R, where R is a random number between 0 and 1. A probability density distribution is then constructed by counting the number of accepted thermal histories passing through each grid cell, which is generated by dividing the time–temperature space (0–100 ka and −50 to 50 ^∘C) into 100×100 cells. This approach is commonly used in different thermochronometric studies (Biswas et al., 2018; Braun et al., 2012; Gallagher et al., 2009; King et al., 2016b). It should be noted that the misfit function (Eq. 4) used here is different to the one used in previous studies (Biswas et al., 2018; King et al., 2016b) but is the same as that used in King et al. (2020). We find that a log misfit enables us to better fit data that vary across orders of magnitude, as trapped charge populations vary greatly for different TL signals.

3.3 Synthetic approach 1

We choose three arbitrary periodic thermal histories, Path1 (T_mean=10 ^∘C, T_amp=10 ^∘C and P=25.77 kyr), Path2 (T_mean=20 ^∘C, T_amp=10 ^∘C and P=25.77 kyr) and Path3 (T_mean=10 ^∘C, T_amp=20 ^∘C and P=25.77 kyr), as shown in Fig. 6a. For each thermal history, the present-day trapped charge concentrations (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) are calculated for four TL thermometers (210–250 ^∘C, 10 ^∘C interval) as described in Sect. 3.1 and represented in Fig. 6b–d. We then invert the TL data (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) into a thermal history as described in Sect. 3.2. For the inverse modelling, we first generate a large number of random periodic histories (300 000), with T_mean and T_amp randomly varying from 0 to 50 ^∘C and P randomly varying between three cycles, 25.77, 41 and 100 kyr. We do not vary P in a completely random fashion because ${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$ is less sensitive to P; i.e. it is difficult to resolve neighbouring periods (as discussed in Sect. 2.2.2 and Fig. 5). The results are shown in Fig. 6e–j. Although this approach predicts the very recent temperature well (up to max 5 ka) it loses the periodic information (25.77 kyr) because of the significant number of accepted thermal histories with different periods (41 and 100 kyr). The same exercise was repeated but the period was fixed to 25.77 kyr for the inversion. The results are shown in Fig. 6k–p. Interestingly, this approach enables the actual solution to be recovered within 1σ uncertainty. This shows that a periodic thermal history can be predicted well if the period is known a priori; it enables T_mean and T_amp to be constrained satisfactorily. To circumvent the limitation (period) of this method we use synthetic approach 2, in which the δ¹⁸O data impose the shape of the thermal histories as a priori information. This is typically done for inversion problems when appropriate. We then constrain T_mean and T_amp of the spectrum (as discussed in the next section).

Figure 6 Results of synthetic experiment of approach 1. Panels (a–d) are for the forward modelling and (e–p) are for the inverse modelling. Panel (a) shows three arbitrary periodic thermal histories with different mean temperature (T_mean) and amplitude (T_amp) but fixed period P. Panels (b–d) are the evolution of trapped charge population ($\stackrel{\mathrm{‾}}{n}$) for four different thermometers (210–250 ^∘C, 10 ^∘C interval) of the corresponding three thermal histories, respectively. The present-day trapped charge populations are considered as observed values (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) for the inverse modelling. Panels (e–g) are the inferred probability density plots when T_mean, T_amp and P are randomly varied. Panels (h–j) depict the fit between the observed TL (obtained through forward modelling). The solid red lines show the predicted median; white lines and black lines show the 1σ and 2σ confidence intervals in the probability density distribution. Panels (k–m) are the inferred probability density functions when T_mean and T_amp are randomly varied but P is fixed. Panels (n–p) depict the fit between the observed TL (obtained through forward modelling). These results are obtained using the kinetic parameters of sample MBTP9.

3.4 Synthetic approach 2

Here we report the result of three thermal histories assuming that the temperature follows the measured δ¹⁸O records from Greenland for the past 60 kyr, which is based on various records from the DYE-3, the GRIP and the NorthGRIP ice cores (Svensson et al., 2008). For our purpose, we scale the δ¹⁸O records to thermal histories and assume a constant temperature prior to 60 ka (Fig. 7a). For each thermal history, the present-day trapped charge concentrations (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) are calculated for four TL thermometers (210–250 ^∘C, 10 ^∘C interval). The results reported in Fig. 7b–d show a clear depletion of the trapped charge population during the last 20 kyr, for all investigated scenarios. However, high-frequency temperature variations are dampened, implying that TL thermometry is insensitive to short-term variations. This result is consistent with the results shown in Fig. 4. For the inverse modelling, we first generate a large number of random periodic histories (300 000), assuming that they all follow the Greenland ice core δ¹⁸O record (Svensson et al., 2008), which we scaled randomly by varying the amplitude of the temperature oscillation (i.e. the difference between minimum the temperature, which is at ∼20 ka, and the maximum temperature at present) between 0 and 40 ^∘C and the minimum temperature (i.e. temperature at ∼20 ka) from −20 to 30 ^∘C (see Sect. S2). Note that making this assumption is somewhat equivalent to assuming a prior estimated on the inferred thermal history (Tarantola, 2005). The inversion results for three tested thermal histories are shown in Fig. 7e–p. The probability density functions (Fig. 7e, f and g) show it is possible to recover all three thermal histories within the 1σ confidence level using this inversion approach.

Figure 7 Results of synthetic experiment of approach 2. Panels (a–d) are for the forward modelling and (e–p) are for the inverse modelling. Panel (a) shows three arbitrary temperature histories obtained by scaling the Greenland δ¹⁸O ice core record (Svensson et al., 2008). Panels (b–d) are the evolution of trapped charge population ($\stackrel{\mathrm{‾}}{n}$) for four different thermometers (210–250 ^∘C, 10 ^∘C interval) corresponding to three histories respectively. The present-day trapped charge populations are considered as observed values (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) for the inverse modelling. Panels (e–g) are the inferred probability density functions. The dashed red lines show the actual temperature history as used in forward modelling, the solid red lines show the predicted median, and white lines and black lines show the 1σ and 2σ confidence intervals in the probability density distribution. Panels (h–j) depict the fit between the observed TL (obtained through forward modelling) and modelled TL (obtained through inverse modelling). Panels (k–m) represent histograms of the parameter, base temperature, T_base (which is temperature at 20 ka). Panels (n–p) represent histograms of present temperature, T_present (which is T_base+T_amp as shown in Eq. S9). These results are obtained using the kinetic parameters of sample MBTP9.

4 Proof of concept

The following sections explore the potential of multiple TL thermometers in the lower-temperature region of the TL glow curve (210–250 ^∘C) to infer the rock temperature histories for two samples collected in the French Alps.

4.1 Sample location

Two bedrock samples (MBTP1 and MBTP9) were collected at the Mer de Glace glacier (Mont Blanc massif, European Alps) at an altitude of 2545 and 2133 m. The rock surfaces were exposed since the last glacial maximum (LGM); with exposure ages younger than the LGM of about 20 kyr, based on ¹⁰Be terrestrial cosmogenic nuclide and OSL surface exposure dating (Lehmann et al., 2019, 2020).

4.2 Sample preparation

The sample preparation followed the method reported previously (King et al., 2016b). The light-exposed outer layer (>2 cm from the surface) was removed using a diamond saw under subdued red-light conditions with constant water flow to avoid frictional heating. The interior part of the sample was gently crushed with a mortar and pestle and sieved to separate the 150–250 µm grain size. The samples were sequentially treated with 10 % HCl and 30 % H₂O₂ to remove carbonate and organic matter respectively. Once dried, the magnetic fractions were removed using a hand magnet. The K-feldspar fraction was separated by density separation (<2.58 gm cm⁻³) using sodium polytungstate. The grains were mounted on stainless steel discs using Silkospray. Small aliquots of 2 mm diameter (containing ∼100 grains) were prepared as these feldspars were highly luminescent.

4.3 Experimental procedure

The TL luminescence measurements were made using a Risø TL/OSL reader (TL/OSL DA-20; Bøtter-Jensen et al., 2010) equipped with a ⁹⁰Sr∕⁹⁰Y irradiation source (∼0.24 Gy s⁻¹) at the University of Lausanne. A heating rate of 1 ${}^{\circ }\mathrm{C}{\mathrm{s}}^{-\mathrm{1}}$ was used, under constant flow of N₂ gas. The TL emission was restricted to violet-blue (395±30 nm) using a filter combination of BG3 and BG39. The measurement details are discussed below. Typically, the minimum detectable limit for the present instrument is ∼300 photon counts per second (cps), considering the signal should be 3 times the background level, which is ∼100 cps. The present highly luminescent feldspar has a maximum photon count of ∼10⁶ cps. This restricts the use of the TL signals up to $\sim {\mathrm{10}}^{-\mathrm{3}}$ % of maximum TL signals.

4.3.1 Measurements

Following Biswas et al. (2018), three sets of experiments were performed to constrain the growth parameters (D₀,a), thermal decay parameters ($E,s,b$) and athermal decay parameters (ρ^′). The athermal frequency factor ($\stackrel{\mathrm{̃}}{s}$) is taken as $\mathrm{3}\times {\mathrm{10}}^{\mathrm{15}}{\mathrm{s}}^{-\mathrm{1}}$ (Huntley, 2006).

The growth parameters and the natural TL level, i.e. the trapped charge population (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$), are estimated using the multiple aliquot regeneration dose (MAR) protocol (Aitken, 1985) with post-glow normalization (Tang and Li, 2017). Eight regeneration doses (0, 24, 47, 118, 236, 472, 944 and 1888 Gy) were given and three aliquots were used for each dose point. A cut heat of 200 ^∘C was applied to remove traps that are unstable over laboratory timescales. We observed a significant sensitivity change (decrease) during the very first measurement of natural TL, which means that natural and regenerative TL signals were not measured under identical TL sensitivity conditions. To circumvent this sensitivity change, we adopted the natural correction factor method (NCF; Chauhan and Singhvi, 2019; Singhvi et al., 2010, 2011). However, the NCF was initially developed for quartz OSL (Singhvi et al., 2011) and should be adapted for feldspar.

The NCF method for quartz relies on the fact that the 110 ^∘C TL peak and the blue stimulated OSL are correlated (Singhvi et al., 2011). In contrast, the TL of feldspar does not exhibit a distinct 110 ^∘C TL peak. The luminescence process in feldspar is more complicated because it arises from a continuous distribution of trapping energies and the dose response characteristics (D₀) vary along the TL glow curve (Fig. 1). To circumvent these issues, we proceeded as follows. We first give a small dose (i.e. <100 Gy) in addition to the natural dose and subsequently measure the TL signal up to 200 ^∘C (TL₁). Then the sample is annealed by heating it to 450 ^∘C, which is followed by a dose of the same amount and measurement of the TL glow curve to 200 ^∘C (TL₂). We observe that the TL sensitivity of the natural measurement is higher than the post-natural-regeneration measurement (Fig. 8a). We then calculate the NCF at different temperatures between 90 and 150 ^∘C, similar to the use of the 110 ^∘C TL peak for quartz. We find that the NCF decreases with increasing temperature during the TL measurement (Fig. 8c). Since there is no direct way to measure the NCF beyond 150 ^∘C, we then extrapolate the NCF value at the region of interest to higher temperatures, i.e. 210–250 ^∘C (Fig. 8c), which we call var-NCF. In turn, the trapped charge population (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) is corrected with the corresponding factor which is between 1 and 2 for samples MBTP1 and MBTP9 (Table 1 and Fig. 8). Finally, we investigated the effects of variable doses and the NCF and found it had no effect for doses below 100 Gy (Fig. 8b).

Figure 8 (a) The lower-temperature TL before (TL₁) and after (TL₂) natural TL measurement (up to 450 ^∘C) of samples MBTP1 and MBTP9. (b) NCF for TL signal (integrated over 90–120 ^∘C) for different given doses. The data are scattered and vary within ±10 %, which possibly suggests the NCF is dose independent. (c) Plot of var-NCF ($={\text{TL}}_{\mathrm{1}}/{\text{TL}}_{\mathrm{2}}$) at different temperature (90–150 ^∘C) along TL glow (solid circles), and extrapolated to calculate NCF in the region of interest (210–250 ^∘C; empty circles). The values were calculated by taking the average of three measurements.

The thermal decay parameters are estimated using the T_m−T_stop method (McKeever, 1980) and analysed by subtraction and fitting of sub-peaks (Pagonis et al., 2014). For the athermal decay parameter, a fading experiment (Huntley and Lian, 2006) was performed for different delay times; aliquots were preheated to 200 ^∘C prior to storage.

4.3.2 Estimating the kinetic parameters

The kinetic parameters of growth (D₀,a), thermal decay (E, s, b) and athermal decay (ρ^′) were inferred using the approach of Biswas et al. (2018) for all thermometers (210–250 ^∘C, 10 ^∘C interval). The results are summarized in Table 1 and shown in S3. It can be noted that with increasing the TL temperature (or thermometer), the activation energy (E) increases and athermal fading (ρ^′) decreases (Table 1). The dose rate ($\dot{D}$) values, another growth parameter, were taken from Lehmann et al. (2020). Since a cut heat of 200 ^∘C was applied for the MAR growth analysis and fading experiments, we focus on the 210–250 ^∘C TL thermometers (i.e. four thermometers). We did not use TL signals beyond 250 ^∘C, as they are insensitive to typical surface temperature fluctuations, as discussed in Sect. 2.2.

Table 1 List of kinetic parameters that describe growth ($\dot{D},D_{\mathrm{0}},a$), thermal decay ($E,s,b$) and athermal decay (ρ^′) and natural TL (observed) trapped charge populations (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) for the four thermometers (210–250 ^∘C, 10 ^∘C interval) for samples MBTP1 and MBTP9. The dose rate ($\dot{D}$) values were taken from Lehmann et al. (2020). Note that in the thermal decay parameters (E,s and b) no errors are mentioned. The mean values of these parameters are calculated from the distribution and an arbitrary error of 5 % was considered.

4.4 Predicting the surface temperature

The measured TL signals (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) are then inverted to infer the thermal history as described in Sects. 3.2 and 3.4 (synthetic approach 2). For the thermal histories, we again use the Greenland ice core δ¹⁸O record (Svensson et al., 2008), which we scaled as described in Sect. 3.2. We assume that the atmospheric temperatures of the Mont Blanc massif followed the trend observed for the Greenland ice core data over the last 60 ka. Note that the temperature increase during the last glacial cycle was synchronous with the temperature anomalies observed in Greenland (e.g. Heiri et al., 2014; Schwander et al., 2000; van Raden et al., 2013). The rationale here is that all temperatures in the Mont Blanc massif follow the Greenland ice core δ¹⁸O data but the amplitude of temperature oscillation (minimum temperature at ∼20 ka to maximum temperature at the present day) and mean temperature are unknown. We pick the amplitude of temperature oscillation randomly between 0 and 40 ^∘C, and the base temperature (temperature at ∼20 ka) between −20 and 30 ^∘C. By generating a large number of random thermal histories (300 000), the probability density function is constructed as discussed in Sect. 3.2. The results of two samples, MBTP1 and MBTP2, are shown in Fig. 9 and suggest that the temperature rose from $-{\mathrm{4.6}}_{-\mathrm{4.1}}^{+\mathrm{3.7}}$ to ${\mathrm{6.2}}_{-\mathrm{3.5}}^{+\mathrm{3.1}}$ ^∘C for sample MBTP1 and $-{\mathrm{2.0}}_{-\mathrm{4.1}}^{+\mathrm{3.9}}$ to ${\mathrm{7.9}}_{-\mathrm{3.1}}^{+\mathrm{3.0}}$ ^∘C for sample MBTP9, since 20 ka, considering 1σ uncertainty. The inferred median suggests an increase of ∼10–11 ^∘C for the rock surface temperature over the last 20 ka.

Figure 9 Inferred rock surface temperature history for sample MBTP1 and MBTP9 collected in the Mont Blanc massif at an altitude of ∼2.0–2.5 km, obtained through inverse modelling of TL data as described in Sect. 4.4. Panels (a–d) are for the sample MBTP1. Panel (a) is the probability density function. The red line, white lines and black lines are the predicted median, 1σ and 2σ confidence intervals. Panel (b) is the plot of observed TL (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) and modelled TL (predicted TL through inverse modelling). Panels (c) and (d) are histograms of temperature at 20 ka and present respectively. Panels (e–h) are the result of same analysis for sample MBTP9.

Figure 10 The impact of initial sensitivity correction for past temperature prediction to sample MBTP1 for three different scenarios: (1) no initial sensitivity correction, i.e. NCF=1; (2) initial sensitivity corrected with fixed NCF at 100 ^∘C ($=\mathrm{1.64}\pm \mathrm{0.08}$) for all thermometers; and (3) initial sensitivity corrected to all TL thermometers with var-NCF for all TL thermometers (i.e. the selected method). Panels (a), (e) and (i) are the probability distributions, and (b), (f) and (j) are observed and predicted TL plots for all three scenarios respectively. Panels (c), (g), (k) (d), (h) and (l) are the histograms of predicted temperature for all accepted paths at 20 ka and the present day for all three scenarios respectively.

Figure 11 Result of synthetic test for annual fluctuation of temperature. Panel (a) shows a probability density plot. The dotted red line is the actual thermal history. The solid red line is the predicted median of the isotherms. White and black lines are the 1σ and 2σ confidence intervals. Panel (b) shows the plot of observed (${\stackrel{\mathrm{‾}}{n}}_{\text{obs}}$) and modelled TL values (predicted TL through inverse modelling) for all TL thermometers.

5 Discussion

The theoretical model for the rate equation of trapped charge population in feldspar has been described in several ways; first-order kinetics (Brown et al., 2017; Yukihara et al., 2018), general-order kinetics (Biswas et al., 2018; Guralnik et al., 2015b), charge transport through sub-conduction band-tail states (King et al., 2016a; Li and Li, 2013), Gaussian distribution of trapped energies (Lambert et al., 2020) or localized recombination in randomly distributed defects (Jain et al., 2012). What is common to all these models is that luminescence of feldspar is complicated and exhibits a non-linear non-first-order kind of behaviour due to either the presence of sub-conduction band-tail states (Morthekai et al., 2019; Poolton et al., 2002) or a complex charge transport mechanism. TL in feldspar is even more complicated because it shows continuous distribution of trapping energies (Biswas et al., 2018; Duller, 1997; Grün and Packman, 1994; Pagonis et al., 2014; Strickertsson, 1985) and TL is a more diffusive process than OSL; OSL of feldspar has resonant energy levels (Hütt et al., 1988). Different models were reviewed and tested by Guralnik et al. (2015b) who suggested that the general-order kinetics, a mathematically simplified model, could be used to explain the luminescence phenomenon well. We adopted this model for the TL of feldspar where the power terms (a, b) accounts for the nonlinearity involved in the TL of feldspar. The efficacy of using general-order kinetics has been demonstrated on samples with known thermal history (KTB borehole samples) for OSL of feldspar (Guralnik et al., 2015a) and the TL of feldspar (Biswas et al., 2018).

Here we investigate the difference in temperature sensitivity of different TL thermometers, which correspond to individual TL temperature or TL signals. On the basis of the kinetic parameters derived for our sample, and our sensitivity tests (Sect. 2.2.1), we recommend using TL thermometers with a temperature range of 200 to 250 ^∘C for a typical surface temperature fluctuation, e.g. ∼10 ^∘C. If the temperature fluctuations are larger, higher-temperature TL (>250 ^∘C TL) can be used. The multiple TL signals (200 to 250 ^∘C, 10 ^∘C interval) can constrain thermal history of ∼50 kyr. A higher-temperature fluctuation can be better constrained with a greater number of thermometers (as discussed in Sect. 2.2.1).

For periodic oscillations, when the period is comparable to the lifetime of the trapped electron for a given thermometer, it may be used to infer temporal variation of surface temperature (see Sect. 2.2.2). Typically, tens of thousands of years of temperature oscillation can be detected using TL thermometers with peak temperatures higher than 200 ^∘C (210 to 250 ^∘C). Periodic oscillation with shorter periods (<1 kyr) will exhibit a similar effect to isothermal temperature conditions, yielding a temperature higher than the mean of oscillation.

One outstanding issue when using TL is the sensitivity change during the very first measurements up to 450 ^∘C, which cannot be corrected by post-glow normalization (Tang and Li, 2017). Here we show that sensitivity changes during natural measurements can be monitored for lower-temperature TL (<150 ^∘C) following the same method adopted for the OSL of quartz, which is called the natural correction factor (NCF; Singhvi et al., 2011). Because there is no direct method to track the TL sensitivity change in the region of interest (210–250 ^∘C TL), we simply extrapolate the sensitivity change observed in the lower temperature TL peaks (i.e. 90–150 ^∘C) to the region of interest (i.e. 210–250 ^∘C). This is new and it will need further investigation. However, we find that the effect of the initial sensitivity change on the amplitude of the inferred temperature histories is small. In Fig. 10, we compare inversion results for three different scenarios for sample MBTP1: (1) there is no initial sensitivity correction, i.e. NCF = 1; (2) the initial sensitivity correction is done using the value obtained at 100 ^∘C (NCF${}_{\mathrm{100}}=\mathrm{1.64}\pm \mathrm{0.08}$); and (3) the var-NCF approach described in the Sect. 4.3.1 is used. Although the results show that the sensitivity correction has a significant impact on the absolute inferred temperature, we do not observe much difference between using a constant value and the extrapolated value. Furthermore, and more importantly, the difference between the present-day temperature and temperature at about 20 ka remains about 10–11 ^∘C in the three tested cases (median of the prediction).

The estimated constant erosion rates in these two sample locations, MBTP1 and MBTP9, are $\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{2}}$ mm yr⁻¹ with maximum possible times of erosion of 20.9 and 19.5 ka respectively (Lehman et al., 2020). This translates to maximum erosion depths in these two locations of 0.07 and 0.62 m respectively. At those depths, mean temperature should be in equilibrium with atmospheric temperature (e.g. Hasler et al., 2011). Based on the rock temperature measurement of borehole samples in Mont Blanc massif, Magnin et al. (2017) suggest that up to a 2 m depth the rock temperature is nearly constant with depth or a maximum variation of up to 1 ^∘C (Fig. 6c of Magnin et al., 2017). No effect of erosion on surface temperature is considered here.

For the inverse modelling of natural samples, δ¹⁸O data are used as a prior on the shape of the thermal histories, but we leave two scaling parameters free – minimum temperature at 20 ka, and amplitude (temperature difference between at 20 ka and present) – and we did not include the role of ice on setting the rock temperature during glaciation. Lehman et al. (2020) provided a range of solutions for deglaciation in the present location: either thinning of 450 m of the glacier occurred progressively between ∼17 and ∼12 ka or it was instantaneous. In the absence of a clear scenario, we took two samples from two extreme altitudes of 2545 m a.s.l. (MBTP1) and 2133 m a.s.l. (MBTP9). The top-most sample, MBTP1, had a very thin or no glacier coverage during LGM, and the bottom-most sample, MBTP9, was exposed or covered by ice. Thus, it is expected that the during LGM, the rock temperature of MBTP1 would have been in equilibrium with atmospheric temperature (Hoelzle et al., 1999), whereas it is less clear for MBTP9. If it was covered by ice, it was likely temperate ice and the basal temperature would be close to 0 ^∘C. Interestingly, we find similar results from the inversion; predicted rock surface temperatures of MBTP1 and MBTP9, during LGM, are $-{\mathrm{4.6}}_{-\mathrm{4.1}}^{+\mathrm{3.7}}$ and $-{\mathrm{2.0}}_{-\mathrm{4.1}}^{+\mathrm{3.9}}$ respectively, and the temperature of the bottom-most sample during LGM is close to 0 ^∘C, at least within error.

The application of the introduced method predicts that the final rock surface temperatures at the locations of Mont Blanc massif are 6.2±3.2 ^∘C for MBTP1 (2545 m a.s.l.) and 7.9±3.0 ^∘C for sample MBTP9 (2133 m a.s.l.). Although these temperatures have large uncertainty, they are higher than the mean annual atmospheric temperature in this location. The mean annual temperature of Chamonix (1035 m), a nearby city, is 7.3 ^∘C. Considering an adiabatic lapse rate of 5 ${}^{\circ }\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$, the expected mean annual atmospheric temperatures at the sample locations of MBTP1 and MBTP9 are ∼0 and 2 ^∘C respectively. The offset between the predicted and expected temperature (∼6 ^∘C) can be explained by two main factors: (1) the rock surface temperature is always higher than atmospheric temperature and the temperature difference can be up to 10 ^∘C (Magnin et al., 2019), and (2) seasonal temperature fluctuations may lead to an overestimation of the mean annual temperature (as discussed in Sect. 2.2.2). To quantify this latter offset, we performed a simple synthetic test, with annual oscillation of +10 ^∘C (summer) to −10 ^∘C (winter) with mean at 0 ^∘C, up to 20 ka (before that temperature was set to a 0 ^∘C isotherm), and predicted the equivalent isothermal temperature using the inverse approach. The result suggests a mean annual temperature that is 2.7±0.7 ^∘C higher than the mean temperature of the periodic signal (Fig. 11), confirming the results of Guralnik and Sohbati (2019).

The inverse modelling results show an increase in rock surface temperature in the Mont Blanc Massif of ∼10–11 ^∘C (considering median of the prediction) from 20 ka to today. The median of the distribution of possible thermal histories of the two samples follows Greenland ice core δ¹⁸O anomalies with missing low frequencies. Climate reconstructions in Europe using fossil pollen suggest that the mean annual temperature anomaly (the difference between the temperature at the LGM and today) is 12±3 ^∘C in the north of Pyrenees–Alps line (Peyron et al., 1998). Wu et al. (2007) inferred that LGM temperatures in Europe were ∼10–15 ^∘C lower than the present-day temperature based on pollen analysis. Although there are large uncertainties associated with pollen data and its methodological constraints, the overlap in temperature estimates between the two proxies suggests that TL may be a reliable paleothermometer.

6 Conclusions

A new approach to reconstruct the temporal variation of rock surface temperature using the TL of feldspar is introduced. Forward modelling of different TL signals suggests that TL signals in the range of 210 to 250 ^∘C are sensitive to typical surface temperature fluctuations, which we define as TL thermometers. Multiple TL thermometers (210–250 ^∘C, 10 ^∘C interval) can then be used to constrain thermal histories of rocks over ∼50 kyr for temperature fluctuations of ∼10 ^∘C. The sensitivity of the periodic forcing on trapped charge populations suggest that natural TL is sensitive enough to mean temperature and amplitude of periodic forcings. Typically, tens of thousands of years of temperature oscillation can be predicted using this approach. Finally, we show that it is possible to recover thermal histories of rocks when one assumes that the temperature followed observed the Greenland ice core δ¹⁸O record.