## Abstract

Thermochronology can provide information on palaeotopography and its evolution. We present a new method to reconstruct the shape of the palaeotopography from the palaeoisotherm derived from low-temperature thermochronological data. While deriving palaeoisotherm and reconstructing palaeotopography, the exhumation rate may also be constrained. The proposed method is independent of the relationship between the ancient and modern topography; however, it requires that the palaeotopography maintains its shape and the subsurface thermal field is invariant with time during the period when a set of samples passed through the closure temperature. It is shown that the inherent uncertainties in sample age and the limited sampling density will inevitably induce errors in the reconstructed topography, and these errors can be reduced by eliminating the noise in the derived isotherm. Increasing the number of samples helps to reduce the noise, and if, before sampling, we can make rough estimates about the geothermal gradient, the exhumation rate, the uncertainties in sample age, and the wave components and the relief of the palaeotopography to be reconstructed, the number of samples can be tentatively decided. It is also shown that the reconstructed topography is sensitive to exhumation rate and geothermal gradient, and a concrete sensitivity analysis is needed for a given dataset. By reinterpreting the (U–Th)/He data from the Sierra Nevada, California, we show that our method can reveal palaeorelief as well as other valuable information on the palaeotopography of the studied area.

The Earth's surface and its evolution are among the most important factors that affect the thermal structure of the uppermost crust, to which low-temperature thermochronological systems, such as (U–Th)/He or fission-track dating in apatite and zircon, are most sensitive (e.g. Turcotte & Schubert 1982; Stüwe *et al.* 1994; Mancktelow & Grasemann 1997; Braun 2002*a*, *b*). Although this complicates the interpretation of low-temperature thermochronological data, it enables some of the important attributes of the palaeotopography to be constrained. These attributes may include topographic relief or the form, location and scale of topographic relief in the past, but do not include palaeoelevation because thermochronology is only concerned with vertical movement of rocks towards the Earth's surface, which is a moving boundary relative to sea level, etc. (Reiners 2007).

In the last decade, there have been several studies to constrain palaeotopography using data from low-temperature thermochronology (e.g. House *et al.* 1998, 2001; Braun 2002*a*, *b*; Persano *et al.* 2002; Ehlers *et al.* 2003; Reiners *et al.* 2003; Braun & Roberts 2005). Some of these studies are based on numerical modelling and require a topography evolution model that links the palaeotopography with the present-day topography. For example, the model by Braun (2002*b*) hypothesized that the ancient and the current topography had the same shape but different relief (or unevenness). However, so far no model can offer a detailed description of the long-term topography evolution because many factors that may influence the evolution of topography are difficult to constrain.

Here we present a new method to reconstruct the shape of palaeotopography from the thermochronological data. This method is independent of the relationship between the ancient and the present-day topography, but requires that the palaeotopography maintains its shape and the subsurface thermal field is invariant with time during a period in the past; for example, during the period when a set of samples are passing through the closure temperature.

## Reconstructing palaeotopography from palaeoisotherm

