Quantifying the rate and depth dependence of bioturbation based on optically‐stimulated luminescence (OSL) dates and meteoric 10Be

Both the rate and the vertical distribution of soil disturbance modify soil properties such as porosity, particle size, chemical composition and age structure; all of which play an important role in a soil's biogeochemical functioning. Whereas rates of mixing have been previously quantified, the nature of bioturbation's depth dependence remains poorly constrained. Here we constrain, for the first time, the relationship between mixing rate and depth in a bioturbated soil in northeast Queensland, Australia using a novel method combining OSL (optically‐stimulated luminescence) ages and meteoric beryllium‐10 (10Be) inventories. We find that the best fit mixing rate decreases non‐linearly with increasing soil depth in this soil and the characteristic length scale of 0.28 m over which the mixing coefficient decays is comparable to reported rooting depth coefficients. In addition we show that estimates of surface mixing rates from OSL data are highly dependent on erosion rate and that erosion rate must be constrained if accurate mixing rates are to be quantified. We calculate surface diffusion‐like mixing coefficients of 1.8 × 10−4 and 2.1 × 10−4 m2 yr−1 for the studied soil for two different estimates of soil erosion. © 2014 The Authors. Earth Surface Processes and Landforms Published by John Wiley & Sons Ltd.


Introduction
Soil is the living biogeochemical reactor that supports the majority of terrestrial life. In recent years there has been renewed interest in modelling of soil processes and functioning due to concern about the soil resource in the face of uncertain future climate and land use (Minasny et al., 2008). Models of pedogenesis are only as good as the mathematical descriptions of soil processes they contain, and the functional form of many soil processes remains poorly constrained. Constraining soil processes is often inhibited by the long timescales over which they operate. Here we focus on bioturbation, or the mixing and turnover of soil resulting from biological processes. The rate and depth dependence of mixing can fundamentally alter soil functioning, both by modifying the pore and particle size structure of the soil (Wilkinson et al., 2009) and by modifying the age structure of particles contained within the soil, which influences soil geochemistry (Mudd and Yoo, 2010). Some studies have proposed that soil disturbance declines non-linearly with depth (Humphreys and Field, 1998;Wilkinson and Humphreys, 2005;Wilkinson et al., 2009), whereas other field studies assume a constant mixing rate with depth (Cousins et al., 1999;Bunzl, 2002;Kaste et al., 2007). Recent studies have incorporated depth dependent mixing into modelling attempts (Vanwalleghem et al., 2013) and it has also been employed to analyse observations (Yoo et al., 2011). Yet no study has yet attempted to discriminate between these two end member scenarios using field data.
We show here how information from two tracers can be exploited in novel ways to estimate and constrain the process of bioturbation. First, we use optically-stimulated luminescence (OSL) dates in combination with a bioturbation model to evaluate both the model and the depth dependence of mixing. Secondly, we calculate a long-term erosion rate from a measured inventory of meteoric beryllium-10 ( 10 Be) and incorporate this into the bioturbation model in order to understand the importance of accounting for eroding profiles when interpreting observed soil properties.

Study site
These two geochemical tracers have been measured in a basaltic soil profile of the Hughenden region in northeast Queensland, Australia (Figure 1). The mean annual rainfall for this region is approximately 500 mm, with 80% falling in the summer monsoon months (October to March), and the climate is classed as semi-arid tropical. Vegetation is open savanna woodland, dominated by Eucalyptus and Acacia trees.
The soil profile is located at an elevation of 550 m, on a gently sloping (< 2°) plateau surface, called the Denna Plain (Coventry et al., 1985). The Denna Plain is underlain by Mesozoic quartz-rich sandstones, well exposed in gorges along nearby Porcupine Creek (Figure 1b), and is partly mantled by late Cenozoic sediments and basalt flows with K/Ar ages in the range 5.7 to 5.3 Ma (Coventry et al., 1985).
The sampled soil profile is located on a basalt flow, named the Wongalee Flow, dated at 5.6 Ma (Coventry et al., 1985). It is the oldest and deepest (130 cm) of several soils forming a chronosequence on basalt flows that were studied by Pillans (1997), who, using a simple relationship between soil age and profile depth, estimated a very low net soil production rate of 0.3 m Myr À1 . The soil is a Ferralsol (IUSS Working Group WRB, 2006) or a Lateritic Krasnozem (Stace et al., 1968). The soils of this region are colonized by termites, which actively transport soil to the surface in order to build mounds (Coventry et al., 1988). Termites are a dominant invertebrate in tropical soils (Wood and Sands, 1978); we thus believe this soil to be representative for a substantial fraction of semi-arid tropical and tropical soils.

OSL sampling
Optical dating of sediment is a technique which utilizes the OSL from constituent mineral grains to measure the time elapsed since the grains were last exposed to sunlight (the 'optical age'). The ages obtained from the luminescence signal can therefore provide information on the rate of mixing within the soil. The basis of the optical dating technique is that absorption of radiation leads to the accumulation of light-sensitive trapped charges in quartz crystals when they are stored in the dark, i.e. when buried. The population of trapped charges grows with time and is proportional to the total radiation dose accrued since the crystal was last exposed to light, i.e. when it was 'reset' during a residence period at the ground surface. In this work, OSL is used for grain-by-grain measurements of 'equivalent dose' (ED) to calculate the time since the individual quartz grains were last resident at the surface. The sensitivity of the OSL signal to light (Spooner, 1994) gives confidence that effectively complete bleaching of the OSL occurs for any grain exposed at the ground surface for in excess of several seconds. We chose the 'improved' single-aliquot regenerative-dose (SAR) protocol (Murray and Wintle, 2000) to measure the exposure history of each individual grain. It is important to note that because the underlying basalt of the studied soil profile does not contain quartz, quartz grains can only enter the soil via the surface and not from the underlying parent material. This greatly simplifies our OSL age profiles because chemical weathering at the soil/saprolite interface will not introduce infinite OSL ages to the soil.

Sample collection and preparation
Samples were collected in steel coring cylinders manually hammered into freshly-exposed profiles; with each tube containing approximately 0.1 kg sediment. Coarse grains of quartz were available in each sample, and as our aim was to utilize individual grains as tracers to reveal individual particle motion within the sediment, we chose grains of 350 to 425 μm in size rather than more conventionally used 200 μm grains. Both grain fractions approximated the median sediment particle size. The greater mass of quartz per grain provided a higher OSL output, which given the extreme faintness of typical single-grain OSL ensured better counting statistics and correspondingly better precision of the measured ages. Quartz grains of 350 to 425 μm were extracted from each sample in the laboratory under low-intensity red light using a procedure including, sequentially: hydrochloric acid (HCl) digestion, sieving, heavy liquid flotation (retaining < 2.68 g cm À3 fraction), and etching in 48% hydrofluoric acid (HF) for 40 minutes to remove both surface discoloration and the 6 À 8 μm α-particle irradiated outer layer.

Equivalent dose (ED) measurement
Measurements were made using a Type TL-DA-15 Risø TL/OSL reader, with optical stimulation of 420 to 560 nm provided by an optically-filtered halogen lamp. OSL emissions were detected by an EMI 9235QA photomultiplier tube optically filtered by two 2 mm Hoya U 340 filters. One single 350-425 μm fraction was attached to the centre of each of 48 stainless steel discs using silicone oil, prior to analysis by the SAR protocol (Murray and Wintle, 2000). Samples were illuminated for 125 seconds at 125°C. The OSL signal was extracted by subtracting an integral of the final 20 seconds of exposure from an integral of the first 20 seconds. Irradiations were performed using a 1.5 GBq 90 Sr/ 90 Y beta plaque incorporated in the apparatus and delivering a dose-rate of 0.097 ± 0.003 Gy s À1 . A 10 seconds preheat at 240°C was applied before each regenerative-dose OSL measurement. The dose-response curves were constructed for five different regenerative-dose levels ranging from 0.5 to 60 Gy, and a test dose of 3.0 Gy was used to track and allow correction for sensitivity changes occurring in the procedure. Uncertainties on ED are based on counting statistics and curve fitting, and also incorporate calibration uncertainties for the laboratory beta source. A total of about 42-49 individual grains (see Table I, number of grains sampled varies) were measured for each sample.