Surface topography can be expanded in an infinite Fourier series. In practice, the sum of a finite number of cosine and sine is adequate to approximate surface topography:
1
where *h*(*x*) is the height of ground surface at the horizontal position *x*, λ is the length of a given transect, and *a*_{n} and *b*_{n} are amplitudes of cosine and sine functions whose wavelength are λ/*n*. Equation (1) can be expressed in another form:
2
where *c*_{n}=√[(*a*_{n})^{2} + (*b*_{n})^{2}], *c*_{n} and ϕ_{n} are amplitude and phase, respectively, of the *n*-order harmonic whose wavelength is λ/*n*.

Because heat conduction in solids can be approximated as a linear physical process, the total perturbation of the thermal field caused by the topography can be adequately approximated by the sum of the perturbations caused by each of the components of the Fourier series. The perturbation caused by a periodic topography with small amplitude has been shown to be proportional to the amplitude of the topography, and to decay exponentially with depth (Turcotte *et al.* 1982). The steady-state solution for topography of constant geometry can therefore be expressed as (Mancktelow & Grasemann 1997):
3
4
where *u* is exhumation rate, κ is thermal diffusivity and ϕ is the lapse rate of surface temperature with altitude (−4.5 °C km^{−1} altitude). *T*_{0}(*y*) describes the geotherm in the absence of any surface perturbations, and can be fixed according to the temperature at the surface and at the lower boundary of a slab (see equation 23 of Mancktelow & Grasemann 1997).

During the period when a set of samples are passing through their closure temperature (*T*_{c}) – if the palaeotopography maintains its shape and the subsurface thermal field is invariant with time – the isotherm of *T*_{c} maintains its shape and depth. Assuming the exhumation rate is constant, the difference in the depth of palaeoisotherm of *T*_{c} at the horizontal position of two samples, Sample_{1} (Age_{1}, Elevation_{1}) and Sample_{2} (Age_{2}, Elevation_{2}), is:
5

The shape of palaeoisotherm is determined by equation (5) and the average depth of palaeoisotherm (*y*) can be estimated through the following equation:
6

Thus, the depth of palaeoisotherm at the horizontal position of each sample can be obtained, that is for each sample we have:
7
where *x*_{i} and *y*_{i} are the horizontal position and the depth of palaeoisotherm at the position of the *i*-th sample, respectively. Expanding equation (7), a system of linear equations with unknown values of *a*_{n} and *b*_{n} (0≤*n*≤*N*) as variables is obtained. These unknown values can be fitted by multiple linear regression. Then *a*_{n} and *b*_{n} (0≤*n*≤*N*) can be used in equation (1) and the palaeotopography is reconstructed.

A vital issue is how to select a suitable value of *N*. Obviously, the greater the value of *N*, the more detailed the topography to be reconstructed, so we have called λ/*N the precision of topography reconstruction* (PTR). According to the theory of linear regression, the maximal value of *N* cannot exceed half of the number of samples. In reality, the value of *N* should be smaller than its maximum owing to data uncertainties.

## Uncertainties in the calculated palaeoisotherm and selection of the PTR for a given dataset

### A synthetic example

Here, a synthetic example is used to test the validity of the method and to illustrate how uncertainties in the calculated isotherm affect the reconstruction of palaeotopography. In our experiment, we use a 192 km-long topographic profile from a 1 km-resolution DEM (GTOPO30) of Dabie Shan, eastern China.

We first calculated the steady-state thermal field using the finite-difference method and accepted it as the ‘real’ thermal field under the topography. In this forward model, the parameters were set as: the surface temperature *T*_{s}=0 °C, the thickness of the layer being exhumed *L*=100 km, the constant exhumation rate *u*=50 m Ma^{−1}, the heat diffusivity κ=10^{−6} m^{2} s^{−1} and the lapse rate of surface temperature with altitude ϕ=−4.5 °C km^{−1}. Because heat production is not considered in this model, the temperature at the low boundary is set as *T*_{L}=2500 °C to ensure a near-surface geothermal gradient of approximately 25 °C km^{−1}. This basal temperature is clearly much too large, but it is the selection of the near-surface geothermal gradient, instead of the basal temperature, that affects the result.

The isotherms obtained by the finite-difference method are then sampled at equal intervals to form a dataset (DS1) of 128 data pairs. Each data pair contains the horizontal position and the depth of isotherm at that position. Finally, the topography is reconstructed with this dataset using the method that we have described above, with the value of the PTR set as 6 km.

As shown in Figure 1, taking the 75 °C isotherm as an example, the reconstructed topography is almost identical to the original topography. Therefore, if the isotherm could be accurately derived from the low-temperature thermochronological data, the topography can be accurately reconstructed given that the parameters, such as the exhumation rate and the near-surface geothermal gradient, are known.

However, uncertainties are inherent in analysing the cooling age of rocks. They will lead to uncertainties in the calculated isotherm and errors in the reconstructed topography. According to equation (5), the error of the calculated isotherm at the horizontal position of the *i*-th sample is:
8
where ε_{ai} is the error in the *i*-th sample age. If the error in all sample ages obeys the normal distribution with a mean of zero and one standard error (1 SE) of σ, then ε_{i} obeys the normal distribution with a mean of zero and a deviation of *u*σ.

To illustrate how uncertainties in the calculated isotherm affect the reconstruction of palaeotopography, random error is added to the depth of isotherm of every data pair of DS1 to form a new dataset (DS2). The random errors are generated from a normal distribution with a mean of zero and 1 SE of 150 m to simulate the situation in which 1 SE of sample age is 3 Ma and the exhumation rate is 50 m Ma^{−1}.

Using DS2, the wave components whose wavelengths exceed 8, 32, 48, 96 and 192 km (i.e. λ=192 km, and *n*=24, 6, 3, 2 and 1, respectively, see equation 2) are regressed respectively, and the reconstructed topographies are shown in Figure 2. In order to quantify the ‘goodness’ of topography reconstruction, average deviations between the reconstructed and the original topography are calculated at every 1.5 km.

As shown in Figure 2, the topography reconstructed with long wave components (Fig. 2a) exhibits a similar overall trend as the original topography, but lacks detail. More details are revealed with decreasing wavelength (Fig. 2b). However, with a further decrease in wavelength, drastic oscillations of short wavelength appear in the reconstructed topography (Fig. 2c). Correspondingly, with a decrease in wavelength, the deviation reduces (Fig. 2b), but this deviation gradually widens with a further decrease in wavelength.

Obviously, owing to uncertainties in sample ages and errors in calculated isotherm, only the long wave component of the topography can be generally reconstructed. However, the problem is how to select an appropriate wavelength and to reconstruct the topography given the dataset and the uncertainties in sample ages.

### The PTR for a given dataset

To understand how uncertainties in the calculated isotherm affect the reconstruction of palaeotopography, isotherms obtained by the above two datasets (DS1 and DS2) are broken down into different wave components by Fast Fourier Transform (FFT) procedure. The amplitude of each wave components is shown in Figure 3. Where *c*_{tn} (the actual amplitude) is for DS1 and *c*_{n} (the calculated amplitude with errors) is for DS2.

As shown in Figure 3, for wave components whose wavelength are longer than 27.4 km (i.e. *n*≤6), the calculated amplitudes are near the actual amplitudes and so the corresponding wave components in the topography are well reconstructed, as shown in Figure 2a and b. For the wave components whose wavelength are shorter than 32 km (i.e. *n* > 6), the actual amplitudes are near zero; however, the corresponding calculated amplitudes are much greater than zero and thus lead to short-wavelength oscillations in the reconstructed topography, as shown in Figure 2c.

In fact, if the actual amplitude is zero, the uncertainties of sample ages may lead to calculated amplitudes much greater than zero. It can be proven that, under this circumstance:
9
where χ^{2}(2) is the Chi-square distribution with a mean of 2. Therefore, if a calculated *c*_{n} is small, they may be noise induced completely by uncertainties in the sample ages and limited sampling density. However, it is very difficult to tell whether a value of *c*_{n} is small enough to be regarded as noise and should be eliminated from the calculated isotherm while reconstructing palaeotopography.

Here, based on our model experiments and considering the statistic characteristics of the calculated amplitudes, we select a threshold value of: 10

This threshold value (TV) is used to help us discern whether a wave component of the isotherm is noise or not. In this example the threshold value is 0.046 km, and the amplitudes for wave components of 192, 96, 64, 38.5 and 32 km exceed the threshold value (Fig. 3). Therefore, these wave components of the isotherm can be discerned from the noise, and the corresponding wave components of the topography should be reconstructed.

In practice, we reconstruct all the wave components of the topography whose wavelength exceeds that of the shortest discernible wave component of the isotherm. In this example the calculated amplitude of the 48 km-wave component of the isotherm is less than the threshold value. Nonetheless, for the sake of convenience, the 48 km-wave component of the topography is reconstructed.

Therefore, in our paper, the PTR is defined as the wavelength of the shortest discernible wave component of the isotherm and all wave components whose wavelengths exceed the PTR are reconstructed. In this example, the PTR is 32 km. It should be mentioned that:

the threshold value is a statistical concept, so the noise may exceed it; in fact, other threshold values may be selected and it is suggested its value should not be larger than 5

*u*σ/√*N*_{s};in most of our model experiments the deviation between the original topography and the topography reconstructed using the PTR thus selected is the least.

## Some factors that affect the isotherm and the selection of sampling density

Whether a wave component of the palaeotopography can be reconstructed is dependent on two facts:

the wave component can cause remarkable perturbation on the isotherm;

this perturbation can be discerned from the noise given the number of the samples and the uncertainties in the sample ages.

These two factors can be linked by the following inequality: 11

That is to say, if the wave component of the topography could cause an undulation whose amplitude exceeds 3.2*u*σ/√*N*_{s} in the isotherm, the corresponding calculated amplitude would have a good chance (50%) of exceeding the threshold value of 3.2*u*σ/√*N*_{s} and to be discerned from the noise. Then what factors affect the magnitude of the perturbation on the isotherm induced by topography?

It has been shown that the perturbation is proportional to the amplitude of the topography. As shown in Figure 4, for a periodic topography, the isotherm follows the shape of the topography but with different amplitude, and the ratio α (Braun 2002*b*) between the amplitude of the topography and that of the isotherm is almost constant, only weakly sensitive to topographic amplitude (Reiners 2007).

However, the amplitude of the topography is just one of the factors that affect the magnitude of the perturbation. Other factors include the geothermal gradient, the wavelength of the topography and the temperature of the isotherm. Figure 5 shows the ratios (α) for the topographies of different wavelengths and the isotherms of different temperatures when the geothermal gradient is 10, 25 or 40 °C km^{−1}, where the heat diffusivity is 10^{−6} m^{2} s^{−1} and heat production is not considered. It is shown in Figure 5 that, as the wavelength of the topography and the near-surface geothermal gradient decrease, the ratio α decreases. In contrast, the ratio α increases as the temperature of the isotherm decreases.

For topographies of the same amplitude, a smaller α means that the perturbation, induced by the topography, is more difficult to discern from the noise in the calculated isotherm. Nonetheless, the corresponding topography could be reconstructed, if we were able to reduce the noise by increasing the number of samples.

However, it is difficult to decide how many samples are sufficient for our purpose, because such a decision relies on our knowledge of the geothermal gradient, the exhumation rate, and, in particular, the wave components and the relief of the topography to be reconstructed. If we do have some vague knowledge about these factors, then Figure 5 can help us tentatively to decide the number of the samples. In the case shown in Figure 5a, assuming the main wave component of the palaeotopography is 40 km, its amplitude is 1 km and the temperature of the isotherm is 105 °C, then the amplitude of isotherm is approximately 0.13 km. If the 1 SE of the sample age is 3 Ma, then, according to inequality (11), Therefore, 14 samples are needed on a 40 km-long transect. It should be mentioned that inequality (11) is only applicable in the situation where the number of the samples is a multiple of 2 and the samples are sampled at equal intervals.

In practice, the topography contains a series of wave components, and numerical modelling is needed to calculate the ‘conjectural’ isotherm and break it into different wave components. Then the number of samples can be decided according to the amplitudes of the wave components in the isotherm. In the case of a moderate exhumation rate of 50 m Ma^{−1} and a geothermal gradient of 25 °C km^{−1} and one standard error of 3 Ma or so in apatite (U–Th)/He ages, and assuming the topography in Figure 1 is a rough estimation of the topography to be reconstructed, the wave components of the isotherm are shown in Figure 3. According to the amplitude of each component and inequality (11), four samples are needed for the reconstruction of the 192 km-wave component of the topography, 32 samples are needed for the reconstruction of the 96 km-wave component, 64 samples are needed for the reconstruction of the 38.5, 48 and 64 km-wave components, and 128 samples needed to reconstruct the 32 km-wave component.

## Sensitivity to exhumation rate and geothermal gradient

For a given dataset, the sensitivity of reconstructed topography to the exhumation rate is related to the age and height distribution of the samples. According to equation (5), if the samples come from locations of the same height, the isotherm is determined completely by the age distribution of the samples. In such case, the relief of the isotherm and then the reconstructed topography are proportional to the exhumation rate. However, if all the sample ages are identical, the isotherm is determined completely by the height distribution of the samples, and the exhumation rate has no influence on the reconstructed topography.

In practice, a concrete analysis is needed. As shown in Figure 6a, different exhumation rates are used to reconstruct the topography in Figure 1 with DS2, and the reconstructed topography is not sensitive to exhumation rate.

The near-surface geothermal gradient is another parameter that has a vital influence on the reconstructed topography. According to Figure 5, this influence varies with the wavelength of the topography to be reconstructed. Supposing the exhumation rate is 50 m Ma^{−1} and the actual geothermal gradient is 40 °C km^{−1}, an observed amplitude of 0.5 km for the 30 km-wave component of the 100 °C isotherm means a relief of 1 km or so in the topography (see Fig. 4c). However, if we reconstruct the topography with a geothermal gradient of 10 °C km^{−1}, we will get a relief of 8 km or so in the topography (see Fig. 4a).

In practice, limited samples are only sufficient to reconstruct the long wavelength of the topography. Under this circumstance, the reconstructed topography is not as sensitive to geothermal gradient. As shown in Figure 6b, different geothermal gradients were used to reconstruct the topography in Figure 1 with DS2, where the ‘real’ geothermal gradient is 25 °C km^{−1}.

## Constraining exhumation rate while reconstructing palaeotopography

Exhumation rate can be constrained by the slope of the ‘age–elevation relationship’. This method requires that all the samples have undergone the closure temperature at constant depth below a mean surface. If the samples are not come from a vertical profile, the influence of topography on the isotherm should be taken into account, and a topographic correction is needed (Stüwe *et al.* 1994). However, the palaeotopography is unknown. Here, we propose a new method to constrain the exhumation rate while reconstructing palaeotopography.

According to our model experiments, if the samples do not come from locations of the same height, an exhumation rate other than the real value is *inclined* to boost up the short wave components much more than that of long wave components of the isotherm. This will lead to a ‘rougher’ calculated isotherm than the real one, as shown in Figure 7, where the actual exhumation rate is 50 m Ma^{−1}.

Therefore, a smoother isotherm is usually preferred. For a given dataset and different exhumation rates (e.g. *u*_{1}, *u*_{2}), the corresponding PTRs (e.g. PTR_{1}, PTR_{2}) can be determined by the threshold value. After eliminating from the calculated isotherms the wave components whose wavelengths are shorter than the PTR, using equations (5) and (6), smoothed isotherms are obtained.

If PTR

_{1}> PTR_{2}, the smoothed isotherm with*u*_{1}contains less short wave components and is therefore smoother, and so*u*_{1}is preferred.If PTR

_{1}=PTR_{2}=PTR, reconstructing the topography with*u*_{1}and*u*_{2}, respectively, and calculating the corresponding root mean square (RMS) of the difference between the predicted and the observed ages:

where Age_{pi} and Age_{oi} are the predicted and observed ages for the *i*-th sample, respectively. Note that RMS is the ‘distance’ between the calculate isotherm and the smoothed isotherm (not the real isotherm); a small RMS means that the short wave components (<PTR) in the calculated isotherm are less remarkable and the calculated isotherm is smoother, so the corresponding exhumation rate is more preferred.

Here we provide a synthetic example. In this example, the sampling interval is 12 km, exhumation rate is 50 m Ma^{−1}, and 1 SE of sample age is 3 Ma. Through a forward modelling, a synthetic dataset was obtained, as shown in Figure 8a. The isotherm is then calculated using different exhumation rates between 10 and 100 m Ma^{−1}, and the PTR was determined. The PTR is 192 km for exhumation rates of between 40 and 100 m Ma^{−1}, and is shorter than 192 km for exhumation rates of between 10 and 30 m Ma^{−1}, so an exhumation rate larger than 30 m Ma^{−1} is preferred. For different exhumation rates from 40 to 100 m Ma^{−1}, topographies are reconstructed and RMSs are calculated. As shown in Figure 8b, the RMS is least when the exhumation rate is 50 m Ma^{−1}.

Because the sampling interval is too large and the height difference between samples is too small, it is not suitable to use the slope of the ‘age–elevation relationship’ to constrain the exhumation rate with this synthetic dataset, as shown in Figure 8c. In fact, under such circumstances, the slope of the ‘age–elevation relationship’ tends to underestimate the exhumation rate because of the influence of the topography.

## The example of the Sierra Nevada, California

One of the most successful attempts to directly constrain the palaeotopography by low-temperature thermochronology was undertaken by House *et al.* (1998, 2001) in the Sierra Nevada of northern California. This region is a relatively simple west-tilted block with no evidence of significant folding or faulting. These features make it an ideal place to study the evolution of topography.

In the study by House *et al.* (1998), 36 samples were collected at the approximately same height of 2 km along a 200 km range-parallel transect. As shown in Figure 9a, the average sampling interval is about 6 km, and most of the mean apatite (U–Th)/He ages range from 50 to 80 Ma. The 1 SE of sample age is approximately 3 Ma. House *et al.* (1998) showed that if the geothermal gradient is close to 20 °C km^{−1} and radioactive heat production is considered for (U–Th)/He samples taken from the same elevation, the age distribution is mainly determined by the relief and less affected by the wavelength of the topography. According to their calculations, the observed 20–30 Ma variation in (U–Th)/He ages correspond to a long-wavelength relief of 2–4 km.

We reinterpreted the data of House *et al.* (1998) using our method. The parameters were set as: *L*=100 km, *T*_{s}=0 °C, *T*_{L}=2500 °C, κ=10^{−6} m^{2} s^{−1}, ϕ=−4.5 °C. Therefore, the near-surface geothermal gradient is 27 °C km^{−1}, a little larger than that suggested by House *et al.* (1997). We assumed that the closure temperature of apatite (U–Th)/He thermochronological system is 75 °C (e.g. Wolf *et al.* 1996).

The exhumation rate in this area is between 40 and 100 m Ma^{−1} (House *et al.* 1998; Braun 2002*a*). We first selected an exhumation rate of 40 m Ma^{−1} and calculated the 75 °C isotherm, as shown in Figure 9b. Under this circumstance, the threshold value is 0.074 km, so the PTR is 38.5 km. The same value of the PTR is also obtained when the exhumation rate is set anywhere between 10 and 100 m Ma^{−1}.

We tried to constrain the exhumation rate while reconstructing the palaeotopography. As shown in Figure 9c, the exhumation rate should not be less than 40 m Ma^{−1}. However, the RMS is not affected significantly by the change in exhumation rate when the rate is greater than 40 m Ma^{−1}. Therefore, the exhumation rate cannot be well constrained with this dataset.

The topography at 50–80 Ma was reconstructed for different exhumation rates, as shown in Figure 9d. If the exhumation rate is 40 m Ma^{−1}, the palaeorelief would be 2–3 km, and an exhumation rate of 100 m Ma^{−1} would give a relief of 4–6 km. Nevertheless, the reconstructed topography has a larger relief compared to the present-day topography. This suggests that the current topography is probably the result of long-term erosion of the mountain belts. In addition to palaeorelief, other valuable information on the palaeotopography of the studied area can be revealed. For example, by comparing the current topography and the reconstructed topography, we note that the ancient mountain ridge between valley S and K is narrower than its current state. This may be related to a lateral migration of the drainage system in this area after 50 Ma.

## Conclusions

In a tectonically stable region, assuming the palaeotopography maintains its shape and the thermal field is invariant with time during the period when a set of samples is passing through the closure temperature, the shape of palaeotopography can be reconstructed from the isotherm derived from the low-temperature thermochronological data. By reinterpreting the (U–Th)/He data from the Sierra Nevada, California, we show that this method can reveal more details, in addition to the relief (i.e. the unevenness), of the palaeotopography.

Because of the inherent uncertainties in sample age and the limited sampling density, errors are inevitable in the calculated isotherm. If the amplitudes of the wave components of the calculated isotherm are small, they may be noise induced completely by uncertainties in the sample ages and should be eliminated from the calculated isotherm while reconstructing palaeotopography. A threshold value, based on model tests and the statistic characteristics of the calculated amplitudes, is needed to help us discern whether a wave component of the isotherm is noise or not.

It is difficult to decide how many samples are sufficient for our purpose because such a decision relies on our knowledge of the geothermal gradient, the exhumation rate, the uncertainties in sample ages, and, in particular, the wave components and the relief of the topography to be reconstructed. However, if we do have some knowledge about these factors, then the number of the samples can be tentatively decided.

The reconstructed topography is sensitive to exhumation rate and geothermal gradient. The sensitivity varies with many factors and cannot be generalized. Therefore, a concrete sensitivity analysis is needed for a given dataset.

The exhumation rate can be constrained while deriving palaeoisotherm and reconstructing palaeotopography. This method can be used in the situation where the slope of the ‘age–elevation relationship’ is not applicable to constrain the exhumation rate when the sampling interval is too large and the height difference between samples is too small.

Owing to the presence of many factors that are hard to constrain, it is difficult to establish models of long-term topography evolution. Our method is independent of the relationship between the ancient and the present-day topography, and can be used to constrain models of long-term topography evolution.

## Acknowledgments

This work has been supported by Chinese National Science Foundation (Programs 40572075 and 40621063). We are grateful to J. Braun and C. Persano for their constructive reviews that greatly improved the manuscript.

## Appendix

### Appendix: A statistic analysis for the uncertainties in the calculated isotherm

The depth of the isotherm in the horizontal position of the *i*-th sample can be denoted as:
A1
where *y*_{i} is the calculated depth from equations (5) and (6), and *y*_{ti} and ε_{i} are the true value and the error of *y*_{i}, respectively. According to equation (5), the error of the isotherm in the position of the *i*-th sample is:
A2
where *u* is the exhumation rate and ε_{ai} is the error in the *i*-th sample age. Supposing the error in all sample ages obeys the normal distribution with a mean of zero and a deviation of σ^{2}, and then ε_{i} obeys the normal distribution with a mean of zero and a deviation of *u*^{2}σ^{2}.

If the number of the samples is a multiple of 2 and the samples are sampled at equal intervals, then the calculated isotherm can be broken down into a series of cosine and sine functions (see equation 1) by Fast Fourier Transform. The amplitudes of the cosine and sine functions, whose wavelength are λ/*n*, consist of two elements:
A3
A4
where *a*_{tn} and *b*_{tn} are Fourier transformations of *y*_{ti}, *a*_{en} and *b*_{en} are Fourier transformations of ε_{i}, and *N*_{s} is the number of samples. The observed value and the true value of the amplitude of the *n*-order wave component of the isotherm are:
A5
A6
If errors in different sample ages are independent of each other, it can be proven that *a*_{en} and *b*_{en} obey the normal distribution, and:
A7
A8

Supposing *a*_{tn} and *b*_{tn} are zero, then *C*_{n}^{2}=*a*_{en}^{2} + *b*_{en}^{2}. Under this circumstance, the errors in the calculated isotherm result completely from the uncertainties of the sample ages. Because *a*_{en} and *b*_{en} obey the normal distribution and are independent of each other, then:
A9
where χ^{2}(2) is Chi-square distribution with a mean of 2. The conditional probability of *P*(*C*_{n} < *C*|*C*_{tn}=0) can be readily obtained, as shown in Figure 10a. According to Figure 10a, *P*(*C*_{n}>*u*σ/√*N*_{s}|*C*_{tn}=0)=1−*P*(*C*_{n} < *u*σ/√*N*_{s}|*C*_{tn}=0)=0.8. This is to say, if the actual amplitude is zero, then the calculated amplitude would exceed *u*σ/√*N*_{s}, with a probability of 80%. Therefore, if a calculated *c*_{n} is small (e.g. smaller than *u*σ/√*N*_{s}) it may be induced completely by uncertainties in the sample ages and should be regarded as noise and eliminated from the calculated isotherm while reconstructing the palaeotopography. However, it is not easy to tell whether a value of *c*_{n} is small enough to be regarded as noise.

Based on our model experiments, and considering *P*(*C*_{n}<3.5*u*σ/√*N*_{s}|*C*_{tn}=0)=95%, we selected a threshold value of 3.5*u*σ/√*N*_{s} to help us to discern whether a wave component of the isotherm is noise or not. It should be mentioned that the noise may exceed this threshold value and this might cause intolerable errors in the reconstructed topography. In fact, a larger threshold value can be selected to avoid such a situation and it is suggested that its value should not be greater than 5*u*σ/√*N*_{s}.

Then, a related question is how large the actual amplitude of the wave components in the isotherm should be to ensure that the calculated amplitude exceeds the threshold value of 3.5*u*σ/√*N*_{s}. In order to understand this question the probability of *P*(*C*_{n}>3.5*u*σ/√*N*_{s}|*C*_{tr}=*C*) was calculated using a Monte Carlo sampling of equation (5). As shown in Figure 10b, if the actual amplitude is 3.2*u*σ/√*N*_{s} then the calculated amplitude would exceed 3.5*u*σ/√*N*_{s} with a probability of 50%.

- © Geological Society of London 2009