Dose-rate determination
Radioisotope concentrations in splits of the sampled sediment were measured by neutron activation analysis (NAA) for thorium (Th) and potassium (K), and delayed neutron activation (DNA) for uranium (U), at the Becquerel Laboratories, Lucas Heights Science and Technology Centre, Australia. Potassium concentration was also measured using X-ray fluorescence (XRF), by CSIRO, Canberra, Australia. The radioisotope concentrations within HF-etched quartz grains was assumed to be 10% of the bulk sediment concentration following Aitken (1985), and the efficiency of the internally-emitted alphaparticles in producing OSL was assumed as 0.05 ± 0.02, based on previous results for quartz-rich sediment (Questiaux, 1990). The cosmic ray dose-rate was calculated using the relevant environmental data and the equations of Prescott and Hutton (1994). Dose-rates are corrected for attenuation by soil water using the water content measured from the as-collected samples. The radioisotope data, along with the data used for cosmic ray calculation are shown in Table S1 (Supplementary Material), along with the calculated dose-rate components and total dose-rates for each sample.
Age evaluation was performed with the AGE program of Grün (1999), incorporating dose-rate conversion factors from Adamiec and Aitken (1998). Beta-particle attenuation factors used were 0.76 ± 0.01 for U, 0.70 ± 0.01 for Th and 0.86 ± 0.01 for K.
Estimating erosion with meteoric 10 Be Measurements of meteoric 10 Be in soils provide a robust method to estimate rates of soil loss by measuring the accumulated inventory of 10 Be within the soil profile (Willenbring and von Blanckenburg, 2010). The inventory method is particularly attractive because it does not assume any prior knowledge of the mixing rates and mechanisms in the soil or requires the presence of endogenous quartz. To estimate long-term erosion rates with meteoric 10 Be we must know the period of time that 10 Be has been accumulating in the soil (i.e. the age of the soil), the total inventory of 10 Be within the soil profile and an estimate of the rate of deposition on to the surface of the soil. A detailed review on the methods of estimating the flux of 10 Be to soils and calculations for estimating erosion rates can be found in Willenbring and von Blanckenburg (2010).
Because we know the age of the Queensland soil, we can calculate the rate of surface erosion from the total inventory of meteoric 10 Be. Meteoric 10 Be concentrations were measured by accelerator mass spectrometry (AMS), full details of sample preparation and AMS measurements can be found in Fifield et al. (2010). The procedure of deriving a long-term erosion rate from meteoric 10 Be measurements is explained in the Appendix. We tested two meteoric 10 Be flux scenarios: (i) a dilution scenario, where the flux of 10 Be to the soil surface is constant regardless of the rate of rainfall and (ii) an additive scenario where the flux of 10 Be is a function of the rate of precipitation.
In estimating the flux of 10 Be to the Hughenden soil profiles, we assume that there have been no major shifts in stratospheric-tropospheric exchange (STE) since the parent material was first exposed to meteoric 10 Be and that there is no longterm trend in solar intensity. Because we do not know the exact flux of 10 Be to these soils we must estimate the flux using the literature. For the dilution scenario where the 10 Be flux remains constant regardless of the rate of precipitation, the flux of 10 Be to the soil surface is obtained from Willenbring and von Blanckenburg (2010) who calculated the average flux from Field et al.'s (2006) and Heikkila et al.'s (2008) modelled fluxes. For northeast Queensland this is 5.0 × 10 5 atoms cm À2 yr À1 . Additions due to dust deposition are also accounted for. Mahowald et al. (1999) simulate dust deposition for northeast Queensland at 0.5 g m À2 yr À1 under the current climate. However, it is acknowledged that rates of dust deposition have likely changed over the course of soil development at these sites. Lal (2007) estimate a global average 10 Be concentration for dust of 6.5 (±3.5) × 10 9 10 Be atoms g À1 dust, based on a number of approaches including, ice core measurements and inland dust samples. Using the rate of dust deposition over the period of soil development from Mahowald et al. (1999) and Lal's (2007) estimated 10 Be concentration per gram of dust an extra 3.25 × 10 5 atoms cm À2 yr À1 is deposited onto the soil surface. For the additive scenario where the 10 Be flux is proportional to the rate of precipitation, the same rate of 1.5 × 10 4 atoms 10 Be cm À3 of rainfall used by Fifield et al. (2010) for sites in northwest and southeast Australia is adopted. The assumption that the rate of precipitation (500 mm yr À1 ) at the site has remained constant over the course of soil development is supported with results from palaeoclimatic studies of the region. Pollen analyses at sites on the Atherton tableland, 300 km north of Hughenden, suggest that glacial periods during the last 230 000 years were drier than present and inter-glacial periods experienced similar rainfall . Pollen analyses of offshore sediment cores yield similar conclusions (Moss and Kershaw, 2007). In the same Atherton Tableland region but at an older site (1-5 Ma), Kershaw and Sluiter (1982) concluded from pollen results that rainfall has remained similar to present. The sediment record at Lake Buchanan, some 200 km southeast of Hughenden, indicates fluctuating lake levels over the past 1.8 Ma, consistent with wetter and drier periods, but with no overall long-term trend towards wetter or drier conditions (Chivas et al., 1986). We thus reason from these few palaeoclimatic reconstructions that although there are likely inter-glacial fluctuations in rainfall (drier in glacial periods, and amounts similar to present during inter-glacials), there is no evidence to suggest a long-term trend to wetter or drier conditions over the last 5 Ma.
For the additive flux scenario the rate of delivery of 10 Be for these sites is 1.075 × 10 6 atoms cm -2 yr -1 and the constant flux rate is 0.825× 10 6 atoms cm -2 yr -1 . We acknowledge that there are some large assumptions associated with this method of calculating erosion. First, it is assumed that the soil density of the profiles is uniform (ρ in Equation (A6), see Appendix). Secondly, it is also assumed that the surface concentrations of 10 Be (N surf in Equation (A6), see Appendix) are at a steady state through time. Like Fifield et al. (2010) we estimate that the largest uncertainty is associated with the flux of 10 Be to the soil surface (i.e. assumptions of constant rainfall and dust additions) and allocate 25% uncertainty to this, a further 10% uncertainty is allocated for the assumptions associated with soil density and the concentration of 10 Be in surface layers lost to erosion.

Modelled age tracer
Optical dates effectively provide a clock which records the time since a soil parcel was last exposed to sunlight and we show here how this feature can be exploited in a novel way with a modelled age tracer.
A number of authors have reasoned that the probability of soil disturbance decreases with increasing depth because the number and activity of biological agents in the soil (e.g. roots, worms, termites, burrowing mammals) should decline with depth (Humphreys and Field, 1998;Cousins et al., 1999;Wilkinson and Humphreys, 2005;Vanwalleghem et al., 2013). Following Cousins et al. (1999) we propose that bioturbation can be described as a diffusive process with depth dependent diffusivity: where y is a soil property such as a constituent element, D is the diffusion coefficient (in m 2 yr À1 ) and z is depth (in metres).
To determine if a diffusive parameterization is in fact a realistic mechanism for modelling bioturbation and to estimate the depth dependence of the diffusivity, an age tracer is introduced to this model which mimics OSL age, A, directly, by accumulating time of a soil particle since it was last at the soil surface. The time evolution equation for A is the same as the transport equation for any soil property apart from one additional term which takes into account that once a soil particle is beneath the soil surface the soil particle starts to age. To reproduce this ageing mathematically, a constant of one is added to the diffusion term of the transport equation: If there were no diffusion, the solution of the equation would simply be A = t. In the presence of diffusion, mixing leads to a redistribution of ages represented by the first term on the right of Equation (2). The technique we describe here is the same as that used to model the age of interior ocean water in ocean general circulation models (e.g. England, 1995).
The equation is solved for a boundary condition of A(0) = 0. The boundary condition of zero is set for the middle of the first depth level instead of the top because the upper surface of a soil profile is constantly subject to rain splash and scour and therefore the soil grains in the upper 50-100 mm are regularly exposed to sunlight (Wilkinson and Humphreys, 2005). The 1191 QUANTIFYING THE RATE AND DEPTH DEPENDENCE OF BIOTURBATION bottom boundary condition allows no mixing into the bedrock so at D(n), ∂A/∂z is equal to zero, where n is the total number of model soil layers. We discretize the equation with dz = 0.1 m so that there are h/dz soil layers where h is the soil thickness of the sampled soil profile. The soil in the surface level retains an age of zero and the rest of the profile increments in age by the model time step. The diffusion process means that indefinite growth of the tracer age is prevented because the zero age signal at the surface is continuously mixed into the interior. More generally porosity would need to be included in Equation (2), i.e. D(z) would be replaced by D(z).ϕ(z), where ϕ is fractional porosity; here we assume constant porosity with depth which is a good approximation in these well weathered soils. We assume here that mixing in the soil profile has reached a stationary state (i.e. time invariant) and thus solve for 0 (1) and (2) may be either constant with increasing soil depth or decline with soil depth to reproduce the assumed decrease in biological activity with increasing soil depth. In this study we test two forms of the relationship between D and depth. The first considers an exponential decrease in the diffusion coefficient with depth: where z b is the e-folding length scale (in metres). The second adopts a linear function: where a is the gradient of the slope, which includes the case of a constant diffusion coefficient (= 0). The diffusion equation is solved numerically. Vertical profiles of the modelled ages for both formulations of the distribution of D are compared with the measured OSL ages from quartz grains in the Hughenden soil profile. Optimized values are fitted to D(0), z b and a based on minimum root mean square error. So far, the prediction of the diffusivity rate and distribution is based on the assumption that the soil is not eroding. However, it is likely that an eroding surface will introduce bias to the estimates neglecting erosion.
Equation (2) is valid in the limiting case of no erosion. However, erosion can have a strong control of the age structure of mixed soils. Erosion samples particles from the top layers of the soil, and mixing determines the age population of these sampled particles (Mudd and Yoo, 2010). We account for this effect by introducing an erosion term to Equation (2). The erosion formulation acts by advecting the age profiles up towards the profile surface at a rate T, effectively removing the surface layers (Kirkby, 1985). Thus, Equation (2) becomes (see Appendix for derivation of the advection term). To determine the influence of erosion on the interpretation of OSL age profiles, the long-term erosion rate for the soil profile estimated from the measured meteoric 10 Be inventory is used.

Results
The mean OSL ages of the individual quartz grains at each sample depth are shown in Table I. We have made the assumption that the mean is representative of the data (histograms of the OSL ages at each sample depth are provided in the Supplementary Material, Figure S1). The luminescence ages of the quartz grains increase with increasing soil depth. If mixing was absent from the soil the OSL dates would not be dateable below the surface and only the surface grains would show any bleaching. The vertical distribution of the OSL ages implies that down to the depth measured there is no barrier to biological activity or mechanical mixing. Figure 2 shows the observed ages versus the optimized model ages for the exponential and linear decreases in D with soil depth. The optimal D(0) is equal to 9.81 × 10 À5 m 2 yr À1 for the exponential function with the characteristic length z b equal to 0.28 m. For the linear function D(0) is equal to 6.9 × 10 -5 m 2 yr À1 and a is 1.06 × 10 À4 m yr -À1 . The normalized mean bias (NMB), defined is the observed value and n is the number of vertical observations, of the modelled age versus the fractional error-weighted mean OSL ages is À1.54% for the exponential function and 12.33% for the linear decrease in D. Therefore, the modelled ages agree best when an exponential decline in the diffusion coefficient with increased soil depth is adopted. Figure 3 shows the age profiles achieved for a range of erosion rates with the same exponential mixing parameters esti- Figure 2. Boxplots of observed optically-stimulated luminescence (OSL) ages at sampled depths and simulated soil ages for the Hughenden soil profile. The crosses are outliers and the filled squares are the fractional error-weighted means of the OSL ages. Both linear and exponential formulations of the modelled diffusion coefficients are displayed. Figure 3. Steady-state age profiles simulated with the modified age tracer model (Equation (5)) for a range of erosion rates (T). mated from the OSL ages and Equation (2). Eroding profiles result in older soil particles at equivalent depths. This is due to the removal of the younger surface layers which brings the deeper, older soil particles closer to the soil surface. This means that the mixing rate in an eroding soil will be higher than that suggested by the previous model in order to match the observed OSL ages. We constrain the rate of mixing in the Hughenden soil profile using Equation (5), which combines OSL and 10 Be data.
The long-term erosion rates for this soil calculated from the measured meteoric 10 Be inventory for two flux scenarios are shown in Table II. We demonstrate that for this soil, comparable rates of erosion are achieved for both scenarios.
These estimated rates of erosion rely on the assumption that the only losses of 10 Be are via radioactive decay and soil erosion and that there has been no losses from leaching in solution. Figure 4 shows the distribution of the meteoric 10 Be in the Hughenden soil profile. The distribution suggests that the 10 Be has been transported within the profile. If soil mixing is the only mechanism of 10 Be transport in the soil, the expected distribution is a steady decline with increasing soil depth, whereas the distribution observed suggests that 10 Be must also move through the soil by another mechanism. Figure 4 illustrates the close association of 10 Be concentration with the clay content of the soil; the Pearson correlation coefficient is 0.89 (significant at p < 0.05). Beryllium-10 is highly particle-reactive at neutral to alkali pH (You et al., 1989) and the close to neutral pH of the Ferralsol soil profile suggests that any 10 Be which is transported in solution will immediately readsorp onto the fine clay particles, therefore minimum losses of 10 Be from the soil profile are expected. Another assumption made when calculating these erosion rates is that the whole inventory of 10 Be has been measured. However, Fifield et al. (2010) found 10 Be present in saprolite at one of their sites, which indicates that there was some transport of 10 Be in solution in their soil. They propose that the 10 Be moves in solution by repeated adsorption on and desorption from soil grains as the water percolates through the soil and eventually slowly infiltrates into the saprolite rather than being lost from the profile in groundwater flow. Saprolite was not sampled in this study so it is possible that the whole inventory of 10 Be has not been measured. Thus, accounting for the possibility that some 10 Be is lost by leaching and that the whole inventory of 10 Be has not been captured, the rates of erosion presented in this study are likely overestimates.
To determine the new diffusivity parameters needed to match the OSL ages we assign the values from Table II to T in Equation (5). We assume that erosion does not alter z b , the e-folding length scale, so this remains equal to 0.28 m. The new values of D(0) for soil eroding at rates of 5.0 and 6.8 m Myr À1 are 1.8 × 10 À4 and 2.1 × 10 À4 m 2 yr À1 , respectively. Table I shows the diffusion coefficients for each depth interval of the soil profile. These surface diffusion rates are approximately a factor of 2 higher than the rate calculated when we did not account for erosion, demonstrating the importance of accounting for erosion rates when constraining mixing rates.

Discussion
We find that a diffusion model with depth dependent mixing rates is capable of simulating realistic soil mixing within the Queensland soil profile. The Ferralsols of this region are colonized by termites. These mound builders bring grains to the surface from deeper in the profile. Eventual collapse and subsidence of the termite mounds returns the grains to the soil profile. For this soil dominated by biotic activity an exponential decline (e Àz=z b ) in the diffusion coefficient fits the data best.
The modified surface mixing rates fall within the range 1-2 × 10 À4 m 2 yr À1 estimated by Kaste et al. (2007) for soils in both south-eastern Australia and California. The rates of   Kaste et al. (2007) were calculated from measurements of short-lived nuclides which suggests that the short-timescale processes of termite mixing may be the dominant mixing agent compared with longer timescale processes such as tree-throw. The z b value of 0.28 m is very similar to the vertical root scaling for tropical evergreen biomes (0.26 m) obtained from the global study of tree distributions by Jackson et al. (1996). Although it is not possible to say with certainty from just one site, this study suggests that interestingly the distribution of biological activity in soils mirrors the vertical distribution of roots. This could be because the termites are distributed relative to the abundance of a food source (dead root material for example). Thus, mixing in models could be parameterized to correspond to rooting distributions as there is more data available with which to constrain this trend (e.g. Jackson et al., 1996).
We have shown that OSL ages in combination with a modelled age-tracer provide a robust means of evaluating the mechanism and rates of bioturbation. This study introduced two novel methods to achieve these results. First, the use of an age-tracer in a model to replicate OSL ages of quartz grains and secondly the combination of meteoric 10 Be measurements to understand the importance of accounting for erosion in such studies. We demonstrate that a diffusive process of mixing with a non-linear decline in the diffusion coefficient with increasing soil depth can reproduce the OSL age profile. We also highlight the influence that erosion has in altering these age profiles and how this must be accounted for in order to estimate realistic mixing in soils.
We are confident that bioturbation is the dominant process responsible for the observed variability in OSL ages. Previous work by one of us (NAS), reported in Heimsath et al. (2002), showed a significant proportion of grains to have no detectable OSL. On subsequent analysis using thermoluminescence measured at the 'hard-to-bleach' 375°C peak, these were found to be thermoluminescence-saturated ('infinite age'), which was attributed to their saprolitic origin and not having yet received exposure to light at the ground surface. However, in the present work, almost all of the 48 grains measured per sample gave usable OSL; those that did not were two discs from 27 to 29 cm, three discs from 46 to 50 cm and six discs from 57 to 60 cm. No further work was performed on these unresponsive discs, given their very low proportion, and the possibility that some or all of these null results may reflect absence of grains through loss from those discs, In this study a fundamental assumption is that the measured ages of quartz grains actually represent the time elapsed since last surface exposure. The validity of this assumption depends on the occurrence of a uniform dose rate down the soil profile. Inspection of Table I shows that the variation in total dose-rate between samples is negligible in comparison with the magnitude of the age variations seen within each individual sample. The maximum and minimum dose rates are 1.00 and 0.85 Gy ka À1 , respectively, with a mean dose rate of 0.93 ± 0.06 Gy ka À1 . The cosmic ray dose rate was calculated for each sample using the burial depth at the time of sampling, and hence does not account for changes in this dose-rate component over time due to the grain's trajectory in the sediment. However, the cosmic ray dose rate for the uppermost sample is only 0.01 Gy ka À1 greater than for the lowermost sample, i.e. barely 1% of the total dose rate, and hence the variation of this dose-rate component down the profile is insignificant here (see Supplementary Material, Table S1).
We conclude that the range of ages seen in any given sample here is not a consequence of varying dose rates, regardless of the trajectories that grains may have followed through the soil. This potentially complicating issue was also noted by Bateman et al. (2003, p. 150), and in a more recent study by Stockmann et al. (2013), who found similar dose rate variations to the present work, reporting that 'total dose rates … did not differ by more than 10% down the profile', and similarly concluded that this removed this 'potential complication'.
The results of this study are relevant to a number of Earth surface processes. Besides the contribution of bioturbation to geomorphological processes such as hillslope evolution (Kirkby, 1985;Heimsath et al., 2002), bioturbation also plays an important role in biogeochemical processes, especially nutrient and carbon cycling. The soil body, for example, is the largest terrestrial reservoir of organic carbon (Jobbágy and Jackson, 2000) so it is important that we can quantify the transport processes within soil in order to predict carbon cycle climate feedbacks. The importance of accounting for buried carbon is recognized and attempts are increasingly being made to resolve the vertical distribution of carbon in processoriented, predictive models (Elzein and Balesdent, 1995;Jenkinson and Coleman, 2008;Koven et al., 2013). However, the burial and redistribution of carbon within the soil due to bioturbation is either ignored as in the RothC model (Jenkinson and Coleman, 2008) or when bioturbation is included (Elzein and Balesdent, 1995;Koven et al., 2013) the mixing rate is assumed constant with depth due to lack of evidence suggesting otherwise. The results of this study are particularly significant for modelling processes of soil carbon distribution with depth in the Tropics, because termites are a dominant mixing agent in these ecosystems. Tropical forests and savannas account for about 25% of global soil carbon stocks (IPCC et al., 2000). The ability to realistically model soil carbon in these systems is therefore important for predicting future land-climate interactions.
It is difficult to apply this method to evaluate mixing in a wider range of soils due to the lack of OSL data from soil profiles. However, we have demonstrated how extra information can be obtained from OSL age profiles using this powerful method and we hope that this will encourage further research into this field.

Conclusions
We have, for the first time, tested two end member scenarios for the functional relationship between mixing rates and soil depth, and found that an exponentially decreasing diffusion coefficient explains the data better than a linearly decreasing or a constant diffusion coefficient. We were able to quantify the length scale of the decline in these diffusion-like mixing coefficients, and this length scale is similar to the length scale describing the decline in root density with depth, estimated from a global dataset. We have also outlined examples of problems that can be addressed with these results.

∂A ∂z
Àδz ð Þþ ∂A ∂t δt ¼ 0 ( A 2 ) or finally Estimating rates of erosion from meteoric 10 Be The total inventory of meteoric 10 Be is calculated from observations of 10 Be in a soil profile by where I is the total inventory of 10 Be atoms (in atoms cm À2 ) in a soil from the surface to the measured soil depth z m (in centimetres), N(z) is the measured number density of 10 Be at depth z (in atoms g À1 ) and ρ is the soil density (in g cm À3 ). The total inventory of 10 Be in a soil profile at time t is the difference between the total flux into the soil by deposition from the atmosphere minus and losses from radioactive decay and surface erosion. Thus where F is the flux of 10 Be to the soil surface (in atoms cm À2 yr À1 ) at time t (here, F, is assumed to be constant over time), ε is the rate of surface erosion (in cm yr À1 ), N surf is the concentration of 10 Be in surface layers lost to erosion (in atoms g À1 ) (which is assumed to be at steady-state) and λ (in yr À1 ) is the radioactive decay constant. The value of λ is equal to ln2/t ½ where t ½ is the half-life of 10 Be (1.39 × 10 6 , Korschinek et al., 2010). The solution of Equation (A6), for constant F and N surf , solved for the long-term, average erosion rate, ε, is ε ¼ 1 ρN surf F À Iλ 1 À e Àλt (A6) (Willenbring and von Blanckenburg, 2010).
Acknowledgements-This work was funded by a NERC studentship awarded to MOJ. The authors are very grateful to Jon Lloyd for his constructive comments on earlier stages of the paper. The invaluable assistance of Mariana di Tada with sample preparation for the 10 Be samples is also gratefully acknowledged. The authors also thank two anonymous reviewers for their very helpful and constructive comments.