The rarefied (non-continuum) conditions of tracer particle transport in soils, with implications for assessing the intensity and depth dependence of mixing from geochronology

We formulate tracer particle transport and mixing in soils due to disturbance-driven particle motions in terms of the Fokker–Planck equation. The probabilistic basis of the formulation is suitable for rarefied particle conditions, and for parsing the mixing behavior of extensive and intensive properties belonging to the particles rather than to the bulk soil. The significance of the formulation is illustrated with the examples of vertical profiles of expected beryllium-10 (¹⁰Be) concentrations and optically stimulated luminescence (OSL) particle ages for the benchmark situation involving a one-dimensional mean upward soil motion with nominally steady surface erosion in the presence of either uniform or depth-dependent particle mixing, and varying mixing intensity. The analysis, together with Eulerian–Lagrangian numerical simulations of tracer particle motions, highlights the significance of calculating ensemble-expected values of extensive and intensive particle properties, including higher moments of particle OSL ages, rather than assuming de facto a continuum-like mixing behavior. The analysis and results offer guidance for field sampling and for describing the mixing behavior of other particle and soil properties. Profiles of expected ¹⁰Be concentrations and OSL ages systematically vary with mixing intensity as measured by a Péclet number involving the speed at which particles enter the soil, the soil thickness, and the particle diffusivity. Profiles associated with uniform mixing versus a linear decrease in mixing with depth are distinct for moderate mixing, but they become similar with either weak mixing or strong mixing; uniform profiles do not necessarily imply uniform mixing.

1 Introduction

Soils on Earth's surface are granular materials consisting of polymineralic clasts and individual mineral grains, organic matter and live biota. These materials experience patchy, intermittent mixing motions associated with disturbances due to bioturbation (Darwin, 1881; Shaler, 1891; Gabet, 2000; Reichman and Seabloom, 2002; Meysman et al., 2006; Wilkinson et al., 2009; Covey et al., 2010; Astete et al., 2015), the effects of frost and ice growth and thaw (Branson, 1992; Matsuoka and Moriwaki, 1992; Auzet and Ambroise, 1996; Branson et al., 1996; Harris et al., 1997; Matsuoka, 1998; Anderson, 2002), and the swelling and shrinking of certain minerals with wetting and drying (Eyles and Ho, 1970; Fleming and Johnson, 1975). In addition, these soil materials may undergo mixing motions in relation to the chronic creation and relaxation of disordered granular structures (Hsiau and Hunt, 1993; Utter and Behringer, 2004; Fan et al., 2015) associated with granular creep (Houssais et al., 2017; Ferdowsi et al., 2018).

Soil particle mixing is a key process in soil formation (Shaler, 1891; Birkeland, 1984; Wilkinson et al., 2009) and in its associated ecological role of “modifying geochemical gradients, redistributing food resources, viruses, bacteria, …and eggs” (Meysman et al., 2006), as well as being responsible for redistributing substances, including contaminants, attached to particles (Cousins et al., 1999; Covey et al., 2010; Astete et al., 2015). Moreover, the idea of disturbance-driven transport and mixing of soil particles is central to current treatments of soil creep (Culling, 1963; Roering et al., 1999, 2002; Gabet, 2000; Anderson, 2002; Gabet et al., 2003; Furbish, 2003; Roering, 2004; Furbish et al., 2009b, 2018a), the slow but steady bulk motion of soils on hillslopes, where the influence of gravity gives a downslope bias to particle motions. Because of the significance of soil particle mixing in numerous problems spanning ecological to geomorphic timescales, there is a continuing, compelling need to fully clarify the kinematics, and eventually the mechanical basis, of soil particle motions during transport and mixing (Furbish et al., 2009b, 2018a, b; BenDror and Goren, 2018; Ferdowsi et al., 2018).

Currently it is not possible to directly measure disturbance-driven particle motions and associated mixing in the setting of a natural soil (although this is entirely possible in experiments and numerical simulations of granular creep, Utter and Behringer, 2004; Kamrin and Koval, 2012; Fan et al., 2015). Moreover, we do not yet have a mechanical theory to describe these motions given the complexity – notably the biotic complexity – of phenomena involved in disturbances and associated particle displacements (Furbish et al., 2009b, 2018a). Thus, as in studies of particle mixing associated with marine bioturbation (Boudreau, 1986a, b; Boudreau and Imboden, 1987; Teal et al., 2008; Lecroart et al., 2010), a key strategy to clarify the nature of particle motions and mixing in soils involves using tracer particles identified by specific physical or chemical properties. Two tracer properties have emerged in the field of geomorphology as being of particular interest: in situ cosmogenic radionuclide (CRN) concentrations and optically stimulated luminescence (OSL) particle ages (Granger and Riebe, 2014; Heimsath et al., 2002; Johnson et al., 2014). Cosmogenic nuclides continually accumulate within minerals due to cosmic ray interactions with mineral atom nuclei, for example, producing ¹⁰Be from spallation of oxygen nuclei. Using luminescence systematics, the time elapsed since luminescence-sensitive particles were last exposed to light or heat at the soil surface is estimated from the luminescence signal that accumulates within the crystal lattice in response to a combination of ionizing radiation emitted from the decay of radioactive elements in the surrounding soil and cosmic radiation (Rhodes, 2011). Particles that accumulate CRN atoms or luminescence signals during their complex motions within soils – thereby serving as tracer particles – are naturally occurring (as opposed to being “seeded”) and therefore behave mechanically the same as other soil particles. As a consequence, CRN and OSL tracer particles are particularly relevant in assessing particle mixing over timescales of soil formation and transport in the context of landform and landscape evolution.

Building on the pioneering work of Lal (1991) concerning the relation between rock erosion rates and the in situ production of cosmogenic radionuclides, vertical profiles of CRN concentrations in soils and underlying saprolite are now used to calculate soil production rates (e.g., Heimsath et al., 1997, 2000, 2005, 2012; Small et al., 1999; Anderson, 2002; Wilkinson et al., 2005) as well as to infer the intensity of soil particle mixing in the presence of mechanical and chemical erosion (Small et al., 1999; Schaller et al., 2009; Granger and Riebe, 2014; Furbish et al., 2018b). Similarly, profiles of particle OSL ages are used to assess particle mixing (Heimsath et al., 2002; Wilkinson and Humphreys, 2005; Johnson et al., 2014; Furbish et al., 2018b). Because profiles of CRN concentrations and OSL ages inform descriptions of soil transport and interpretations of the delivery of CRNs to channels (Heimsath et al., 2002; Anderson, 2015; Furbish et al., 2018b), and associated interpretations of erosion rates at catchment scales (e.g., Brown et al., 1995; Bierman and Steig, 1996; Granger et al., 1996; Granger and Riebe, 2014; Granger and Schaller, 2014; Lukens et al., 2016), there is merit in further clarifying what these profiles reveal about particle mixing in soils.

It is now conventional to conceptualize certain soil particle mixing motions as a diffusion-like process (Furbish et al., 2009b, 2018a, b; Campforts et al., 2016), building on the pioneering work of Culling (1963), who first pointed to the idea that soil particles undergo Gaussian diffusion in response to small disturbances. Various studies have thus appealed to some form of a diffusion equation or an advection–diffusion equation (Cousins et al., 1999; Covey et al., 2010; Stang et al., 2012; Johnson et al., 2014; Furbish et al., 2009b, 2018a, b; Astete et al., 2015; Campforts et al., 2016; Gray, 2018) to describe transport and mixing for comparison with measured vertical profiles of tracer particles in soils, notably including in situ CRN concentrations and particle OSL ages. But herein arises a need for caution, and clarity.

As described in Sect. 2, natural tracer particles – quartz particles in particular – occur under rarefied conditions, where it is unclear that a description of particle mixing based on a diffusion or advection–diffusion equation formulated for continuum conditions is satisfactory. Moreover, we often are interested in the transport of quantities that are associated with the particles and are not in themselves subject to advection and diffusion as normally envisioned to occur in a continuum. This includes particle CRN concentrations and OSL ages. Rather, such quantities might experience advection and diffusion, but only indirectly via the motions of the particles with which the quantities are associated. Within this context, our objectives in this mostly theoretical contribution are fivefold.

First, we illustrate why quartz tracer particles in soils experience transport and mixing under rarefied (non-continuum) conditions, and why it therefore becomes important to treat transport and mixing probabilistically, in a manner that formally appeals to the statistical mechanics idea of ensemble-expected (average) quantities. Our focus on quartz particles is purposeful, as these are ideal targets for in situ production of ¹⁰Be atoms, and for accumulating OSL signals. Second, we illustrate how the probabilistic basis of the Fokker–Planck equation, versus an “ordinary” continuum-like advection–diffusion equation, is well suited to the problem of rarefied conditions. Third, because extensive and intensive properties such as particle volume, ¹⁰Be concentration and OSL age “belong” to individual particles, not to the bulk soil, we illustrate why the probabilistic basis of the Fokker–Planck equation is suitable for parsing the mixing behavior of these properties – as opposed to assuming de facto a continuum-like mixing behavior in which these properties are assigned to the bulk soil. Fourth, we provide complementary numerical analyses that reveal important information not readily apparent in the analytical formulations, including an illustration of the variability in ¹⁰Be concentrations and OSL ages of individual particles in soils, with implications for interpreting field-based measurements. This part of the paper highlights a benchmark situation involving a one-dimensional mean upward soil motion with nominally steady surface erosion in the presence of either uniform or depth-dependent particle mixing, and varying mixing intensity. Fifth, we use the results for this benchmark case in relation to published field-based measurements to suggest constraints on assessing the intensity and depth dependence of mixing.

Note that, in the formulations presented below, we use full functional notation throughout. This provides clarity in how random variables, parameters, and moments of random variables depend on position and time, as well as how random variables might covary.

2 Rarefied versus continuum particle conditions in soils

Quartz particles targeted in sampling for ¹⁰Be analysis typically are within the range of 0.25–0.50 mm, but sometimes grains as small as 0.125 mm and as large as 0.85 or 1 mm are sampled from quartz-poor source materials (Gosse and Phillips, 2001; Morgan et al., 2011; Shakun et al., 2018). Quartz particles targeted for single-grain OSL analysis typically are within the range of 0.35–0.425 mm (e.g., Heimsath et al., 2002; Johnson et al., 2014), but smaller grains sometimes are used. Thus, neglecting aeolian inputs, target grains represent a subset of the total population of quartz grain sizes in soils released from parent bedrock during soil formation. In the following discussion we consider for illustration a single particle size, with recognition that the ideas extend to other particles.

Figure 1 Definition diagram of soil-mantled hillslope with mechanically active soil thickness $h=\mathit{\zeta }-\mathit{\eta }$, and cutout soil element with dimensions XYh. Bedrock material is continually transformed into soil by chemical and mechanical processes, and soil particles are transported downslope by creep or surface erosion.

Consider a soil element with dimensions XYh, where $X=Y=h=\mathrm{1}$ m, residing on a soil-mantled hillslope (Fig. 1). If in the ideal this element contains uniform particles of diameter d=1 mm that are approximately closely packed, then the total number of particles in the soil element is O(10⁹). Each cubic centimeter contains O(10³) particles. The average spacing is approximately equal to one particle diameter, and the geometrical mean free path λ (Furbish et al., 2009b) is a fraction of the particle diameter. For comparison, the number of molecules in a cubic centimeter of air, a continuum material at ordinary pressure–temperature conditions, is O(10¹⁹). The mean free path, which varies inversely with the molecular collision frequency or number density, is O(10⁻⁷) m, approximately 10³ larger than the effective molecular diameter. Assuming the continuum hypothesis is satisfied for a value of the Knudsen number $\mathit{\text{Kn}}=\mathit{\lambda }/L\le \mathrm{0.01}$, then the averaging length scale L defining a continuum physical point for air is O(10⁻⁵) m, far smaller than most scales of interest in treating particle transport and mixing in air.

For a soil developed from granitic bedrock, 20 % to 60 % of the volume of particles are quartz particles, some larger than 1 mm in diameter and many smaller. Per unit volume, the number of quartz particles targeted in sampling for ¹⁰Be analysis thus is generally smaller than the close-packed value of O(10³) cm⁻³ estimated above, and the average spacing may be on the order of millimeters to a centimeter or more. For example, in practical terms, ¹⁰Be analysis requires about 10 g of quartz. Assuming 0.5 mm grains, this represents ∼ 60 000 grains. For soils formed on granitic bedrock, one typically samples at least 1 L of soil for 10 g of quartz. This translates to n∼10–100 grains per cubic centimeter. The associated geometrical mean free path is about −10 cm, and the average spacing ${\mathit{\lambda }}_{\mathrm{s}}\sim ({\mathrm{cm}}^{\mathrm{3}}/n{)}^{\mathrm{1}/\mathrm{3}}$ is 0.2–0.5 cm. Although the behavior of tracer particles is unlike gas particle kinetics, we nonetheless can use these quantities in analogy with the mean free path. Conservatively using the average spacing as a suitable measure of the particle number density, then to satisfy the Knudsen condition of Kn≤0.01 in order to appeal to a continuum description of particle behavior, the averaging length scale L may approach the soil thickness. This condition is exacerbated if the parent bedrock is quartz poor. In addition, a small fraction – only a few percent – of quartz particles initially released from bedrock are sensitive to OSL and develop a less-than-saturated luminescence signal following exposure to sunlight or heat at the soil surface. Thus, tracer particles identified as those possessing a finite OSL age (Heimsath et al., 2002; Johnson et al., 2014) may be highly rarefied.

We therefore must admit at the outset that the number concentration of target quartz particles does not necessarily satisfy the continuum hypothesis. Nonetheless, we wish to use continuum-like formulations of transport and mixing of particle concentrations and associated quantities, that is, where particle concentrations, ¹⁰Be concentrations and particle OSL ages may be viewed as continuously differentiable functions of position and time. In order to justifiably do this, we therefore appeal to the idea of an ensemble of particle configurations, a statistical mechanics idea designed to treat rarefied particle conditions.

For an element of soil with dimensions XYh (Fig. 1), let f_z(z,t) denote the probability density function of particle positions z within the element. Thus, f_z(z,t)dz represents the probability that a particle is located within the small interval z to z+dz at time t. This represents an ensemble-expected value, as follows. We envision, as did Gibbs (1902), a great number (an ensemble) of nominally identical but independent systems, each containing a large number N of particles and behaving in a statistically similar manner with respect to transport and mixing. The expected number of particles within the interval z to z+dz in any system (realization) at time t may vary from one system to another. However, we then imagine taking the expected value over the ensemble (Kittel, 1958), akin to ensemble Reynolds averaging (Monin and Yaglom, 1971). This represents the expected number of particles within dz, namely, Nf_z(z,t)dz, where f_z(z,t) now is interpreted as the ensemble-expected density. Moreover, we may assume that f_z(z,t) is a smooth, continuous function. Further details regarding rarefied versus continuum conditions and ensemble averaging are provided in Appendix A.

In the developments below, we also consider joint probability density functions, for example, the joint density $f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)$ of individual particle volumes V_p, ¹⁰Be atom number concentrations n_p and positions z. We similarly assume that these represent ensemble-expected densities with respect to z. In principle, therefore, we are considering the expected concentration of particles and associated properties within any small interval z to z+dz in a soil element with dimensions XYh (Fig. 1), where averaging is over an ensemble of nominally identical but independent systems. In practical terms, one hopes to sample over an area XY such that the number of particles within any small interval z to z+dz is sufficiently large to provide reasonable estimates of ensemble-averaged values, where these estimates vary approximately smoothly over z and average over the effects of patchy, intermittent particle motions. However, this cannot be known a priori, a point to which we return below.

3 Formulation

3.1 Tracer particles

Consider a set of tracer particles that are undergoing transport and mixing within a soil. Here we initially restrict this set to chemically resistant quartz particles. Nonetheless, this set could consist of particles defined by other mineralogies, or it could be defined as the subset of quartz particles of a given size that possess a specified ¹⁰Be concentration or finite OSL age. For simplicity, and in anticipation of further analyses below, we focus on one-dimensional motions parallel to the z axis.

As above, let f_z(z,t) denote the probability density function of tracer particle positions z. Following Furbish et al. (2009b, 2018a, b), this density satisfies a Fokker–Planck equation of the form

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial f_z(z,t)}{\partial t}}\\ \text{(1)}& {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)f_z(z,t)-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial f_z(z,t)}{\partial z}}\right],\end{array}$$

where w_p(z,t) denotes the ensemble-averaged particle velocity (sometimes referred to as the “drift speed”) and κ_z(z,t) denotes the ensemble-averaged particle diffusivity. Specifically, let $r=z(t+\mathrm{d}t)-z\left(t\right)$ denote a particle displacement during the small interval of time dt. Then let $f_r(r;z,t)$ denote the probability density function of displacements r. The particle velocity w_p(z,t) is then defined kinematically as

$$\begin{array}{}\text{(2)}& w_{\mathrm{p}}(z,t)=\underset{\mathrm{d}t\to \mathrm{0}}{lim}{\displaystyle \frac{a(z,t)}{\mathrm{d}t}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}r{f}_{\mathrm{r}}(r;z,t)\mathrm{d}r,\end{array}$$

and the particle diffusivity κ_z(z,t) is defined as

$$\begin{array}{}\text{(3)}& {\mathit{\kappa }}_z(z,t)=\underset{\mathrm{d}t\to \mathrm{0}}{lim}{\displaystyle \frac{a(z,t)}{\mathrm{2}\mathrm{d}t}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{r}^{\mathrm{2}}{f}_{\mathrm{r}}(r;z,t)\mathrm{d}r,\end{array}$$

where a(z,t) denotes the particle activity probability, effectively the proportion of time that particles are in motion (Furbish et al., 2009a, b, 2016, 2018a).

This formulation assumes Gaussian diffusion of particles. Interestingly, Culling (1963) first pointed to the idea that soil particles undergo Gaussian diffusion in association with particle concentration gradients, in response to small disturbances. Culling developed his ideas from kinetic theory and statistical mechanics, borrowing the description of Brownian motion due to Einstein (1905) and the formulation of a particle diffusion-like equation due to Chandrasekhar (1943), both of which start from the master equation (Risken, 1984; Ebeling and Sokolov, 2005; Furbish et al., 2009a, b). Culling's formulation has for decades provided the inspiration for conceptualizing what now are referred to as “disturbance-driven” particle motions associated with bioturbation, freeze–thaw cycles, etc. (Darwin, 1881; Shaler, 1891; Eyles and Ho, 1970; Fleming and Johnson, 1975; Matsuoka and Moriwaki, 1992; Auzet and Ambroise, 1996; Harris et al., 1997; Matsuoka, 1998; Gabet, 2000; Anderson, 2002; Reichman and Seabloom, 2002; Meysman et al., 2006; Wilkinson et al., 2009; Covey et al., 2010; Astete et al., 2015), and numerous authors have applied some form of a diffusion equation to describe transport and mixing of soil particles (Cousins et al., 1999; Furbish et al., 2009b, 2018a, b; Covey et al., 2010; Johnson et al., 2014; Astete et al., 2015; Campforts et al., 2016; Gray, 2018).

We emphasize that Eq. (1) is basically an advection–diffusion equation. As written, it is purely kinematic, as nothing is specified mechanically about the velocity w_p(z,t) or the diffusivity κ_z(z,t). In this view, the ideas of particle advection and diffusion are purely probabilistic constructs based on the first and second moments of the particle displacements r (Furbish et al., 2016, 2018a), as in Eqs. (2) and (3). As a description of the time evolution of the probability density f_z(z,t) of particle positions z, advection and diffusion in Eq. (1) refer to fluxes of probability. This means that, for a great number of particles within the soil element XYh, each particle “carries” a small, finite amount of probability with it as it moves over z. Moreover, despite the fact that Eq. (1) has the continuous form of a continuum advection–diffusion equation, Eq. (1) does not necessarily imply a continuum behavior. Only if conditions satisfy the continuum hypothesis can Eq. (1) be reinterpreted as an ordinary advection–diffusion equation describing transport and mixing in an individual (continuum) realization. For rarefied conditions, however, Eq. (1) represents the ensemble-expected behavior, not necessarily what happens in an individual realization (Appendix A). We elaborate on this point below in relation to expected particle positions z, ¹⁰Be concentrations and OSL ages.

We reemphasize a point made above: currently it is not possible to directly measure particle displacements r and the associated probability density $f_{\mathrm{r}}(r;z,t)$ in the setting of a natural soil. Nor is it possible to directly calculate the activity probability a(z,t). Thus, in the absence of a mechanical theory to describe these displacements, indirect measures of particle mixing behavior as reflected by profiles of ¹⁰Be concentrations and particle OSL ages are particularly valuable. Namely, any kinematic formulation of particle motions and mixing, specifically the underlying assumptions of the formulation, must be judged by its consistency with these profiles. In this vein, assuming Gaussian mixing is parsimonious, as an initial step, and in the absence of evidence of non-Gaussian behavior (Furbish et al., 2018a). This is essentially the same strategy adopted in early statistical mechanics: the veracity of the fundamental assumption of equally probable microstates (Gibbs, 1902) only can be “tested” against experimental outcomes (Tolman, 1938). Moreover, we suggest that a Gaussian formulation of mixing possesses the right granularity to accommodate uncertainty that goes with field sampling of soils. That is, this formulation captures the essence of particle mixing behavior that can be tested within the current capabilities of field-based sampling and measurements of ¹⁰Be concentrations and OSL ages.

3.2 Expected ¹⁰Be concentrations

3.2.1 Conservation of ¹⁰Be atoms

Let $f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)$ denote the joint probability density function of particle volumes V_p, ¹⁰Be atom number concentrations n_p and positions z. For particles with a given volume V_p and concentration n_p, and momentarily neglecting the production and decay of ¹⁰Be atoms, the density $f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)$ satisfies a Fokker–Planck equation of the form

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial t}}\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}\right(z,t)f_{V_{\mathrm{p}},n_{\mathrm{p}},x}\left(V_{\mathrm{p}},n_{\mathrm{p}},x,t\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial z}}].\end{array}$$

We now define a conditional joint probability density function of volumes V_p and concentrations n_p, namely,

$$\begin{array}{}\text{(5)}& f_{V_{\mathrm{p}},n_{\mathrm{p}}|z}\left(V_{\mathrm{p}},n_{\mathrm{p}}|z,t\right)={\displaystyle \frac{f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{f_z(z,t)}}.\end{array}$$

Multiplying both sides of Eq. (5) by the product V_pn_p, rearranging, and integrating with respect to V_p and n_p,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_z(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}}|z}\left(V_{\mathrm{p}},n_{\mathrm{p}}|z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z,t}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}.\end{array}$$

Note that the product V_pn_p is equal to the number of ¹⁰Be atoms within a particle of volume V_p.

The double integral on the left side of Eq. (6) defines the expected value of the product V_pn_p, that is, $\stackrel{\mathrm{‾}}{V_{\mathrm{p}}{n}_{\mathrm{p}}}(z,t)$. Thus,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_z(z,t)\stackrel{\mathrm{‾}}{V_{\mathrm{p}}{n}_{\mathrm{p}}}(z,t)\\ \text{(7)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z,t}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}.\end{array}$$

If, however, V_p and n_p are independent, then $\stackrel{\mathrm{‾}}{V_{\mathrm{p}}{n}_{\mathrm{p}}}(z,t)={\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)$. More formally, if the particles are small and within a limited size range, we may assume that V_p and n_p are independent. In this case, $f_{V_{\mathrm{p}},n_{\mathrm{p}}|z}(V_{\mathrm{p}},n_{\mathrm{p}}|z,t)=f_{V_{\mathrm{p}}|z}(V_{\mathrm{p}}|z,t)f_{n_{\mathrm{p}}|z}\left(n_{\mathrm{p}}\right|z,t)$, and we rewrite Eq. (6) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_z(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{f}_{V_{\mathrm{p}}|z}\left(V_{\mathrm{p}}|z,t\right)\mathrm{d}{V}_{\mathrm{p}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{n}_{\mathrm{p}}{f}_{n_{\mathrm{p}}|z}\left(n_{\mathrm{p}}|z,t\right)\mathrm{d}{n}_{\mathrm{p}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z,t}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}.\end{array}$$

Evaluating the integrals on the left side of Eq. (8) then yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\\ \text{(9)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z,t}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t)$ is the expected (average) particle volume and ${\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)$ is the expected particle ¹⁰Be concentration. We use these results momentarily.

We now multiply Eq. (4) by the product V_pn_p and integrate with respect to V_p and n_p, namely,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial t}}\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\\ {\displaystyle }& {\displaystyle }=-\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}\right(z,t)f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\\ \text{(10)}& {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial z}}]\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}.\end{array}$$

Noting that the random variables V_p and n_p are not functions of time t or position z, and using Leibniz's rule, Eq. (10) may be written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\right]\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_z\right(z,t)\\ \text{(11)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\partial }{\partial z}}\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z,t}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\right]).\end{array}$$

Using Eq. (7), this becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t)\stackrel{\mathrm{‾}}{V_{\mathrm{p}}{n}_{\mathrm{p}}}(z,t)\right]=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)f_z(z,t)\stackrel{\mathrm{‾}}{V_{\mathrm{p}}{n}_{\mathrm{p}}}(z,t)\right]\\ \text{(12)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[f_z(z,t)\stackrel{\mathrm{‾}}{V_{\mathrm{p}}{n}_{\mathrm{p}}}(z,t)\right]\right),\end{array}$$

and using Eq. (9),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\right]\\ \text{(13)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\right]\right).\end{array}$$

We now turn to the production and decay terms to be added to Eqs. (12) or (13).

3.2.2 Production and decay of ¹⁰Be atoms

In the absence of advection and diffusion, the joint probability density $f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)$ satisfies a statement of conservation of probability having the form of an advection equation with respect to the n_p domain, namely,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial t}}\\ \text{(14)}& {\displaystyle }& {\displaystyle }=-P(z,t){\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial n_{\mathrm{p}}}},\end{array}$$

where the advective speed $P(z,t)=\mathrm{d}{n}_{\mathrm{p}}/\mathrm{d}t$ is the rate of production of ¹⁰Be atoms per unit particle volume. Multiplying Eq. (14) by the product V_pn_p and using the product rule leads to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{V}_{\mathrm{p}}{n}_{\mathrm{p}}{\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}(V_{\mathrm{p}},n_{\mathrm{p}},z,t)}{\partial t}}\\ {\displaystyle }& {\displaystyle }=-P(z,t)V_{\mathrm{p}}{n}_{\mathrm{p}}{\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}(V_{\mathrm{p}},n_{\mathrm{p}},z,t)}{\partial n_{\mathrm{p}}}}\\ {\displaystyle }& {\displaystyle }=-P(z,t){\displaystyle \frac{\partial }{\partial n_{\mathrm{p}}}}\left[V_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}(V_{\mathrm{p}},n_{\mathrm{p}},z,t)\right]\\ \text{(15)}& {\displaystyle }& {\displaystyle }+P(z,t)V_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}(V_{\mathrm{p}},n_{\mathrm{p}},z,t).\end{array}$$

Because the product $V_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}(V_{\mathrm{p}},n_{\mathrm{p}},z,t)$ represents a proportion of all ¹⁰Be atoms in the soil column, we may at this point add the effect of radioactive decay, so that Eq. (15) becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{V}_{\mathrm{p}}{n}_{\mathrm{p}}{\displaystyle \frac{\partial f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)}{\partial t}}\\ {\displaystyle }& {\displaystyle }=-P(z,t){\displaystyle \frac{\partial }{\partial n_{\mathrm{p}}}}\left[V_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\right]\\ {\displaystyle }& {\displaystyle }+P(z,t)V_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\\ \text{(16)}& {\displaystyle }& {\displaystyle }-\mathit{\lambda }{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right),\end{array}$$

where λ denotes the decay constant.

In turn, integrating Eq. (16) with respect to V_p and n_p, and using Eq. (9),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\right]\\ {\displaystyle }& {\displaystyle }=-P(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial }{\partial n_{\mathrm{p}}}}\left[V_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\right]\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\\ {\displaystyle }& {\displaystyle }+P(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\\ \text{(17)}& {\displaystyle }& {\displaystyle }-\mathit{\lambda }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{V}_{\mathrm{p}}\mathrm{d}{n}_{\mathrm{p}}.\end{array}$$

Assuming that $f_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},\mathrm{\infty },z,t\right)=\mathrm{0}$, the first double integral on the right side of Eq. (17) is equal to zero. We then write Eq. (17) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\right]\\ {\displaystyle }& {\displaystyle }=P(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}\mathrm{d}{V}_{\mathrm{p}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{n}_{\mathrm{p}}\\ {\displaystyle }& {\displaystyle }-\mathit{\lambda }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{V}_{\mathrm{p}}\mathrm{d}{n}_{\mathrm{p}}\\ {\displaystyle }& {\displaystyle }=P(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{f}_{V_{\mathrm{p}},z}\left(V_{\mathrm{p}},z,t\right)\mathrm{d}{V}_{\mathrm{p}}\\ \text{(18)}& {\displaystyle }& {\displaystyle }-\mathit{\lambda }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{V}_{\mathrm{p}}{n}_{\mathrm{p}}{f}_{V_{\mathrm{p}},n_{\mathrm{p}},z}\left(V_{\mathrm{p}},n_{\mathrm{p}},z,t\right)\mathrm{d}{V}_{\mathrm{p}}\mathrm{d}{n}_{\mathrm{p}}.\end{array}$$

Using $f_{V_{\mathrm{p}},z}\left(V_{\mathrm{p}},z,t\right)=f_z(z,t)f_{V_{\mathrm{p}}|z}\left(V_{\mathrm{p}}|z,t\right)$ then leads to the result that

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)\right]\\ \text{(19)}& {\displaystyle }& {\displaystyle }=P(z,t)f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t)-\mathit{\lambda }{f}_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t).\end{array}$$

Thus, the production and decay terms to be added to Eqs. (12) or (13) are given by the right side of Eq. (19).

3.3 Expected particle OSL ages

3.3.1 Conservation of OSL age

In principle, the experimentally determined OSL burial age of a particle is independent of its size. In addition, as previously mentioned, quartz particles targeted for single-grain OSL analysis have a relatively narrow range of sizes (0.35–0.425 mm). For these reasons we may neglect particle volume in the following formulation.

Let $f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)$ denote the joint probability density function of particle OSL ages A_p and positions z. For particles with a given age A_p, and momentarily neglecting the production of age, the density $f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)$ satisfies a Fokker–Planck equation of the form

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)}{\partial t}}\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}\right(z,t)f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)\\ \text{(20)}& {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)}{\partial z}}].\end{array}$$

We now define a conditional joint probability density function of ages A_p, namely,

$$\begin{array}{}\text{(21)}& f_{A_{\mathrm{p}}|z}\left(A_{\mathrm{p}}|z,t\right)={\displaystyle \frac{f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)}{f_z(z,t)}}.\end{array}$$

With Eqs. (20) and (21) in place, we multiply both by A_p, integrate with respect to A_p, then follow the same steps as presented in Sect. 3.2.1 above to give

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)\right]=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)f_z(z,t){\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)\right]\\ \text{(22)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)\right]\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{A_{\mathrm{p}}}(z,t)$ is the expected particle OSL age. We now turn to the production term to be added to Eq. (22).

3.3.2 Production of OSL age

In the absence of advection and diffusion, the joint probability density $f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)$ satisfies a statement of conservation of probability having the form of an advection equation with respect to the A_p domain, namely,

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\partial f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)}{\partial t}}=-S{\displaystyle \frac{\partial f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)}{\partial A_{\mathrm{p}}}},\end{array}$$

where the advective speed $S=\mathrm{d}{A}_{\mathrm{p}}/\mathrm{d}t=\mathrm{1}$ is the rate at which particles accumulate OSL age. We then multiple Eq. (23) by A_p, integrate with respect to A_p, then follow the same steps as presented in Sect. 3.2.2 above to give

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial }{\partial t}}\left[f_z(z,t){\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)\right]=S{f}_z(z,t).\end{array}$$

Thus, the production term to be added to Eq. (22) is given by the right side of Eq. (24). We elaborate below in practical terms on the relation between the rate S and the radiation dose rate during particle motions within the soil.

3.3.3 Variance of OSL ages

Because of its significance for sampling of particles for OSL analysis, here we consider the variance of particle OSL ages. Let $f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)$ denote the joint probability density function of ages A_p and positions z. We now form the conditional probability density function,

$$\begin{array}{}\text{(25)}& f_{A_{\mathrm{p}}|z}\left(A_{\mathrm{p}}|z,t\right)={\displaystyle \frac{f_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)}{f_z(z,t)}}.\end{array}$$

Multiplying by ${\left(A_{\mathrm{p}}-{\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}\right)}^{\mathrm{2}}$, rearranging and integrating with respect to A_p,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_z(z,t)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(A_{\mathrm{p}}-{\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}\right)}^{\mathrm{2}}{f}_{A_{\mathrm{p}}|z}\left(A_{\mathrm{p}}|z,t\right)\mathrm{d}{A}_{\mathrm{p}}\\ \text{(26)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(A_{\mathrm{p}}-{\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}\right)}^{\mathrm{2}}{f}_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)\mathrm{d}{A}_{\mathrm{p}}.\end{array}$$

This yields

$$\begin{array}{}\text{(27)}& f_z(z,t)m_{\mathrm{2}}(z,t)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(A_{\mathrm{p}}-{\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}\right)}^{\mathrm{2}}{f}_{A_{\mathrm{p}},z}\left(A_{\mathrm{p}},z,t\right)\mathrm{d}{A}_{\mathrm{p}},\end{array}$$

where m₂(z,t) denotes the variance of particle OSL ages.

In turn, multiplying Eq. (20) by ${\left(A_{\mathrm{p}}-{\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}\right)}^{\mathrm{2}}$ and integrating with respect to A_p – recognizing that ${\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}$ is a function of position and time and therefore judiciously applying the product rule and Leibniz's rule – then leads to the conclusion that

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[f_z\right(z,t\left)m_{\mathrm{2}}\right(z,t\left)\right]=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)f_z(z,t)m_{\mathrm{2}}(z,t)\right]\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[f_z\right(z,t\left)m_{\mathrm{2}}\right(z,t\left)\right]\right)\\ \text{(28)}& {\displaystyle }& {\displaystyle }+\mathrm{2}{\mathit{\kappa }}_z(z,t)f_z(z,t){\left[{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)}{\partial z}}\right]}^{\mathrm{2}}.\end{array}$$

Thus the variance m₂(z,t) satisfies a Fokker–Planck equation with a source-like term involving the average age ${\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)$. Because this term depends on the structure of ${\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)$, it therefore is indirectly associated with the production of OSL age. However, it is straightforward to show that direct production of the variance m₂(z,t) of OSL ages is zero.

3.4 Advection and diffusion

The Fokker–Planck equation is basically an advection–diffusion equation. But here we reemphasize that the ¹⁰Be concentration n_p and the OSL age A_p are intensive properties of individual particles, and the volume V_p is an extensive property of individual particles. These quantities do not experience advection and diffusion as normally envisioned as occurring in a continuum. To be clear, the particles experience advection and diffusion, and the quantities n_p, V_p and A_p are merely carried with the particles.

With respect to a soil column with dimensions XYh, let us assume a great number N of particles. Then the product $N{f}_z(z,t)=c(z,t)$ may be interpreted as the expected number concentration of particles. That is, c(z,t)XYdz represents the expected number of particles within the volume XYdz between z and z+dz at time t. We may then rewrite Eq. (1) as

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\partial c(z,t)}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)c(z,t)-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial c(z,t)}{\partial z}}\right],\end{array}$$

which looks like a familiar advection–diffusion equation.

Similarly, $N{f}_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}(z,t){\stackrel{\mathrm{‾}}{n}}_{\mathrm{p}}(z,t)=n(z,t)$ represents the expected number concentration of ¹⁰Be atoms, and $f_z(z,t){\stackrel{\mathrm{‾}}{V}}_{\mathrm{p}}$ represents the volumetric particle concentration. We may then rewrite Eq. (13) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial n(z,t)}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left[w_{\mathrm{p}}(z,t)n(z,t)-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial n(z,t)}{\partial z}}\right]\\ \text{(30)}& {\displaystyle }& {\displaystyle }+P(z,t)-\mathit{\lambda }n(z,t),\end{array}$$

where P(z,t) now is interpreted as the production rate per unit volume of soil.

We write the product $N{f}_z(z,t){\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)$ as $c(z,t){\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)$, noting that c(z,t) now specifically refers to particles with finite OSL age A_p. To simplify the notation, we denote the first moment of particle OSL ages as $m_{\mathrm{1}}(z,t)={\stackrel{\mathrm{‾}}{A}}_{\mathrm{p}}(z,t)$. We may then rewrite Eq. (22) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[c\right(z,t\left)m_{\mathrm{1}}\right(z,t\left)\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left(w_{\mathrm{p}}\right(z,t\left)c\right(z,t\left)m_{\mathrm{1}}\right(z,t)\\ \text{(31)}& {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[c\right(z,t\left)m_{\mathrm{1}}\right(z,t\left)\right])+Sc(z,t).\end{array}$$

For the variance m₂(z,t),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[c(z,t)m_{\mathrm{2}}(z,t)\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left(w_{\mathrm{p}}\right(z,t\left)c\right(z,t\left)m_{\mathrm{2}}\right(z,t)\\ {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[c\right(z,t\left)m_{\mathrm{2}}\right(z,t\left)\right])\\ \text{(32)}& {\displaystyle }& {\displaystyle }+\mathrm{2}{\mathit{\kappa }}_z(z,t)c(z,t){\left[{\displaystyle \frac{\partial m_{\mathrm{1}}(z,t)}{\partial z}}\right]}^{\mathrm{2}}.\end{array}$$

Figure 2 Schematic diagram of (a) soil element with dimensions XYh. Particles move from the soil–saprolite interface (z=0) into the element at a steady rate W and are eroded from the surface (z=h). Particles experience a mean motion (gray arrows) with superimposed mixing motions. (b) in situ ¹⁰Be production rate P(z). (c) idealized luminescence dose rate D as the sum of the external rate D_e(z) and the contribution from cosmic rays D_c(z). Compare with Fig. 1 in Mudd and Yoo (2010).

4 The steady one-dimensional problem

We now turn to a benchmark situation inspired by the pioneering work of Lal (1991) and Lal and Chen (2005) concerning CRN profiles within rock, and within well-mixed soils above rock, undergoing steady surface erosion. With reference to Fig. 2, we imagine the idealized situation involving a one-dimensional vertical mean motion of particles through a soil column, where steady surface erosion plus any chemical mass losses matches the rate of soil production at the base of the column (e.g., Mudd and Yoo, 2010; Dixon and Riebe, 2014; Granger and Riebe, 2014). Although idealized, given that surface erosion rates generally are not steady (e.g., Small et al., 1997; Parker and Perg, 2005; Schaller et al., 2009), this benchmark nonetheless represents a valuable starting point for assessing actual conditions in field settings, including the possibility of a sudden change in surface erosion (Granger and Riebe, 2014), and as a contrast for two-dimensional transport by soil creep (Small et al., 1999; Anderson, 2015; Furbish et al., 2018b). With respect to cosmogenic nuclides – ¹⁰Be in particular – previous formulations of this problem have focused on two end-member cases: absence of soil particle mixing, and the so-called “well-mixed” case (or “complete” mixing) (e.g., Lal and Chen, 2005; Granger and Riebe, 2014), without reference to partial mixing or to the possible significance of the vertical structure of mixing, that is, whether particle mixing is uniform or depth dependent. This contrasts with the idea that soil disturbances and associated mixing likely involve a systematic depth dependence (Humphreys and Field, 1998; Cousins et al., 1999; Roering, 2004; Wilkinson et al., 2009). No analogous benchmark formulation exists for particle OSL ages.

We note that quartz enrichment (Small et al., 1999; Granger and Riebe, 2014) due to chemical weathering and mass loss may occur during any transient approach to steady conditions, but under steady conditions this enrichment does not impact the mechanical transport and mixing of quartz particles. In addition, we are for simplicity neglecting the vertical variation in soil bulk density that can occur with bioturbation (e.g., Furbish et al., 2009b; see Fig. 4 therein).

In this steady problem, note that $w_{\mathrm{p}}(z,t)=W$ and ${\mathit{\kappa }}_z(z,t)={\mathit{\kappa }}_z\left(z\right)$. We consider two forms of the particle diffusivity κ_z(z). In the first case we consider uniform mixing such that κ_z(z)→K_z. In the second case we consider a linear variation in mixing such that ${\mathit{\kappa }}_z\left(z\right)=K_{z}z/h$. This represents the first-order structure of a depth dependency in mixing which, although currently not well constrained, appeals to the idea that disturbances leading to particle mixing systematically decline with depth (Humphreys and Field, 1998; Cousins et al., 1999; Roering, 2004; Wilkinson and Humphreys, 2005; Wilkinson et al., 2009; Johnson et al., 2014). These two cases provide a straightforward contrast for considering how the form of κ_z(z) might influence the profiles of ¹⁰Be concentration and particle OSL age. Following Furbish et al. (2018a, b), we define a Péclet number as $\mathit{Pe}=Wh/K_z$. This provides a measure of the overall intensity of mixing. A large value of Pe represents weak mixing, whereas a small value of Pe represents strong mixing.

Following Furbish et al. (2018b), we assume that particles experience a constant radiation dose rate D (Fig. 2) during their motions within the soil column. Indeed, single-grain OSL systematics require assuming a constant natural dose rate in order to calculate a burial age A_p from the measured particle luminescence and a regeneration curve created by subjecting the particle to varying experimental “equivalent dose” values (Duller, 2008). But the natural dose rate that a particle experiences may vary with its position, and therefore with time, as the particle moves up and down within the soil column. This means that a particle will yield a luminescence signal, and thus an OSL age, that depends on its history of exposure to different dose rates, but this particle history cannot be inferred in the experimental determination of its OSL age.

With respect to the source $S=\mathrm{d}A/\mathrm{d}t=\mathrm{1}$ of particle OSL aging in Eq. (31), we are essentially assuming, as described above, that particles experience a uniform radiation dose rate during their motions within the soil column. Namely, assuming homogeneous soil material and moisture content, the external dose rate D_e(z) supplied by the radioactive decay of elements within the surrounding soil is uniform below ∼30 cm (or less according to Aitken et al., 1985) and declines toward the soil surface because of the incomplete gamma dose field at shallow depths (Fig. 2C). The dose rate D_c(z) due to cosmic rays (varying with latitude and altitude) declines nonlinearly below the soil surface (Prescott and Hutton, 1988, 1994). The total dose rate D(z) equals the sum of the external and cosmic rates. In general, the cosmic contribution tends to offset the decline of the external rate. But this depends on the relative magnitudes of these two contributions, where the magnitude of the external rate is determined by the mineral content of the soil, and the associated concentration of radioactive elements.

If the magnitude of the cosmic dose rate is similar to that of the external dose rate near the soil surface, then the total dose rate is approximately uniform (Fig. 2c). If, however, the cosmic rate does not fully offset the decrease in the external rate, we nonetheless suggest that the assumption of a uniform dose rate is a reasonable starting point for comparison with deviations in OSL age profiles that might be expected from a nonuniform dose rate, particularly under conditions of moderate to strong particle mixing, whose effects likely mask spatial variations in the total dose rate (e.g., Furbish et al., 2018b). That is, this is a parsimonious assumption – that the effects of mixing of ages outweigh any consequence of a nonuniform dose field. Previous studies using luminescence to examine soil mixing show relatively uniform total dose rates (e.g., Heimsath et al., 2002; Johnson et al., 2014).

In order to present our results below in a manner that highlights the effects of differences in the intensity and depth dependence of particle mixing, it is convenient to define the following dimensionless quantities denoted by circumflexes:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{z}={\displaystyle \frac{z}{h}},\widehat{c}\left(\widehat{z}\right)={\displaystyle \frac{c\left(z\right)}{c\left(h\right)}},\widehat{n}\left(\widehat{z}\right)={\displaystyle \frac{n\left(z\right)}{n\left(h\right)}}\\ \text{(33)}& {\displaystyle }& {\displaystyle }\mathrm{and}{\widehat{m}}_j\left(\widehat{z}\right)={\left({\displaystyle \frac{W}{h}}\right)}^j{m}_j\left(z\right).\end{array}$$

Here, $\widehat{z}$ denotes the dimensionless height within the soil column above the soil–saprolite interface, $\widehat{n}\left(\widehat{z}\right)$ denotes the dimensionless concentration of ¹⁰Be atoms relative to the concentration at the soil surface, $\widehat{c}\left(\widehat{z}\right)$ denotes the dimensionless number concentration of particles with finite OSL ages relative to the concentration at the soil surface, and ${\widehat{m}}_j\left(\widehat{z}\right)$ denotes the jth moment ($j=\mathrm{1},\mathrm{2}$) of OSL ages relative to the mean residence time, h∕W, of target quartz particles.

The analytical results presented in the next two sections involving $\widehat{n}\left(\widehat{z}\right)$ and ${\widehat{m}}_j\left(\widehat{z}\right)$ are derived in the appendixes of this paper. As described therein, each of the statements of conservation above must satisfy specific boundary conditions that depend on uniform versus nonuniform particle mixing. Here are key constraints. The ¹⁰Be flux across the soil surface equals the flux into the soil column across the soil–saprolite interface plus the total production of ¹⁰Be within the column. The flux of particles with finite OSL age across any surface normal to z is zero, and the concentration of these particles at the surface is equal to the concentration of OSL-sensitive particles entering the base of the column, although these take on finite OSL ages only after they reach the surface and are bleached. The expected OSL age at the soil surface is zero, and a diffusive flux of age across the surface matches the total production of age within the column. Particles with finite OSL age cannot be imported to the soil column. We defer commenting on the results presented next until we present our numerical simulations in Sect. 5.

4.1 Expected ¹⁰Be concentrations

Assuming that ¹⁰Be production is due to spallation (e.g., Gosse and Phillips, 2001), the production rate P(z) in Eq. (30) is (Lal, 1991; Small et al., 1999)

$$\begin{array}{}\text{(34)}& P\left(z\right)=P_{\mathrm{0}}{e}^{-(h-z)/l_{\mathrm{s}}},\end{array}$$

where P₀ is the ¹⁰Be production rate at the surface and l_s is the e-folding attenuation length of the soil. Neglecting the decay of ¹⁰Be with its half-life of ∼10⁶ years, then for steady conditions Eq. (30) becomes

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[Wn\left(z\right)-{\mathit{\kappa }}_z\left(z\right){\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}z}}\right]=P_{\mathrm{0}}{e}^{-(h-z)/l_{\mathrm{s}}}.\end{array}$$

For uniform mixing with κ_z(z)→K_z, the solution of Eq. (35) is (Appendix B)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{n}\left(\widehat{z}\right)=e^{-\mathit{Pe}(\mathrm{1}-\widehat{z})}\\ \text{(36)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{Pe}}_{l_{\mathrm{s}}}}{{\mathit{Pe}}_{l_{\mathrm{s}}}-\mathrm{1}}}\left[e^{-h(\mathrm{1}-\widehat{z})/l_{\mathrm{s}}}-e^{-\mathit{Pe}(\mathrm{1}-\widehat{z})}\right],\end{array}$$

where ${\mathit{Pe}}_{l_{\mathrm{s}}}=W{l}_{\mathrm{s}}/K_z$ is a secondary Péclet number. In turn, for nonuniform mixing with ${\mathit{\kappa }}_z\left(z\right)=K_{z}z/h$, the solution of Eq. (35) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{n}\left(\widehat{z}\right)={\widehat{z}}^{\mathit{Pe}}+\mathit{Pe}[{\left(-{\displaystyle \frac{h\widehat{z}}{l_{\mathrm{s}}}}\right)}^{\mathit{Pe}}\mathrm{\Gamma }\left(-\mathit{Pe},-{\displaystyle \frac{h\widehat{z}}{l_{\mathrm{s}}}}\right)\\ \text{(37)}& {\displaystyle }& {\displaystyle }-{\left(-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right)}^{\mathit{Pe}}\mathrm{\Gamma }\left(-\mathit{Pe},-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right){\widehat{z}}^{\mathit{Pe}}]e^{-h/l_{\mathrm{s}}},\end{array}$$

where Γ denotes the incomplete gamma function. Note that Eq. (37) has real and imaginary parts. Only the real part is physically meaningful in this problem.

4.2 Expected particle OSL ages

Recall that the number concentration c(z,t) in Eq. (31) specifically refers to particles with finite OSL age. For steady conditions Eq. (1) becomes

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[Wc\left(z\right)-{\mathit{\kappa }}_z\left(z\right){\displaystyle \frac{\mathrm{d}c\left(z\right)}{\mathrm{d}z}}\right]=\mathrm{0}.\end{array}$$

Using this, Eq. (31) is simplified to

$$\begin{array}{}\text{(39)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[-{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}\right]=Sc\left(z\right).\end{array}$$

For the variance m₂(z),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[-{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\displaystyle \frac{\mathrm{d}{m}_{\mathrm{2}}\left(z\right)}{\mathrm{d}z}}\right]\\ \text{(40)}& {\displaystyle }& {\displaystyle }-\mathrm{2}{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\left[{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}\right]}^{\mathrm{2}}=\mathrm{0}.\end{array}$$

Note that Eqs. (39) and (40) involve only diffusion, not advection. Advection and diffusion of particles possessing finite OSL ages involve the transport and mixing of OSL ages, thus influencing the age moments. But because upward advection of these particles is balanced by downward diffusion under steady conditions, this balance sets the OSL age structure wherein diffusion maintains the steady, finite values of the age moments in the presence of production of OSL age.

For uniform mixing with κ_z(z)→K_z, the solution of Eq. (38) is (Appendix C)

$$\begin{array}{}\text{(41)}& \widehat{c}\left(\widehat{z}\right)=e^{-\mathit{Pe}(\mathrm{1}-\widehat{z})}.\end{array}$$

For nonuniform mixing with ${\mathit{\kappa }}_z\left(z\right)=K_{z}z/h$, the solution of Eq. (38) is

$$\begin{array}{}\text{(42)}& \widehat{c}\left(\widehat{z}\right)={\widehat{z}}^{\mathit{Pe}}.\end{array}$$

In turn, using these results for $\widehat{c}\left(\widehat{z}\right)$, for uniform mixing the solution of Eq. (39) is (Appendix D)

$$\begin{array}{}\text{(43)}& {\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)=S(\mathrm{1}-\widehat{z})+{\displaystyle \frac{S{e}^{-\mathit{Pe}}}{\mathit{Pe}}}\left[\mathrm{1}-e^{\mathit{Pe}(\mathrm{1}-\widehat{z})}\right],\end{array}$$

and for nonuniform mixing the solution of Eq. (39) is (Appendix D)

$$\begin{array}{}\text{(44)}& {\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)={\displaystyle \frac{S\mathit{Pe}}{\mathrm{1}+\mathit{Pe}}}(\mathrm{1}-\widehat{z}).\end{array}$$

For the variance ${\widehat{m}}_{\mathrm{2}}$ the solution of Eq. (40) for uniform mixing is (Appendix E)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}}{\mathit{Pe}}}(\mathrm{1}-\widehat{z})\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{4}{S}^{\mathrm{2}}}{{\mathit{Pe}}^{\mathrm{2}}}}\left[(\mathrm{1}+\mathit{Pe})e^{-\mathit{Pe}}-(\mathrm{1}+\mathit{Pe}\widehat{z})e^{-\mathit{Pe}\widehat{z}}\right]\\ \text{(45)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{S^{\mathrm{2}}}{{\mathit{Pe}}^{\mathrm{2}}}}\left(e^{-\mathrm{2}\mathit{Pe}}-e^{-\mathrm{2}\mathit{Pe}\widehat{z}}\right),\end{array}$$

and for nonuniform mixing the solution of Eq. (40) is (Appendix E)

$$\begin{array}{}\text{(46)}& {\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)={\displaystyle \frac{S^{\mathrm{2}}{\mathit{Pe}}^{\mathrm{2}}}{(\mathrm{2}+\mathit{Pe})(\mathrm{1}+\mathit{Pe}{)}^{\mathrm{2}}}}\left(\mathrm{1}-{\widehat{z}}^{\mathrm{2}}\right).\end{array}$$

Also for reference below, the column-averaged particle OSL age ${\widehat{M}}_{\mathrm{1}}$ within the soil is

$$\begin{array}{}\text{(47)}& {\widehat{M}}_{\mathrm{1}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\widehat{c}\left(\widehat{z}\right){\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)\mathrm{d}\widehat{z}.\end{array}$$

For uniform mixing, Eqs. (41) and (43) lead to

$$\begin{array}{}\text{(48)}& {\widehat{M}}_{\mathrm{1}}={\displaystyle \frac{S}{{\mathit{Pe}}^{\mathrm{2}}}}\left(\mathrm{1}-e^{-\mathrm{2}\mathit{Pe}}\right)-{\displaystyle \frac{\mathrm{2}S{e}^{-\mathit{Pe}}}{\mathit{Pe}}}.\end{array}$$

For nonuniform mixing, Eqs. (42) and (44) lead to

$$\begin{array}{}\text{(49)}& {\widehat{M}}_{\mathrm{1}}={\displaystyle \frac{S\mathit{Pe}}{\mathrm{1}+\mathit{Pe}}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathit{Pe}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}+\mathit{Pe}}}\right).\end{array}$$

We comment further on the results above after presenting our numerical simulations.

5 Numerical simulations

We now turn to numerical simulations of particles undergoing random-walk motions within the soil column, during which they accumulate ¹⁰Be atoms within the production field and undergo OSL “aging” following their most recent encounters with the soil surface. These simulations have two purposes.

First, the random-walk motions implied by the probabilistic formulations above are in principle straightforward to implement numerically, and it is important to demonstrate that such computational results match the analytical results presented. In doing this, the simulations reveal important information that is not readily apparent in the analytical results. This includes an illustration of the variability in ¹⁰Be concentrations and OSL ages of individual particles, in contrast to expected values at positions z, with important implications for interpreting field-based measurements, and the nature of the terms in Eqs. (14) and (23) describing production of ¹⁰Be atoms and particle OSL age.

Second, numerical simulations of particle motions within soils offer important opportunities to examine phenomena that cannot readily be treated analytically, for example, effects of particle residence times on mineral weathering or effects of a nonuniform radiation dose rate. So, spinning our first objective around, any numerical simulation of random-walk motions must be able to correctly reproduce benchmark (analytical) solutions before being applied to more complex situations, for example, two-dimensional motions and unsteady conditions. The simulations presented here highlight important aspects involved.

Following Furbish et al. (2018a, b), we adopt a straightforward Eulerian–Lagrangian algorithm to simulate particle motions in a mass-conserving manner. Particles are numerically introduced to the base of the soil column (z=0), then undergo a mean upward motion equal to W with superimposed Gaussian fluctuations. For uniform mixing the particle diffusivity is set as κ_z=K_z, and the random walk becomes

$$\begin{array}{}\text{(50)}& z(t+\mathrm{\Delta }t)=z\left(t\right)+W\mathrm{\Delta }t+R_z\left(a\right),\end{array}$$

where Δt denotes the time step and R_z(a) is a Gaussian random variable with argument $a=(\mathrm{2}{K}_z\mathrm{\Delta }t{)}^{\mathrm{1}/\mathrm{2}}$. For nonuniform mixing with $\mathit{\kappa }\left(z\right)=K_{z}z/h$, the random walk becomes

$$\begin{array}{}\text{(51)}& z(t+\mathrm{\Delta }t)=z\left(t\right)+W\mathrm{\Delta }t+R_z\left(a\right)+{\mathit{\kappa }}_z^{\prime }\mathrm{\Delta }t,\end{array}$$

with argument $a=\left[\mathrm{2}{\mathit{\kappa }}_z\right(z+\mathrm{0.5}{\mathit{\kappa }}_z^{\prime }\mathrm{\Delta }t)\mathrm{\Delta }t{]}^{\mathrm{1}/\mathrm{2}}$, where ${\mathit{\kappa }}_z^{\prime }=\partial {\mathit{\kappa }}_z\left(z\right)/\partial z$. This yields a mass-conserving behavior, that is, one that prevents particles from unrealistically drifting from sites with high particle diffusivity to sites with low diffusivity. Moreover, this algorithm has been shown to work for variations in diffusivity that are not linear (e.g., Legg and Raupach, 1982; Hunter et al., 1993; Visser, 1997). The theoretical basis of Eq. (51) and its relation to the Fokker–Planck equation are covered in these references and in Appendix G of Furbish et al. (2018a).

Each particle accumulates ¹⁰Be atoms as a function of its local position z, and it accumulates a numerical OSL age from the time of its last encounter with the soil surface. We spin up each simulation to a steady-state condition, where the rate at which particles exit the soil column is equal to the rate at which they are introduced at the base, and particles within the column are distributed uniformly over the thickness h. The total spin-up time involves at least four e-folding residence times h∕W. At steady state, the total number of particles within the column is $N_T\approx N_{\mathrm{c}}(h/W)/\mathrm{\Delta }t$, where N_c is the number of particles in the cohort introduced at each time step. We use a minimum of N_T≈ 10 000 for the ¹⁰Be simulations.

The lower boundary (z=0) is treated as a reflecting boundary. For each particle reaching the upper boundary (z=h), it either may leave the column with a specified probability that ensures global particle conservation, or it is reflected. In the case of particle OSL ages, the numerical age of an individual particle is set to zero if it is reflected at z=h. The effect of this is to correctly mimic the boundary condition in the formulation above: that ${\widehat{m}}_j=\mathrm{0}$ at $\widehat{z}=\mathrm{1}$. In actuality, however, bleaching of particles can occur just below the soil surface with light penetration (to a few particle diameters) and with heating from fires at the surface (Wilkinson and Humphreys, 2005; Duller, 2008), such that actual values ${\widehat{m}}_j=\mathrm{0}$ occur below the soil surface.

All simulated N_T particles at steady-state possess a ¹⁰Be value. But only a proportion of these N_T particles possess finite OSL ages at steady state, as not all of them reach the surface to subsequently take on finite OSL ages. We cannot know this proportion a priori. Thus, it is important to insist on global particle conservation in the simulations, involving verification of a specified N_T together with a uniform distribution of particle positions z. In addition, we increase N_T (up to 20 000) and the total spin-up time (up to six residence times h∕W) for the OSL simulations to ensure that a sufficiently large number of particles is included in our calculations of expected values. However, this is not entirely possible with large Péclet number Pe, as described below.

Figure 3 Plot of dimensionless ¹⁰Be concentration $\widehat{n}=n\left(z\right)/n\left(h\right)$ versus dimensionless height $\widehat{z}=z/h$ showing simulated particle concentrations ${\widehat{n}}_{\mathrm{p}}$ (gray dots) for Pe=100, 10, 1, 0.1, and estimates of expected concentrations $\widehat{n}$ averaged within 0.1h intervals (black circles) with one standard deviation bars. Simulations represent uniform mixing with κ_z=K_z. Right solid line is the theoretical result, and left solid line represents the absence of mixing.

Figure 4 Plot of dimensionless ¹⁰Be concentration $\widehat{n}=n\left(z\right)/n\left(h\right)$ versus dimensionless height $\widehat{z}=z/h$ showing simulated particle concentrations ${\widehat{n}}_{\mathrm{p}}$ (gray dots) for $\mathit{Pe}=\mathrm{100},\mathrm{10},\mathrm{1},\mathrm{0.1}$, and estimates of expected concentrations $\widehat{n}$ averaged within 0.1h intervals (black circles) with one standard deviation bars. Simulations represent nonuniform mixing with κ_z=K_z z∕h. Right solid line is the theoretical result, and left solid line represents the absence of mixing.

5.1 ¹⁰Be concentrations

The simulated, expected ¹⁰Be concentrations closely match the theoretical results for different values of the Péclet number Pe involving both uniform mixing (Fig. 3) and nonuniform mixing (Fig. 4). These profiles show that with weak mixing (large Pe), the expected concentration approaches the original exponential solution provided by Lal (1991). With strong mixing (small Pe), the expected concentration becomes increasingly uniform over the soil column, approaching the concentration at the soil surface. With uniform mixing (Fig. 3), the concentration $\widehat{n}\left(\mathrm{0}\right)$ may be finite, as diffusion effectively moves particles downward to the soil–saprolite interface. With nonuniform mixing (Fig. 4), the concentration $\widehat{n}\left(\mathrm{0}\right)$ is anchored by the value within the saprolite, as diffusion weakens downward then vanishes at the soil–saprolite interface.

Figure 5 Example histograms representing the distribution $f_{{\widehat{n}}_{\mathrm{p}}}({\widehat{n}}_{\mathrm{p}},\widehat{z})$ using simulated values of ${\widehat{n}}_{\mathrm{p}}$ from Fig. 4 (Pe=1) over the intervals (a) $\mathrm{0.9}\widehat{z}-\mathrm{1.0}\widehat{z}$, (b) $\mathrm{0.5}\widehat{z}-\mathrm{0.6}\widehat{z}$ and (c) $\mathrm{0.1}\widehat{z}-\mathrm{0.2}\widehat{z}$. Analogous histograms associated with uniform mixing show a similar structure.

With both uniform and nonuniform mixing, the distribution $f_{{\widehat{n}}_{\mathrm{p}}}({\widehat{n}}_{\mathrm{p}},\widehat{z})$ of ¹⁰Be concentrations ${\widehat{n}}_{\mathrm{p}}$ of individual particles within any small interval $\mathrm{d}\widehat{z}$ systematically varies with vertical position and the Péclet number Pe (Fig. 5). Notably, this distribution at any $\widehat{z}$ is approximately symmetrical about the expected value for large Pe and becomes increasingly skewed with decreasing Pe. The expected concentration $\widehat{n}\left(\widehat{z}\right)$ at small Pe thus is strongly influenced by the tail of this distribution, that is, by particles possessing concentrations much larger than the modal concentration.

Figure 6 Plot of dimensionless OSL age ${\widehat{A}}_{\mathrm{p}}=(W/h)A_{\mathrm{p}}$ versus dimensionless height $\widehat{z}=z/h$ showing simulated particle ages ${\widehat{A}}_{\mathrm{p}}$ (gray dots) for Pe=100, 10, 1, 0.1, and estimates of expected values ${\widehat{m}}_{\mathrm{1}}$ averaged within 0.1h intervals (black circles) with one standard deviation bars. Simulations represent uniform mixing with κ_z=K_z. Solid line is the theoretical result.

Figure 7 Plot of dimensionless OSL age ${\widehat{A}}_{\mathrm{p}}=(W/h)A_{\mathrm{p}}$ versus dimensionless height $\widehat{z}=z/h$ showing simulated particle ages ${\widehat{A}}_{\mathrm{p}}$ (gray dots) for Pe=100, 10, 1, 0.1, and estimates of expected values ${\widehat{m}}_{\mathrm{1}}$ averaged within 0.1h intervals (black circles) with one standard deviation bars. Simulations represent nonuniform mixing with κ_z=K_z z∕h. Solid line is the theoretical result.

5.2 Particle OSL ages

The simulated, expected OSL ages closely match the theoretical results for different values of the Péclet number Pe involving both uniform mixing (Fig. 6) and nonuniform mixing (Fig. 7), where we note that the simulations yield meaningful results only near the surface for large Péclet number Pe. (Because the concentration $\widehat{c}\left(\widehat{z}\right)$ of particles with finite OSL ages declines rapidly with depth for large Pe (Appendix C), achieving reasonable numerical values of the expected age ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)$ over the entire soil thickness would require unreasonably large computational memory and time.) These profiles show that, with weak mixing (large Pe), the expected particle OSL age increases linearly, or approximately linearly, with depth. With strong mixing (small Pe), the expected age becomes increasingly uniform and close to zero over the soil column. With uniform mixing, the diffusive flux of age must vanish at the soil–saprolite interface, so with finite diffusivity K_z, the slope $\mathrm{d}{\widehat{m}}_{\mathrm{1}}/\mathrm{d}\widehat{z}{|}_{\widehat{z}=\mathrm{0}}=\mathrm{0}$. With nonuniform mixing, the diffusive flux of age likewise vanishes at the soil–saprolite interface as the diffusivity goes to zero. But the magnitude of the slope $\mathrm{d}{\widehat{m}}_{\mathrm{1}}/\mathrm{d}\widehat{z}$ is finite near this interface in order to compensate the decreasing diffusivity.

Figure 8 Example histograms representing the distribution $f_{{\widehat{A}}_{\mathrm{p}}}({\widehat{A}}_{\mathrm{p}},\widehat{z})$ using simulated values of ${\widehat{A}}_{\mathrm{p}}$ from Fig. 7 (Pe=1) over the intervals (a) $\mathrm{0.9}\widehat{z}-\mathrm{1.0}\widehat{z}$, (b) $\mathrm{0.5}\widehat{z}-\mathrm{0.6}\widehat{z}$ and (c) $\mathrm{0.1}\widehat{z}-\mathrm{0.2}\widehat{z}$. Analogous histograms associated with uniform mixing show a similar structure.

With both uniform and nonuniform mixing, the distribution $f_{{\widehat{A}}_{\mathrm{p}}}({\widehat{A}}_{\mathrm{p}},\widehat{z})$ of particle OSL ages within any small interval $\mathrm{d}\widehat{z}$ mostly is highly skewed (Fig. 8). This skew increases with decreasing Péclet number Pe. Particularly with nonuniform mixing, the expected OSL age ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)$ thus is strongly influenced by the tail of this distribution, that is, by particles possessing finite ages much larger than the modal age.

Figure 9 Plot of dimensionless variance ${\widehat{m}}_{\mathrm{2}}=(W/h{)}^{\mathrm{2}}{m}_{\mathrm{2}}$ versus dimensionless height $\widehat{z}=z/h$ showing values obtained from simulations (Pe=1) for uniform mixing (black circles) and nonuniform mixing (gray circles) compared with theoretical values (black and gray lines).

The simulated second moment ${\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)$ of OSL ages reasonably matches the theoretical results for different values of the Péclet number Pe for both uniform and nonuniform mixing. Focusing on the example of Pe=1 (Fig. 9), the variance ${\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)$ rapidly increases with depth from zero at the soil surface, then becomes relatively uniform with increasing depth. With both uniform and nonuniform mixing, the variance at any position $\widehat{z}$ generally decreases with decreasing Pe (Figs. 6 and 7). We note that, whereas in any individual simulation the numerical estimates of the expected OSL ages ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)$ closely match the theoretical values with large N_T for small Pe (Figs. 6 and 7) – a consequence of the central limit theorem – numerical estimates of the variance ${\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)$ may fluctuate about the theoretical values from one simulation to the next (Fig. 9).

Figure 10 Exceedance probability plots of dimensionless particle OSL age ${\widehat{A}}_{\mathrm{p}}=(W/h)A_{\mathrm{p}}$ for (a) uniform mixing and (b) nonuniform mixing for Péclet numbers Pe=10, 1 and 0.1.

Figure 11 Plot of dimensionless column-averaged OSL age ${\widehat{M}}_{\mathrm{1}}=(W/h)M_{\mathrm{1}}$ versus Péclet number $\mathit{Pe}=Wh/K_z$ for uniform and nonuniform mixing.

The simulations suggest that particle OSL ages within the entire soil column are distributed approximately exponentially for both uniform and nonuniform mixing (Fig. 10), where the column-averaged age ${\widehat{M}}_{\mathrm{1}}$ varies systematically with the Péclet number Pe. Interestingly, based on Eqs. (48) and (49), the average ${\widehat{M}}_{\mathrm{1}}$ increases from zero at Pe→0, reaches a maximum of ${\widehat{M}}_{\mathrm{1}}\sim \mathrm{0.1}$ near Pe∼1, then declines again with increasing Pe (Fig. 11), consistent with the simulations (Fig. 10). For Pe→0, small values of ${\widehat{M}}_{\mathrm{1}}$ reflect the idealized condition of complete mixing, where particles that reach the soil surface and are bleached and then move downward rather than being eroded, nonetheless, frequently return to the soil surface due to strong mixing. For large Pe, small values of ${\widehat{M}}_{\mathrm{1}}$ reflect that particles with finite OSL age tend to remain near the soil surface due to the strong effect of upward advection, and thus frequently return to it, many exiting by erosion before accumulating large ages. Relatively large values of ${\widehat{M}}_{\mathrm{1}}$ at intermediate Pe reflect the effects of an approximate balance between upward advection and downward diffusion of particles with finite OSL age such that particles return to the soil surface less frequently. We emphasize that the maximum value of ${\widehat{M}}_{\mathrm{1}}$ is a fraction of the mean residence time h∕W.

6 Discussion and conclusions

6.1 Implications of rarefied transport conditions

We emphasize that, in contrast to continuum formulations of advection and diffusion of material (e.g., mass) measured as an intensive quantity (e.g., concentration) of the continuum, the extensive and intensive particle properties V_p, n_p and A_p “belong” to the particles, not to the bulk soil. For this reason, a formulation of advection and diffusion of ¹⁰Be concentrations and expected particle OSL ages based on the Fokker–Plank equation provides a satisfactory way to parse the behavior of the particle-centric quantities V_p, n_p and A_p. In the case of ¹⁰Be, the formulation describes the behavior of the expected value of individual particle concentrations at a position z. When this is combined with the expected particle volume and number concentration, the expected ¹⁰Be concentration n(z,t) then may be considered an intensive property of the soil at position z. As a consequence, the expected concentration n(z,t) satisfies what looks like an ordinary advection–diffusion equation with production and decay terms – although this does not necessarily imply a continuum behavior (Sect. 3.1).

In the case of particle OSL ages, the formulation similarly describes the behavior of the expected value (and the variance) of individual particle OSL ages at a position z. By definition, our interest is in this expected particle OSL age, as this is what is determined from single-grain OSL measurements. It therefore does not make sense to define OSL age as an intensive property of the soil by combining the expected particle OSL age with the expected particle number concentration (resulting in a “total” OSL age at position z). Moreover, by maintaining this distinction, the formulation reveals that the expected particle OSL age (as well as the variance) satisfies a diffusion-like equation according to Eqs. (39) and (40), not an advection–diffusion equation. This is in contrast to the idea that the “age” of a fluid parcel moving through a continuum domain satisfies an advection–diffusion equation with a production term equal to unity, as described in oceanographic and hydrological applications (England, 1995; Goode, 1996). This is important because, unlike a continuum material, the expected number concentration c(z,t) of particles possessing a finite OSL age generally is not uniform over z (Appendix C). That is, this concentration does not mimic a uniform continuum domain within which particle OSL age is transported.

An essential lesson is this: when the quantity of interest can be expressed as a total value within an interval dz, as with the total number of ¹⁰Be atoms, then this quantity may be treated as an intensive property of the bulk soil. When the quantity of interest is an expected value within dz, as with the moments m_j(z) of particle OSL age, then this quantity cannot be expressed as an intensive property of the bulk soil, and its behavior must be coupled with that of the expected concentration c(z) of the particles possessing the property. Similar quantities include, for example, particle size (in relation to descriptions of vertical sorting; Campforts et al., 2016) and particle age as measured from the time of entry into the mechanically active soil column (in relation to studies of particle weathering; White and Brantley, 2003; Mudd and Furbish, 2006; Almond et al., 2007; Anderson et al., 2007; Yoo and Mudd, 2008; Mudd and Yoo, 2010; Ferrier et al., 2016). In contrast, there is a growing interest in the use and interpretation of the total OSL intensity of bulk soil samples as measured by portable OSL readers (Muñoz-Salinas et al., 2010; Sanderson and Murphy, 2010; Stang et al., 2012; Munyikwa and Brown, 2014; Gray et al., 2017; Gray, 2018; Porat et al., 2018). The luminescence intensities of individual particles – decidedly a random variable (Gray, 2018) – contribute to the total measured intensity. Thus, because the quantity of interest is the total intensity rather than expected moments of individual particle intensities, the total intensity can be formulated as being an intensive property of the bulk soil (Gray et al., 2017; Gray, 2018).

Throughout we have emphasized that ¹⁰Be concentrations and OSL ages are to be considered expected values. Moreover, this expectation is defined with respect to an interval z to z+dz in a soil element with finite areal dimension XY, and it formally is an ensemble average, rather than the expected value associated with an individual realization. The significance of this bears on the practical issue of sampling soil material for measurements of ¹⁰Be and particle OSL age, in view of the fact that disturbance-driven particle motions in soils are patchy and intermittent at many scales, where most particles are at rest most of the time. Namely, vertical profiles of soil properties measured in an individual soil pit (where XY is on the order of 1×1 m) reflect a “snapshot” of possible conditions (Furbish et al., 2009b). This snapshot represents the recent history of transport and mixing, one that is much shorter than the typical soil particle residence time, W∕h.

We cannot avoid this issue of legacy (or “inheritance”), namely, the likelihood that what is being measured reflects only the recent history of transport and mixing as opposed to conditions consistent with an imagined behavior averaged over longer timescales, as represented by the expected profiles in Figs. 3, 4, 6 and 7 above. In the case of measured profiles of ¹⁰Be concentrations and particle OSL ages, this has two parts. Consider a profile that reflects an expected steady-state condition (Figs. 3, 4, 6 or 7). Disturbances that contribute to the mixing motions consistent with the profile may occur at different length scales and with different frequencies, where large disturbances may involve coherent motions whose effects are more akin to stirring than mixing, thus momentarily producing irregularities about the expected state. To the extent that mixing is adequately characterized as being diffusive, then we may define a relaxation timescale as $T=r^{\mathrm{2}}/\mathit{\kappa }\sim r^{\mathrm{2}}/\langle r^{\mathrm{2}}\rangle f$, where now the mixing motion r is used as a measure of the length scale of disturbance, f denotes a characteristic frequency of disturbance, and the angle brackets denote ensemble averaging. With r² in the numerator and 〈r²〉 in the denominator, this expression highlights the duality of disturbances: these provide mixing motions, yet this mixing is responsible for diffusive smoothing of disturbance produced irregularities about the expected profile state. This is in marked contrast to, say, classic molecular diffusion, where molecular motions smooth irregularities, but are not the source of disturbances to the expected state. Thus, for a given ensemble-averaged disturbance magnitude 〈r²〉^1∕2, the relaxation time T goes with the square of the scale of disturbance and inversely with the characteristic frequency of disturbances. For a given frequency f, effects of big disturbances tend to persist whereas effects of small disturbances do not. In either case, this persistence decreases with increasing disturbance frequency f (i.e., decreasing Péclet number Pe).

We now take the ensemble average of relaxation timescales T over all disturbance length scales, namely, $\langle T\rangle \sim \mathrm{1}/f$. This indicates that the overall relaxation in response to a range of disturbance scales goes simply with the reciprocal of the disturbance frequency f. Thus, regardless of the mixture of disturbance scales involved, the disturbance frequency has a dominant role in setting the relaxation timescale. Then, for example, if disturbances and mixing motions are consistently small and relatively uniform in comparison to the size of the soil pit (and the size of individual soil samples), and if the frequency of the disturbances is sufficiently high, then one might anticipate observing at any instant only small variations about the expected steady-state profile. If, however, disturbances are infrequent and patchy at the scale of the soil pit or larger, then one might anticipate a greater likelihood of observing conditions unlike the expected profile. Conversely, frequent and spatially uniform large disturbances likely would lead to wholesale homogenization of tracer particles.

This points to the need to avoid over-interpreting the forms of profiles from individual soil pits in terms of what these forms might reflect about the vertical structure of mixing (e.g., uniform versus depth-dependent mixing). Unfortunately, this issue is exacerbated by the reality that digging soil pits and sampling for ¹⁰Be concentrations and particle OSL ages is quite laborious, and subsequent analytical analyses are prohibitively expensive. In addition, in choosing soil pit sites, we often avoid sites with evidence of recent disturbance. On the one hand, this strategy may obviate the sampling of conditions that likely deviate from averaged conditions, but on the other hand, it neglects observing profile irregularities that reflect the full range of disturbance scales. Connecting sampling strategies (e.g., involving multiple soil pits, choosing sampling intervals within individual pits) with appropriate averaging relative to scales of disturbances and mixing remains an important open question.

Momentarily assuming that mixing conditions are reasonably reflected by the expected particle OSL age profile ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)$, then the results above bear on the practical question of variability in these expected ages as a consequence of small sample sizes. As a point of reference, Heimsath et al. (2002) sampled an average of 41 quartz grains from each of one to three vertical positions within four soil pits. Of the total 10 samples, on average 19 grains had finite OSL ages. Johnson et al. (2014) analyzed 42–49 grains from each of five intervals in a single soil pit. Considering only grains with finite OSL ages, the sample size N_s from each vertical interval is about 20–50 in these examples. Regardless of the form of the distribution of finite particle OSL ages with variance σ² within each interval (Fig. 8), the central limit theorem suggests that the standard error s_e of the estimate of the mean is $s_{\mathrm{e}}\approx \mathit{\sigma }/\sqrt{N_{\mathrm{s}}}$, or, in dimensionless form, ${\widehat{s}}_{\mathrm{e}}\approx \widehat{\mathit{\sigma }}/\sqrt{N_{\mathrm{s}}}$.

Let is assume that within a small interval of $\widehat{z}$, ${\widehat{\mathit{\sigma }}}^{\mathrm{2}}={\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)$ from Eqs. (45) or (46). We may then write

$$\begin{array}{}\text{(52)}& {\widehat{s}}_{\mathrm{e}}\left(\widehat{z}\right)\approx \pm \sqrt{{\displaystyle \frac{{\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)}{N_{\mathrm{s}}}}}.\end{array}$$

This yields an estimate of ${\widehat{s}}_{\mathrm{e}}\left(\widehat{z}\right)$ depending on the intensity and structure of mixing in relation to the sample size N_s, and it represents uncertainty in the mean value ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)$ that is in addition to analytical uncertainty associated with single-grain OSL age estimates. The standard errors ${\widehat{s}}_{\mathrm{e}}\left(\widehat{z}\right)$ for uniform and nonuniform mixing are similar, although nonuniform mixing generally yields smaller values of ${\widehat{s}}_{\mathrm{e}}\left(\widehat{z}\right)$. The well-known formula Eq. (52) suggests that, obtaining a standard error ${\widehat{s}}_{\mathrm{e}}\left(\widehat{z}\right)$ of specified magnitude within a small interval at position $\widehat{z}$ requires that $N_{\mathrm{s}}\approx {\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)/s_{\mathrm{e}}^{\mathrm{2}}$. Because ${\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)$ increases with depth (Fig. 9), uncertainty in the estimate of the expected value ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)$ increases with depth for a given sample size N_s, as directly reflected in the data of Heimsath et al. (2002) and Johnson et al. (2014). Stated another way, there may be value in judiciously varying N_s with depth when faced with a research budget that limits the total number of single-grain OSL age analyses. We note, however, that this uncertainty associated with sample size cannot be distinguished from effects of any legacy of disturbances as described above.

The results of the numerical simulations as depicted in Figs. 3, 4, 6 and 7 provide an important perspective on the nature of production of ¹⁰Be and OSL age in relation to particle transport and mixing, and the associated structuring of the profiles n(z) and m₁(z). We note that the points in Figs. 3 and 4 represent large samples drawn from the joint probability density $f_{n_{\mathrm{p}},z}(n_{\mathrm{p}},z,t)$, and the points in Figs. 6 and 7 represent samples drawn from the joint probability density $f_{A_{\mathrm{p}},z}(A_{\mathrm{p}},z,t)$. With respect to $f_{n_{\mathrm{p}},z}(n_{\mathrm{p}},z,t)$, at any instant a particle within the n_p−z domain only can move in the positive n_p direction due to its accumulation of ¹⁰Be atoms (neglecting decay). Similarly, with respect to $f_{A_{\mathrm{p}},z}(A_{\mathrm{p}},z,t)$, a particle within the A_p−z domain only can move in the positive A_p direction due to its accumulation of OSL age. This means that the distribution $f_{n_{\mathrm{p}}}(n_{\mathrm{p}},z)$ or $f_{A_{\mathrm{p}}}(A_{\mathrm{p}},z)$ at any position z as depicted in Figs. 5 and 8 is at all instants being uniformly advected in the positive n_p or A_p direction. In both cases, particles at any instant may move in either the positive or negative z direction due to their random-walk motions.

Combining Eqs. (4) and (14), neglecting particle volume and the decay of ¹⁰Be, and assuming steady conditions,

$$\begin{array}{}\text{(53)}& -{\displaystyle \frac{\partial }{\partial z}}(q_A+q_D)-P\left(z\right){\displaystyle \frac{\partial f_{n_{\mathrm{p}},z}(n_{\mathrm{p}},z)}{\partial n_{\mathrm{p}}}}=\mathrm{0},\end{array}$$

where q_A and q_D denote the advective and diffusive parts of the flux. Similarly, combining Eqs. (20) and (23),

$$\begin{array}{}\text{(54)}& -{\displaystyle \frac{\partial }{\partial z}}(q_A+q_D)-S{\displaystyle \frac{\partial f_{A_{\mathrm{p}},z}(A_{\mathrm{p}},z)}{\partial A_{\mathrm{p}}}}=\mathrm{0}.\end{array}$$

These highlight how production at any position within the n_p−z or A_p−z domain is exactly balanced by the local, combined effects of particle advection and diffusion. Consider the density $f_{n_{\mathrm{p}},z}(n_{\mathrm{p}},z)$. With reference to Fig. 5, at all locations (n_p,z) where the derivative $\partial f_{n_{\mathrm{p}},z}(n_{\mathrm{p}},z)/\partial n_{\mathrm{p}}<\mathrm{0}$, the effect of production is to increase the ¹⁰Be content at these locations in proportion to the production rate P(z) and the magnitude of this derivative; at locations where the derivative $\partial f_{n_{\mathrm{p}},z}(n_{\mathrm{p}},z)/\partial n_{\mathrm{p}}>\mathrm{0}$ the effect of production is to decrease the ¹⁰Be content at these locations. The variation in q_A and q_D with respect to z must be such that their combined divergence balances these effects of production. In turn, consider the density $f_{A_{\mathrm{p}},z}(A_{\mathrm{p}},z)$. With reference to Fig. 8, at all locations (A_p,z) where the derivative $\partial f_{A_{\mathrm{p}},z}(A_{\mathrm{p}},z)/\partial A_{\mathrm{p}}<\mathrm{0}$ the effect of particle aging is to increase the OSL age content at these locations in proportion to the magnitude of this derivative; at locations where the derivative $\partial f_{A_{\mathrm{p}},z}(A_{\mathrm{p}},z)/\partial A_{\mathrm{p}}>\mathrm{0}$ the effect of particle aging is to decrease the OSL age content at these locations. Variations in q_A and q_D with respect to z must then compensate these effects.

We normally envision that local production of a quantity implies a local increase in the quantity. But this is not necessarily so when viewed in the n_p−z or A_p−z domain. Only when the production is averaged via integration over the n_p or A_p domain, as in Sect. 3.2.2 and 3.3.2, does a production term emerge as normally envisioned. This point further highlights a key idea underlying the formulation: extensive and intensive particle properties are not in themselves subject to advection and diffusion, but rather are merely carried with the particles as these undergo advection and diffusion with respect to z. Indeed, the production terms in Eqs. (53) and (54) represent only advection over the n_p and A_p domains, not diffusion (mixing) over these domains.

The numerical simulations suggest that the overall particle OSL age distribution is approximately exponential (Fig. 10), consistent with field data (see data of Heimsath et al., 2002 as described by Furbish et al., 2018b). This result awaits a theoretical explanation. Meanwhile, as described by Furbish et al. (2018b), the distribution $f_{T_{\mathrm{r}}}\left(T_{\mathrm{r}}\right)$ of the return times T_r between successive encounters of a particle with the soil surface is expected to be a power-law distribution with an undefined mean (Redner, 2001) for the idealized situation involving uniform Gaussian mixing in a vertically unbounded domain, in the absence of upward advection. Because the OSL age of a particle increases at the same rate as its (eventual) return time, the distribution of OSL ages also is likely to be a power-law distribution in this situation. However, upward advection (with surface erosion) combined with a finite soil thickness has the effect of strongly tempering this distribution, yielding an approximate exponential form. Further tempering is provided with nonuniform mixing, where diffusion decreases with depth then vanishes at the soil–saprolite interface. This behavior of particle OSL ages is entirely analogous to the exponential tempering of the power-law distribution of residence times of particles undergoing burial and exhumation in a stream channel, where a finite sediment thickness limits the depth of burial. At long times the particles fully explore the accessible thickness, and a finite (unchanging) average residence time emerges (Voepel et al., 2013).

The emergence of a maximum average OSL age ${\widehat{M}}_{\mathrm{1}}$ at an intermediate Péclet number Pe∼1 (Fig. 11) is in direct contrast with the two-dimensional case involving downslope transport by creep without surface erosion (Furbish et al., 2018b), where the average OSL age monotonically decreases with decreasing Pe. In this case, at large Pe, OSL particles remain near the surface (as in the one-dimensional case), but they can accumulate large ages before exiting the soil mantle downslope. Moreover, in the one-dimensional case, that the average OSL age is a fraction of the mean particle residence time lends support to the idea of defining two distinct populations of OSL tracers (Heimsath et al., 2002; Furbish et al., 2018b) – those with finite age and those that are saturated – having an “infinite” age, inasmuch as the mean residence time is much smaller than the determinable OSL age limit (Murray and Olley, 2002).

That the numerical simulations mimic analytical solutions for the benchmark situation of a one-dimensional mean motion involving both uniform and nonuniform mixing with varying mixing intensities lends confidence in applying the numerics to more complicated situations. Such situations might be motivated by questions concerning consequences of transient conditions of surface erosion and soil production, aeolian inputs to the soil, particle weathering in relation to particle aging, accumulation of luminescence signals with nonuniform dose rates, and the structuring of tracer particles under depositional conditions. Our experience suggests the need to implement the numerics of boundary conditions carefully, ensuring consistency with global particle conservation.

Here we return to our starting point. Our use of the Fokker–Planck equation assumes Gaussian diffusion of tracer particles. As described above, this is a parsimonious choice whose consequences, and veracity, must be judged by its consistency with measurable outcomes of mixing, including profiles of CRN concentrations and OSL ages as emphasized here, but possibly to include other soil properties. We suggest that a Gaussian model of particle mixing is robust inasmuch as this mixing behavior is insensitive to the form of the probability distribution of particle displacements, f_r(r), so long as this distribution is not heavy-tailed. We further emphasize that the effective particle diffusivity may actually represent motions involving a mixture of characteristic length scales and associated frequencies of occurrence in settings involving both biotic and abiotic disturbances. We also acknowledge that it may be more appropriate to consider some disturbances, for example, macro-disturbances by tree throw and fossorial animals, as having the effect of stirring rather than mixing, where homogenization occurs at length scales comparable to the mechanically active soil thickness (see next section). This points to the need for a clearer understanding of the spatiotemporal structure of mixing motions in adopting more sophisticated (i.e., non-Gaussian) models of mixing behavior. The goal is to understand the information content of tracers aimed at constraining mechanical formulations of transport and mixing, notably in relation to soil creep. The one-dimensional benchmark situation described here is a key starting point due to the lessons it offers.

Figure 12 Plot of dimensionless expected concentration $\widehat{n}\left(\widehat{z}\right)=n\left(z\right)/n\left(h\right)$ of ¹⁰Be atoms versus dimensionless height $\widehat{z}=z/h$ with uniform mixing (solid lines) and nonuniform mixing (dashed lines) as these vary with the Péclet number $\mathit{Pe}=Wh/K_z$.

Figure 13 Plot of dimensionless expected particle OSL age ${\widehat{A}}_{\mathrm{p}}\left(\widehat{z}\right)=(W/h)A_{\mathrm{p}}\left(z\right)$ versus dimensionless height $\widehat{z}=z/h$ with uniform mixing (solid lines) and nonuniform mixing (dashed lines) as these vary with the Péclet number $\mathit{Pe}=Wh/K_z$.

6.2 Assessing the intensity and depth dependence of mixing

Here we focus on results for the one-dimensional benchmark case (Sect. 4, Fig. 2) – specifically the profiles of expected ¹⁰Be concentrations and particle OSL ages – to suggest constraints on assessing the intensity and depth dependence of mixing. For ease of comparison, we collect these profiles from Figs. 3 and 4, and from Figs. 6 and 7, and combine them in Figs. 12 and 13.

As described above, these profiles systematically vary with the Péclet number, $\mathit{Pe}=Wh/K_z$. In the case of ¹⁰Be concentrations, the profile converges to the exponential solution provided by Lal (1991) for weak mixing (large Pe), and it converges to a uniform value equal to the surface concentration for strong mixing (small Pe). In the case of expected particle OSL ages, the profiles vary approximately linearly with depth, and converge to a uniform value close to zero for strong mixing.

Not surprisingly, with weak mixing the ¹⁰Be and OSL profiles for uniform and nonuniform mixing are virtually indistinguishable (Figs. 12 and 13), as the profiles in this case are mostly determined by the mean motion. Similarly, with strong mixing the ¹⁰Be profiles are not markedly different except near the base of the soil column, and the OSL age profiles are nearly the same. Significant differences in the profiles appear only in the presence of intermediate mixing intensities. The essence of these differences at intermediate intensities (Pe∼1) arises from how rapidly the particle diffusivity decreases with increasing depth (Sect. 5.1 and 5.2). Thus, the forms of the profiles might change in detail in the presence of a more complicated (e.g., nonlinear) mixing structure. Nonetheless, these results suggest that ¹⁰Be and OSL age profiles may help constrain the mixing structure in the presence of intermediate mixing intensities, albeit depending on the resolution of measurements.

These profiles highlight that uniform particle mixing is not synonymous with the idea of complete mixing, and why a uniform profile of ¹⁰Be concentration or particle OSL age does not necessarily indicate the presence of uniform mixing. The idea of “uniform mixing” refers to the mixing structure wherein the statistical qualities of particle random walks are independent of vertical position. In contrast, “complete mixing” refers to an idea from reservoir theory (Bolin and Rodhe, 1973) that particle mixing within a specified control volume is sufficiently thorough that the probability of a particle exiting the volume is independent of its residence time in the volume (Bolin and Rodhe, 1973; Furbish et al., 2018a). This probability in fact is strongly conditioned by the geometry of particle motions, specifically the proximity of the inflow and outflow locations relative to the particle trajectories, and the degree of mixing between these locations (Bolin and Rodhe, 1973). Both uniform and nonuniform mixing yield uniform ¹⁰Be and OSL profiles in the limit of Pe→0. That said, complete particle mixing within soils is mechanically unlikely, a point that is consistent with available ¹⁰Be and OSL data concerning creeping soils (Furbish et al., 2018a, b), and deserving reexamination in interpreting ¹⁰Be profiles with respect to surface ages and denudation rates (Schaller et al., 2009). This point also is consistent with the idea of depth-dependent mixing (Humphreys and Field, 1998; Cousins et al., 1999; Roering, 2004; Wilkinson and Humphreys, 2005; Wilkinson et al., 2009; Johnson et al., 2014; Gray, 2018), in which the local intensity of mixing declines with depth.

Here we step back and look at published data. We first note that, whereas our benchmark case involves a steady one-dimensional mean motion, available field-based measurements of ¹⁰Be concentrations and OSL particle ages mostly pertain to transient conditions or involve two-dimensional downslope soil transport. One cannot make a direct comparison between tracer profiles sampled on sloping surfaces and the one-dimensional results depicted in Figs. 12 and 13. For example, the upper boundary conditions examined here are quite different from those in the two-dimensional case. One effect of these differences is directly reflected by the column-averaged OSL age as this varies non-monotonically with the Péclet number Pe (Figs. 10 and 11) versus the monotonic variation of this quantity with Pe for two-dimensional particle motions (Furbish et al., 2018b, Fig. 6 therein). Nonetheless, in comparing our results with those presented in Furbish et al. (2018b, Figs. 4 and 5 therein), it is clear that the basic forms of profiles resulting from one-dimensional and two-dimensional transport systematically vary in like manner with the intensity of mixing, as characterized by the Péclet number Pe.

As an important backdrop to the benchmark case examined here, ¹⁰Be profiles from a flight of five marine terraces near Santa Cruz, California, illustrate the continued accumulation of ¹⁰Be atoms with increasing terrace age within the mixed soil and underlying undisturbed material, under the condition of negligible surface erosion (Perg et al., 2001, Figs. 2 and 4 therein; Granger and Riebe, 2014, Fig. 9 therein). Near-surface concentrations are relatively uniform and in three cases (terraces 1, 3 and 5) decline toward the value at the base of the assumed mixing depth, suggesting Pe∼1 and likely an associated decline in mixing intensity. Concentrations are mostly centered about a vertically averaged value that is less than the surface concentration that would occur with steady surface erosion. Similarly, as noted by Furbish et al. (2018b), uniform concentrations of ¹⁰Be in weakly developed soils on the crests of moraines near Pinedale, Wyoming, suggest well-mixed conditions near the surface (Schaller et al., 2009), although there is inconsistency with expected concentrations based on the formulation of Lal and Chen (2005) for the well-mixed case; there also is uncertainty in the calculated lowering rates and mixing depths, and the sites may represent transient conditions.

Relatively uniform ¹⁰Be profiles from hillslopes in the Great Smokey Mountains reflect strongly mixed conditions at the sample locations (Jungers et al., 2009, Fig. 7 therein; reproduced in Anderson, 2015, Fig. 14 therein), likely due to effects of tree throw and other bioturbation events that stir the soil over much of its ∼60 cm thickness. Within the context of the analysis above, these conditions suggest that Pe<1. Similarly, five profiles sampled on hillslopes at Gordon Gulch, Colorado, display a mixture of conditions, varying from relatively uniform concentrations (Pe<1) to an approximately linear variation with depth (Pe>1) (Foster et al., 2015, Fig. 7 therein; reproduced in Anderson, 2015, Fig. 14 therein). In turn, three profiles measured along a 100 m catena flow line in a soil developed from granitic bedrock on Osborn Mountain, Wyoming (Small et al., 1999, Fig. 6 therein; reproduced in Anderson, 2015, Fig. 14 therein), reflect conditions consistent with Pe≳1 (Furbish et al., 2018b), where relatively uniform concentrations in the upper parts of the profiles then decrease in the lower one-third.

An OSL age profile A_p(z) based on single quartz grains collected from a bioturbated soil developed on a basalt flow on the Denna Plain in northeast Queensland, Australia (Johnson et al., 2014, Fig. 2 therein; Furbish et al., 2018b, Fig. 10 therein), most closely matches the benchmark case described here. This profile suggests an approximately linear increase in OSL ages with depth, as in Fig. 13. In addition, the sampled quartz grains likely were added to the soil at its surface, a boundary condition that is consistent with the theoretical formulation (Furbish et al., 2018b; Appendixes C and D). Moreover, the OSL ages are only a fraction of the estimated mean residence time at this site, consistent with moderate to strong mixing (Furbish et al., 2018b). Although mixing at this site likely varies with depth, the similarity between profiles in Fig. 13 suggests that the mixing structure cannot be distinguished. Similarly, the OSL age profiles A_p(z) reported by Heimsath et al. (2002, Fig. 1 therein; Furbish et al., 2018b, Fig. 8 therein) based on single grains of quartz collected from a hillslope with nonuniform soil thickness over granitic bedrock in the Nunnock River catchment, Australia, suggest an approximately linear increase in OSL ages with depth. Although involving downslope transport, the profiles are consistent with strong mixing (small Pe), possibly including macro-disturbances (Heimsath et al., 2002). Moreover, the distribution of all particle OSL ages is approximately exponential with an average age that is much smaller than the calculated mean soil residence time, consistent with strong mixing (Furbish et al., 2018b; Fig. 10).

In all cases summarized above, the profiles suggest moderate (Pe≳1) to strong (Pe<1) mixing. Distinguishing between uniform and non-uniform mixing likely will require higher-resolution sampling than reported in these cases. Our own bias is that many if not most settings with significant mixing by bioturbation or the effects of freezing and thawing likely involve depth-dependent mixing. At least for the one-dimensional case examined here, CRN profiles are capable of revealing mixing intensity and possibly mixing structure for Pe∼1. In contrast, OSL profiles are capable of revealing mixing intensity, but not likely mixing structure.

To our knowledge there are no available measurements of profiles of ¹⁰Be concentrations and OSL ages taken together. We suggest that there is merit in doing just this as a means to provide a more demanding test of formulations of transport and mixing (Furbish, 2003; Roering et al., 2004). We also reiterate that our results provide an analytical benchmark for assessing the veracity of emerging numerical methods aimed at simulating particle transport and mixing, to include Eulerian–Lagrangian descriptions of particle motions that might incorporate individual detrital grain CRN concentrations (Codilean et al., 2010) as well as fully treating the effects of a nonuniform radiation dose field. This includes simulations that start from probabilistic, physically based formulations of total luminescence intensities as measured by portable OSL readers (Gray et al., 2017; Gray, 2018), as an addition to multi-grain and single-grain analyses aimed at extracting particle burial ages. Such measurements may be capable of revealing mixing structure, as well as intensity, from relatively high-resolution sampling since particles involved in accumulating luminescence signals are likely to be more uniformly distributed within the soil column relative to those possessing finite OSL ages (Appendix C).

Appendix A: Rarefied versus continuum conditions

To further illustrate the significance of the probabilistic formulation of conservation in relation to rarefied versus continuum conditions, here we start with the familiar example of Brownian motion, the initial formal description of which is separately attributable to Einstein (1905) and von Smoluchowski (1906). With reference to Fig. A1, let x denote a coordinate along which Brownian particles take one-dimensional random walks, where x extends indefinitely in the positive and negative directions about the origin x=0. Suppose that a particle starts at the origin at time t=0, and with equal probability moves in the positive or negative direction during successive small intervals dt. By the definition of a random walk, the motion of the particle – specifically its expected position x after an interval of time t>0 – can be predicted only in a probabilistic sense. Namely, letting f_x(x,t) denote the probability density function of possible positions x, then this density satisfies a Fokker–Planck equation involving only its diffusion term:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial f_x(x,t)}{\partial t}}={\mathit{\kappa }}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}{f}_x(x,t)}{\partial x^{\mathrm{2}}}},\end{array}$$

where κ_x denotes the particle diffusivity. The solution of Eq. (A1) is the Gaussian distribution with mean μ_x=0, namely,

$$\begin{array}{}\text{(A2)}& f_x(x,t)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{4}\mathit{\pi }{\mathit{\kappa }}_{x}t}}}{e}^{-x^{\mathrm{2}}/\mathrm{4}{\mathit{\kappa }}_{x}t}.\end{array}$$

For this highly rarefied system involving a single particle, we can only offer probabilistic predictions of its position at time t. For example, we may confidently state that with probability $p=\mathrm{1}/\mathrm{2}$ the particle is either at a position x<0 or at a position x>0. Or we may state that with probability p≈0.68 the particle is within the domain defined by plus one and minus one standard deviations about the mean position, namely, $-\sqrt{\mathrm{2}{\mathit{\kappa }}_{x}t}<x<+\sqrt{\mathrm{2}{\mathit{\kappa }}_{x}t}$. For this single-particle system (realization), the actual particle position x_a is represented by a Dirac distribution $\mathit{\delta }(x_a-x,t)$ (Fig. A1), but this cannot be predicted deterministically.

Figure A1 Plot of coordinate position x of particle undergoing a random-walk motion showing Gaussian distribution f_x(x,t) of expected positions at time t as the solution, Eq. (A2), of the Fokker–Planck equation, and the actual (example) particle position x=x_a, and uniform steady-state distribution $f_x\left(x\right)=\mathrm{1}/\mathrm{2}$ for a bounded domain such that $-\mathrm{1}<x<\mathrm{1}$.

Figure A2 Histograms of particle positions x at time t for one system showing that (a) with a great number N of particles representing a continuum condition this histogram converges to the smooth Gaussian distribution in Fig. A1 as dx→0; in this example N=100 000, and (b) with a modest number of particles representing a rarefied condition this histogram is irregular and discontinuous; in this example N=200.

Let us now imagine an arbitrarily great number N of identical, independent particles that start at the origin x=0 at time t=0, each undergoing a random walk during t>0. When viewed together, the distribution of these particles at time t=0 is given by the Dirac distribution, namely, $f_x(x,\mathrm{0})=\mathit{\delta }\left(x\right)$. At any time t>0 these particles are distributed according to Eq. (A2). That is, because N is arbitrarily large, the proportion of particles within any small interval x to x+dx closely matches what is predicted by Eq. (A2), namely f_x(x,t)dx, such that in the limit of dx→0 the actual distribution of positions x varies smoothly (continuously) and converges to Eq. (A2) (Fig. A2). In contrast to the highly rarefied single-particle system in the previous example, we may thus assume that this great number of particles, occurring in one system (realization), satisfies the continuum hypothesis. Nonetheless, upon randomly selecting a single particle from this system, we still can only offer probabilistic predictions of its position at time t – as in the example above involving a system with a single particle. Moreover, note that the continuous distribution of positions x realized at time t for this one system involving a great number N of particles is identical to the distribution that would be realized upon pooling the x positions at time t associated with a great number N of independent systems, each involving a single particle.

Now select a system with a modest number N of particles such that conditions are rarefied. By this we mean that, after some time t, the actual distribution of particle positions x is at best represented by an irregular histogram that roughly appears Gaussian but is decidedly discontinuous (Fig. A2). Moreover, any realization involving N particles possesses a similar irregular form at time t, and no two are the same. In effect, each realization represents a sample of size N drawn from an imagined population represented by Eq. (A2). Also note that each realization involving N particles at time t is the same as N realizations, each involving a single particle, when viewed collectively at time t.

Let us now consider a great number N_e of independent but nominally identical systems – an ensemble – at any fixed time t, where each system contains N particles, large or small. We now wish to describe the ensemble-expected conditions. To envision this, consider any small interval x to x+dx. If N=1 as in the first example above, then f_x(x,t)dx is just the proportion of the N_e systems containing a particle within x to x+dx at time t. Note that this is identical to the result above involving an individual system containing a great number N=N_e of particles. If instead each system involves a great number N of particles, then f_x(x,t)dx simply becomes the expected proportion of the N particles within x to x+dx at time t, where the expectation is calculated over the N_e systems. And note that this outcome is identical to the proportion of N×N_e independent systems, each involving a single particle, which contain a particle within x to x+dx at time t. In either case, the expected proportion within the interval is the same. Moreover, we reach the same conclusion in considering a great number N_e of systems, each involving a modest number N of particles. Thus, when calculated over a great number of systems for all intervals dx, then in the limit of dx→0, the continuous function, Eq. (A2), is retrieved. The key points are these: first, whether N is relatively small (representing a rarefied condition) or N is large (representing a continuum condition), the ensemble-expected behavior represented by Eq. (A2) applies equally to both conditions in a probabilistic sense. Second, if N is small, then Eq. (A2) represents the ensemble-expected behavior, not the actual behavior of any one system (realization); if N is large, then the actual behavior of the system is expected to converge to the smooth ensemble behavior represented by Eq. (A2).

Figure A3 Histogram of particle positions x at time t→∞ for one system showing that with a modest number of particles representing a rarefied condition this histogram is irregular and discontinuous; in this example N=500.

To complete the picture, suppose that the x domain in Fig. A1 is bounded such that $-\mathrm{1}<x<\mathrm{1}$. Particles that reach these boundaries are “reflected” and remain within the domain, continuing their random walks. In the limit of t→∞, the probability density of particle positions x reaches a steady-state form, that is, $\partial f_x(x,t)/\partial t\to \mathrm{0}$ such that $f_x(x,t)\to f_x\left(x\right)$. In this limit, Eq. (A1) becomes ${\mathrm{d}}^{\mathrm{2}}{f}_x\left(x\right)/\mathrm{d}{x}^{\mathrm{2}}=\mathrm{0}$. Moreover, the probability flux $q_x=-{\mathit{\kappa }}_x\mathrm{d}{f}_x\left(x\right)/\mathrm{d}x=\mathrm{0}$ at all positions x, which means that $\mathrm{d}{f}_x\left(x\right)/\mathrm{d}x=\mathrm{0}$. These constraints together with the fact that the distribution f_x(x) must integrate to unity yield the result that $f_x\left(x\right)=\mathrm{1}/\mathrm{2}$ over the bounded domain (Fig. A1). That is, the expected distribution f_x(x) is uniform. As with the unsteady problem described above, a modest number N of particles representing rarefied conditions in any one realization is at best represented by an irregular histogram that roughly appears uniform but is decidedly discontinuous (Fig. A3). Moreover, at an arbitrary later time, the resulting distribution (histogram) would be just as irregular; it does not become smoother with increasing time. As above, the expected continuous steady-state distribution is retrieved when expected values are calculated over a great number N_e of systems.

To place these ideas within the context of soil tracer particles, including the practical assessment of rarefied versus continuum conditions, let us now imagine a soil column of thickness h=1 m. Suppose that we wish to be able to describe the probability density f_z(z) of tracer particle positions z as a smooth function at a resolution of, say, 1 cm. That is, we are aimed at a function that “looks” smooth in proceeding along z from each 1 cm increment to the next. Choosing a finer resolution would not make sense because then we are approaching the scale of individual particles. On the other hand, choosing a significantly larger resolution (e.g., 10 cm) would represent a loss of information at scales of possible interest.

Let us now measure the number of tracer particles within each of 100 cm³ representing a column of soil, that is, a column with horizontal dimensions $X=Y=\mathrm{1}$ cm. For illustration, let us assume that the expected number of particles within each cubic centimeter of soil is 10. The expected number of particles within the entire column therefore is N=1000. In this example the expected proportion f_z(z)dz with dz=1 cm is 0.01.

Figure A4 Plot of proportion $n\left(z\right)/N={\widehat{f}}_z\left(z\right)$ versus position z for the example of N=1000 particles and $X=Y=\mathrm{1}$ cm with expected proportions of f_z(z)dz=0.01 (black dots), then smoothed with boxcar filter with L=10 cm (black line), and proportion for N= 100 000 with $X=Y=\mathrm{10}$ cm and L=4 cm (gray line).

Now imagine that, due to randomness in the particle mixing process, the number of tracer particles n(z) varies from one increment dz=1 cm to the next about the expected average of 10. Assuming for illustration that this variability is spatially random at the centimeter scale, we may formally draw values of n(z) from a binomial distribution (or approximate this using a Poisson distribution or a normal distribution) with a mean of μ=10 to populate each increment in the soil column. The total expected number of particles remains N=1000, so we may calculate the associated proportions n(z)∕N to represent the function ${\widehat{f}}_z\left(z\right)$ (Fig. A4). Notice that ${\widehat{f}}_z\left(z\right)$ fluctuates significantly about the expected value f_z(z)dz=0.01.

Consider first the idea of describing this one realization as a continuum, that is, where we might imagine that ${\widehat{f}}_z\left(z\right)\mathrm{d}z$ is smooth. To do this, we apply a boxcar moving average of width L over the spatial series in Fig. A4. We may think of this width L as the averaging length defining a continuum physical point, as described in the text. Notice that, in this example, L must approach a significant proportion of the soil thickness h before a relatively smooth function emerges – that is, a function that we would not be uncomfortable in taking its derivative to define a gradient having physical significance with respect to particle transport. However, there may be uncertainty in the fidelity of this operation. For example, in the situation where n(z) is not uniform, then because averaging of this sort represents a low-pass filter, averaging may obscure physically meaningful variations of interest.

Consider now the idea of expanding either (or both) the horizontal dimensions X and Y of the soil column. In this situation both n(z) and N increase. Fluctuations ${\widehat{f}}_z\left(z\right)$ about the expected value f_z(z) systematically decrease with increasing N. Using a boxcar moving average, smoothness in ${\widehat{f}}_z\left(z\right)$ emerges with a smaller length L (Fig. A4) relative to the previous situation involving smaller XY and N. That is, ${\widehat{f}}_z\left(z\right)$ looks more like a continuous function at finer resolution. This is the reason for stating in the text that one hopes to sample over a sufficiently large area XY to obtain an approximately smooth distribution in any one realization. Note, however, that this smoothness is scale dependent, as it requires a sufficiently large area XY and assumes that the process of mixing is uniform over x and y. But this also means that variations with respect to x or y in the full three-dimensional probability density $f_{x,y,z}(x,y,z)$ of particle positions x, y and z cannot be resolved in the situation where the mixing process varies with x or y (e.g., two-dimensional transport and mixing).

Finally, consider the idea of taking an ensemble average, where we return to the small 1 cm³ sampling volumes as described above. If at any geometrically similar position ($x,y,z$) within a great number of nominally identical but independent systems the number of particles within the volume dxdydz varies from one system (realization) to the next, the central limit theorem guarantees that the ensemble average of this number converges to its ensemble-expected value. It immediately follows that the ensemble-expected distribution $f_{x,y,z}(x,y,z)$ is a smooth continuous function, or in one dimension, f_z(z) for specified XY is a smooth continuous function.

With respect to developments in the text, the Fokker–Planck equation describes the time evolution of the probability density f_z(z,t). The formulation does not assume either rarefied or continuum conditions. It is indifferent to these conditions, yet equally applicable to both. As described in the text, the Fokker–Planck equation is a special case of the master equation, a general statement of conservation of probability that is independent of scale, and which is the basis of more familiar statements of conservation when probability is reinterpreted in terms of, for example, mass, momentum or energy. If in an individual system (realization) the continuum hypothesis is satisfied (a condition that is independent of the probabilistic basis of the master equation or the Fokker–Planck equation), then the probabilistic formulation based on ensemble-expected conditions and its continuum counterpart are essentially one and the same. If, however, the continuum hypothesis is not satisfied, then one cannot defensibly start with a continuum equation, but instead must appeal to a probabilistic formulation of ensemble-expected conditions (in order to justify the use of continuously differentiable equations), with the proviso that any prediction of the behavior of an individual (rarefied) system is probabilistic in nature.

Conservation of expected number concentration of ¹⁰Be atoms

Assuming ¹⁰Be production is due to spallation (Gosse and Phillips, 2001), the number concentration n(z,t) of ¹⁰Be atoms satisfies a Fokker–Planck-like equation, namely,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial n(z,t)}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left[w(z,t)n(z,t)-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial n(z,t)}{\partial z}}\right]\\ \text{(B1)}& {\displaystyle }& {\displaystyle }+P_{\mathrm{0}}{e}^{-(h-z)/l_{\mathrm{s}}}-\mathit{\lambda }n(z,t),\end{array}$$

where w(z,t) denotes the ensemble-averaged particle velocity and κ_z(z,t) denotes the ensemble-expected particle diffusivity (Furbish et al., 2009b, 2018a, b). For steady conditions involving a one-dimensional mean motion, $\partial n(z,t)/\partial t=\mathrm{0}$, $n(z,t)\to n\left(z\right)$ and $w(z,t)\to W$. Assuming that the half-life of ¹⁰Be is much greater than the mean residence time, h∕W, of target quartz particles, then Eq. (B1) becomes

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[Wn\left(z\right)-{\mathit{\kappa }}_z\left(z\right){\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}z}}\right]=P_{\mathrm{0}}{e}^{-(h-z)/l_{\mathrm{s}}}.\end{array}$$

We first rewrite Eq. (B2) in terms of the flux $q_z\left(z\right)=Wn\left(z\right)-\mathrm{d}n\left(z\right)/\mathrm{d}z$, namely,

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\mathrm{d}q\left(z\right)}{\mathrm{d}z}}=P_{\mathrm{0}}{e}^{-h/l_{\mathrm{s}}}{e}^{z/l_{\mathrm{s}}}.\end{array}$$

Vertically integrating Eq. (B3) over the soil thickness,

$$\begin{array}{}\text{(B4)}& \underset{\mathrm{0}}{\overset{h}{\int }}{\displaystyle \frac{\mathrm{d}q\left(z\right)}{\mathrm{d}z}}\mathrm{d}z=P_{\mathrm{0}}{e}^{-h/l_{\mathrm{s}}}\underset{\mathrm{0}}{\overset{h}{\int }}{e}^{z/l_{\mathrm{s}}}\mathrm{d}z.\end{array}$$

With fixed boundaries, this yields

$$\begin{array}{}\text{(B5)}& q_z\left(h\right)=q_z\left(\mathrm{0}\right)+P_{\mathrm{0}}{l}_{\mathrm{s}}\left(\mathrm{1}-e^{-h/l_{\mathrm{s}}}\right).\end{array}$$

This says that the flux q_z(h) of ¹⁰Be atoms across the soil surface and removed by erosion is equal to the rate at which atoms enter the soil column at its base, q_z(0)=Wn(0), plus the total rate at which they are produced within the column. In this steady problem, no information is available regarding the vertically averaged concentration.

In this problem, the concentration n(0) at the soil–saprolite interface is obtained by solving the purely advective form of Eq. (B2) and using the boundary condition that $n(-\mathrm{\infty })=\mathrm{0}$. The result is

$$\begin{array}{}\text{(B6)}& n\left(\mathrm{0}\right)={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{W}}{e}^{-h/l_{\mathrm{s}}}.\end{array}$$

In turn, the flux q_z(0) is

$$\begin{array}{}\text{(B7)}& q_z\left(\mathrm{0}\right)=P_{\mathrm{0}}{l}_{\mathrm{s}}{e}^{-h/l_{\mathrm{s}}},\end{array}$$

and the flux q_z(h) is

$$\begin{array}{}\text{(B8)}& q_z\left(h\right)=P_{\mathrm{0}}{l}_{\mathrm{s}}.\end{array}$$

As described below, the concentration n(h) at the soil surface depends on what is assumed about the contributions to the flux q_z(h) at this surface. We now consider how the concentration profile n(z) differs with uniform versus nonuniform mixing.

B1 Uniform mixing

With uniform mixing (κ_z=K_z), we start by integrating Eq. (B3) to give

$$\begin{array}{}\text{(B9)}& q_z\left(z\right)=P_{\mathrm{0}}{l}_{\mathrm{s}}{e}^{-h/l_{\mathrm{s}}}{e}^{z/l_{\mathrm{s}}}+C_{\mathrm{1}}.\end{array}$$

The lower boundary condition, Eq. (B7), then gives C₁=0. Using this result together with $q_z\left(z\right)=Wn\left(z\right)-K_z\mathrm{d}n\left(z\right)/\mathrm{d}z$ then leads to

$$\begin{array}{}\text{(B10)}& {\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}z}}-{\displaystyle \frac{W}{K_z}}n\left(z\right)=-{\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{K_z}}{e}^{-h/l_{\mathrm{s}}}{e}^{z/l_{\mathrm{s}}}.\end{array}$$

We now simplify the notation and rewrite Eq. (B10) as

$$\begin{array}{}\text{(B11)}& {\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}z}}-An\left(z\right)=-B{e}^{z/l_{\mathrm{s}}},\end{array}$$

with

$$\begin{array}{}\text{(B12)}& A={\displaystyle \frac{W}{K_z}}\mathrm{and}B={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{K_z}}{e}^{-h/l_{\mathrm{s}}}.\end{array}$$

A general solution of Eq. (B11) is

$$\begin{array}{}\text{(B13)}& n\left(z\right)={\displaystyle \frac{B{l}_{\mathrm{s}}}{A{l}_{\mathrm{s}}-\mathrm{1}}}{e}^{z/l_{\mathrm{s}}}+C_{\mathrm{1}}{e}^{Az}.\end{array}$$

At this point we assume that the upper boundary flux is purely advective. Physically this means we are imagining that the rate of surface erosion E, being externally imposed, removes ¹⁰Be atoms at a rate En(h)=Wn(h). This is consistent with the requirement that the rate of quartz particle removal at the surface is equal to the rate at which quartz particles enter the soil at the soil–saprolite interface, independent of their ¹⁰Be concentration. Then, from Eq. (B8), $q_z\left(h\right)=P_{\mathrm{0}}{l}_{\mathrm{s}}=Wn\left(h\right)$ so that $n\left(h\right)=P_{\mathrm{0}}{l}_{\mathrm{s}}/W$. Using this boundary condition,

$$\begin{array}{}\text{(B14)}& C_{\mathrm{1}}={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{W}}{e}^{-Ah}-{\displaystyle \frac{B{l}_{\mathrm{s}}}{A{l}_{\mathrm{s}}-\mathrm{1}}}{e}^{h/l_{\mathrm{s}}}{e}^{-Ah}.\end{array}$$

Substituting this into Eq. (B13) then doing algebra yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }n\left(z\right)={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{W}}{e}^{-W(h-z)/K_z}\\ \text{(B15)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}^{\mathrm{2}}}{K_z}}{\displaystyle \frac{\mathrm{1}}{W{l}_{\mathrm{s}}/K_z-\mathrm{1}}}\left[e^{-(h-z)/l_{\mathrm{s}}}-e^{-W(h-z)/K_z}\right].\end{array}$$

With dimensionless height $\widehat{z}=z/h$, dimensionless concentration $\widehat{n}\left(\widehat{z}\right)=n\left(z\right)/n\left(h\right)$, primary Péclet number $\mathit{Pe}=Wh/K_z$ and secondary Péclet number ${\mathit{Pe}}_{l_{\mathrm{s}}}=W{l}_{\mathrm{s}}/K_z$, Eq. (B15) becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{n}\left(\widehat{z}\right)=e^{-\mathit{Pe}(\mathrm{1}-\widehat{z})}\\ \text{(B16)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{Pe}}_{l_{\mathrm{s}}}}{{\mathit{Pe}}_{l_{\mathrm{s}}}-\mathrm{1}}}\left[e^{-h(\mathrm{1}-\widehat{z})/l_{\mathrm{s}}}-e^{-\mathit{Pe}(\mathrm{1}-\widehat{z})}\right],\end{array}$$

which is Eq. (36) in the text.

B2 Nonuniform mixing

With nonuniform mixing (${\mathit{\kappa }}_z\left(z\right)=K_{z}z/h$) we start with Eq. (B9) with C₁=0 together with $q_z\left(z\right)=Wn\left(z\right)-K_z(z/h)\mathrm{d}n\left(z\right)/\mathrm{d}z$ to give

$$\begin{array}{}\text{(B17)}& {\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}z}}-{\displaystyle \frac{Wh}{K_{z}z}}n\left(z\right)=-{\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}h}{K_z}}{e}^{-h/l_{\mathrm{s}}}{\displaystyle \frac{e^{z/l_{\mathrm{s}}}}{z}}.\end{array}$$

We now simplify the notation and rewrite Eq. (B17) as

$$\begin{array}{}\text{(B18)}& {\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}z}}-{\displaystyle \frac{A}{z}}n\left(z\right)=-B{\displaystyle \frac{e^{z/l_{\mathrm{s}}}}{z}},\end{array}$$

with

$$\begin{array}{}\text{(B19)}& A={\displaystyle \frac{Wh}{K_z}}\mathrm{and}B={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}h}{K_z}}{e}^{-h/l_{\mathrm{s}}}.\end{array}$$

A general solution of Eq. (B18) is

$$\begin{array}{}\text{(B20)}& n\left(z\right)=B{\left(-{\displaystyle \frac{z}{l_{\mathrm{s}}}}\right)}^A\mathrm{\Gamma }\left(-A,-{\displaystyle \frac{z}{l_{\mathrm{s}}}}\right)+C_{\mathrm{1}}{z}^A,\end{array}$$

where Γ is the incomplete gamma function. Using the advective boundary condition $n\left(h\right)=P_{\mathrm{0}}{l}_{\mathrm{s}}/W$,

$$\begin{array}{}\text{(B21)}& C_{\mathrm{1}}={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{W}}{h}^{-A}-B{\left(-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right)}^A\mathrm{\Gamma }\left(-A,-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right)h^{-A}.\end{array}$$

Substituting this into Eq. (B20) then doing algebra yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }n\left(z\right)={\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}}{W}}{\left({\displaystyle \frac{z}{h}}\right)}^{Wh/K_z}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{P_{\mathrm{0}}{l}_{\mathrm{s}}h}{K_z}}[{\left(-{\displaystyle \frac{z}{l_{\mathrm{s}}}}\right)}^{Wh/K_z}\mathrm{\Gamma }\left(-{\displaystyle \frac{Wh}{K_z}},-{\displaystyle \frac{z}{l_{\mathrm{s}}}}\right)\\ {\displaystyle }& {\displaystyle }-{\left(-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right)}^{Wh/K_z}\\ \text{(B22)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{\Gamma }\left(-{\displaystyle \frac{Wh}{K_z}},-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right){\left({\displaystyle \frac{z}{h}}\right)}^{Wh/K_z}]e^{-h/l_{\mathrm{s}}}.\end{array}$$

With dimensionless height $\widehat{z}=z/h$, dimensionless concentration $\widehat{n}\left(\widehat{z}\right)=n\left(z\right)/n\left(h\right)$ and Péclet number $\mathit{Pe}=Wh/K_z$, Eq. (B22) becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{n}\left(\widehat{z}\right)={\widehat{z}}^{\mathit{Pe}}+\mathit{Pe}[{\left(-{\displaystyle \frac{h\widehat{z}}{l_{\mathrm{s}}}}\right)}^{\mathit{Pe}}\mathrm{\Gamma }\left(-\mathit{Pe},-{\displaystyle \frac{h\widehat{z}}{l_{\mathrm{s}}}}\right)\\ \text{(B23)}& {\displaystyle }& {\displaystyle }{\left(-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right)}^{\mathit{Pe}}\mathrm{\Gamma }\left(-\mathit{Pe},-{\displaystyle \frac{h}{l_{\mathrm{s}}}}\right){\widehat{z}}^{\mathit{Pe}}]e^{-h/l_{\mathrm{s}}},\end{array}$$

which is Eq. (37) in the text. Note that Eq. (B23) has real and imaginary parts. Only the real part is physically meaningful in this problem.

Like the results in Sect. A1 above, the concentration gradient at the soil surface, $\left[\mathrm{d}n\right(z)/\mathrm{d}z{]}_h=\mathrm{0}$, as a consequence of assuming an advective boundary condition.

Appendix C: Conservation of particles with finite OSL age

The number concentration c(z,t) of particles with finite OSL age satisfies a Fokker–Planck-like equation, namely,

$$\begin{array}{}\text{(C1)}& {\displaystyle \frac{\partial c(z,t)}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left[w(z,t)c(z,t)-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial c(z,t)}{\partial z}}\right],\end{array}$$

where w(z,t) denotes the ensemble-averaged particle velocity and κ_z(z,t) denotes the ensemble-expected particle diffusivity. For steady conditions involving a one-dimensional mean motion, $\partial c(z,t)/\partial t=\mathrm{0}$, $c(z,t)\to c\left(z\right)$, $w(z,t)\to W$, and Eq. (C1) becomes

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[Wc\left(z\right)-{\mathit{\kappa }}_z\left(z\right){\displaystyle \frac{\mathrm{d}c\left(z\right)}{\mathrm{d}z}}\right]=\mathrm{0}.\end{array}$$

Moreover, the particle flux must be zero across any surface normal to z, so

$$\begin{array}{}\text{(C3)}& Wc\left(z\right)-{\mathit{\kappa }}_z\left(z\right){\displaystyle \frac{\mathrm{d}c\left(z\right)}{\mathrm{d}z}}=\mathrm{0}.\end{array}$$

Under steady conditions the total number of particles with finite (measurable, non-saturated) OSL age within the soil element remains fixed. A particle entering the soil cannot attain a finite OSL age until it reaches the surface and is bleached, and then it becomes buried and exposed to the dose field. Thus, even though particles that eventually possess a finite OSL age continuously enter the element through its lower boundary, this boundary must be considered a zero flux boundary, as no particle with finite age can be added to the soil. Particles at the soil surface with zero OSL age are removed by erosion. The erosion rate matches W, so the rate of loss of particles is exactly balanced by the rate at which particles reach the surface and become OSL particles (with an OSL age of zero), many of which then take random walks downward. Thus, the upper boundary also must be considered a zero flux boundary with fixed concentration c(h).

With uniform mixing (κ_z=K_z) we integrate Eq. (C3) to obtain

$$\begin{array}{}\text{(C4)}& c\left(z\right)=C_{\mathrm{1}}{e}^{Wz/K_z},\end{array}$$

with constant of integration C₁. The boundary condition $c\left(h\right)=C_{\mathrm{1}}{e}^{Wh/K_z}$ then leads to the solution

$$\begin{array}{}\text{(C5)}& c\left(z\right)=c\left(h\right)e^{-W(h-z)/K_z}.\end{array}$$

The boundary condition c(h) should be equal to the concentration of OSL-sensitive particles entering the base of the soil element, but which only take on finite OSL ages once they reach the surface and are bleached. With nonuniform mixing (${\mathit{\kappa }}_z\left(z\right)=K_{z}z/h$) we integrate Eq. (C3) to obtain

$$\begin{array}{}\text{(C6)}& c\left(z\right)=c\left(h\right){\left({\displaystyle \frac{z}{h}}\right)}^{Wh/K_z}.\end{array}$$

With dimensionless height $\widehat{z}=z/h$, dimensionless concentration $\widehat{c}\left(\widehat{z}\right)=c\left(z\right)/c\left(h\right)$ and Péclet number $\mathit{Pe}=Wh/K_z$, Eqs. (C5) and (C6) become

$$\begin{array}{}\text{(C7)}& \widehat{c}\left(\widehat{z}\right)=e^{-\mathit{Pe}(\mathrm{1}-\widehat{z})}\end{array}$$

and

$$\begin{array}{}\text{(C8)}& \widehat{c}\left(\widehat{z}\right)={\widehat{z}}^{\mathit{Pe}},\end{array}$$

which are Eqs. (41) and (42) in the text. These results (Fig. C1) are used next in obtaining the expected OSL age A(z) of particles.

Figure C1 Plot of dimensionless OSL particle concentration $\widehat{c}\left(\widehat{z}\right)=c\left(z\right)/c\left(h\right)$ versus dimensionless height $\widehat{z}=z/h$ with uniform mixing (solid lines) and nonuniform mixing (dashed lines) as these vary with the Péclet number $\mathit{Pe}=Wh/K_z$.

Appendix D: Conservation of expected particle OSL age

Let m₁(z,t) denote the expected (average) finite OSL age of particles within the small interval z to z+dz in a soil element with dimensions XYh. With a total of N such particles within the element, the product $Nc(z,t)m_{\mathrm{1}}(z,t)XY\mathrm{d}z$ represents the total (collective) OSL age of particles within dz. The product $c(z,t)m_{\mathrm{1}}(z,t)$ satisfies a Fokker–Planck-like equation, namely,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left[c\right(z,t\left)m_{\mathrm{1}}\right(z,t\left)\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial }{\partial z}}\left(w\right(z,t\left)c\right(z,t\left)m_{\mathrm{1}}\right(z,t)\\ \text{(D1)}& {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z(z,t){\displaystyle \frac{\partial }{\partial z}}\left[c\right(z,t\left)m_{\mathrm{1}}\right(z,t\left)\right])+Sc(z,t),\end{array}$$

where w(z,t) denotes the ensemble-averaged particle velocity, κ_z(z,t) denotes the ensemble-expected particle diffusivity, and S is a source term. For steady conditions in both c and m₁ involving a one-dimensional mean motion, $(\partial /\partial t)\left[c\right(z,t\left)m_{\mathrm{1}}\right(z,t\left)\right]=\mathrm{0}$, $c(z,t)\to c\left(z\right)$, $m_{\mathrm{1}}(z,t)\to m_{\mathrm{1}}\left(z\right)$, $w(z,t)\to W$, and Eq. (D1) becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(Wc\left(z\right)m_{\mathrm{1}}\left(z\right)-{\mathit{\kappa }}_z\left(z\right){\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[c\right(z\left)m_{\mathrm{1}}\right(z\left)\right]\right)\\ \text{(D2)}& {\displaystyle }& {\displaystyle }=Sc\left(z\right).\end{array}$$

We now rewrite Eq. (D2) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(m_{\mathrm{1}}\right(z)\left[Wc\left(z\right)-{\mathit{\kappa }}_z{\displaystyle \frac{\mathrm{d}c\left(z\right)}{\mathrm{d}z}}\right]\\ \text{(D3)}& {\displaystyle }& {\displaystyle }-{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}})=Sc(z).\end{array}$$

Using Eq. (B3), this reduces to

$$\begin{array}{}\text{(D4)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[-{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}\right]=Sc\left(z\right).\end{array}$$

Indicating that the flux $q_z\left(z\right)=-{\mathit{\kappa }}_z\left(z\right)c\left(z\right)\mathrm{d}{m}_{\mathrm{1}}\left(z\right)/\mathrm{d}z$ of particle OSL age involves only diffusion. That is, OSL age is not advected. In this problem, particles are advected, carrying their OSL age with them, but particle advection is balanced by particle diffusion.

We now write Eq. (D4) as

$$\begin{array}{}\text{(D5)}& {\displaystyle \frac{\mathrm{d}{q}_z\left(z\right)}{\mathrm{d}z}}=Sc\left(z\right).\end{array}$$

Vertically integrating,

$$\begin{array}{}\text{(D6)}& \underset{\mathrm{0}}{\overset{h}{\int }}{\displaystyle \frac{\mathrm{d}{q}_z\left(z\right)}{\mathrm{d}z}}\mathrm{d}z=S\underset{\mathrm{0}}{\overset{h}{\int }}c\left(z\right)\mathrm{d}z.\end{array}$$

With fixed boundaries, this yields

$$\begin{array}{}\text{(D7)}& q_z\left(h\right)-q_z\left(\mathrm{0}\right)=S\stackrel{\mathrm{‾}}{c}h,\end{array}$$

where the overbar denotes a vertically averaged quantity. Moreover, note that q_z(0)=0, as particles with a finite OSL age cannot be imported to the soil element. Thus,

$$\begin{array}{}\text{(D8)}& q_z\left(h\right)=-{\mathit{\kappa }}_z\left(h\right)c\left(h\right){\left[{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}\right]}_{z=h}=S\stackrel{\mathrm{‾}}{c}h.\end{array}$$

This indicates that OSL age is diffused “through” the soil surface at a rate equal to its total production within the soil column.

D1 Uniform mixing

With uniform mixing (κ_z=K_z) we write

$$\begin{array}{}\text{(D9)}& q_z\left(h\right)=S\underset{\mathrm{0}}{\overset{h}{\int }}c\left(z\right)\mathrm{d}z.\end{array}$$

Using Eq. (C5), this becomes

$$\begin{array}{}\text{(D10)}& q_z\left(h\right)=Sc\left(h\right)e^{-Wh/K_z}\underset{\mathrm{0}}{\overset{h}{\int }}{e}^{Wz/K_z}\mathrm{d}z.\end{array}$$

Evaluating the integral then gives

$$\begin{array}{}\text{(D11)}& q_z\left(h\right)={\displaystyle \frac{Sc\left(h\right)K_z}{W}}\left(\mathrm{1}-e^{-Wh/K_z}\right).\end{array}$$

More generally,

$$\begin{array}{}\text{(D12)}& q_z\left(z\right)=Sc\left(h\right)e^{-Wh/K_z}\int e^{Wz/K_z}\mathrm{d}z.\end{array}$$

Evaluating the integral,

$$\begin{array}{}\text{(D13)}& q_z\left(z\right)={\displaystyle \frac{Sc\left(h\right)K_z}{W}}{e}^{-Wh/K_z}{e}^{Wz/K_z}+C_{\mathrm{1}}.\end{array}$$

Using the boundary condition obtained above for q_z(h),

$$\begin{array}{}\text{(D14)}& C_{\mathrm{1}}=-{\displaystyle \frac{Sc\left(h\right)K_z}{W}}{e}^{-Wh/K_z}.\end{array}$$

Using this result and $q_z\left(z\right)=-K_{z}c\left(z\right)\mathrm{d}{m}_{\mathrm{1}}\left(z\right)/\mathrm{d}z$ with c(z) given by Eq. (C5), we obtain

$$\begin{array}{}\text{(D15)}& {\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}=-{\displaystyle \frac{S}{W}}+{\displaystyle \frac{S}{W}}{e}^{-Wz/K_z}.\end{array}$$

Integrating and using the boundary condition that m₁(h)=0 then yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{m}_{\mathrm{1}}\left(z\right)={\displaystyle \frac{S}{W}}(h-z)\\ \text{(D16)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{S{K}_z}{W^{\mathrm{2}}}}\left(e^{-Wh/K_z}-e^{-Wz/K_z}\right).\end{array}$$

With dimensionless height $\widehat{z}=z/h$, dimensionless OSL age ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)=(W/h)m_{\mathrm{1}}\left(z\right)$ and Péclet number $\mathit{Pe}=Wh/K_z$, Eq. (D16) becomes

$$\begin{array}{}\text{(D17)}& {\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)=S(\mathrm{1}-\widehat{z})+{\displaystyle \frac{S{e}^{-\mathit{Pe}}}{\mathit{Pe}}}\left[\mathrm{1}-e^{\mathit{Pe}(\mathrm{1}-\widehat{z})}\right],\end{array}$$

which is Eq. (43) in the text.

D2 Nonuniform mixing

With nonuniform mixing (${\mathit{\kappa }}_z\left(z\right)=K_{z}z/h$) we use Eq. (C6) and write

$$\begin{array}{}\text{(D18)}& q_z\left(h\right)=Sc\left(h\right)h^{-\mathit{Pe}}\underset{\mathrm{0}}{\overset{h}{\int }}{z}^{\mathit{Pe}}\mathrm{d}z.\end{array}$$

Evaluating the integral then gives

$$\begin{array}{}\text{(D19)}& q_z\left(h\right)={\displaystyle \frac{Sc\left(h\right)h}{\mathrm{1}+\mathit{Pe}}}.\end{array}$$

More generally,

$$\begin{array}{}\text{(D20)}& q_z\left(z\right)=Sc\left(h\right)h^{-\mathit{Pe}}\int z^{\mathit{Pe}}\mathrm{d}z.\end{array}$$

Evaluating the integral,

$$\begin{array}{}\text{(D21)}& q_z\left(z\right)={\displaystyle \frac{Sc\left(h\right)h^{-\mathit{Pe}}}{\mathrm{1}+\mathit{Pe}}}{z}^{\mathrm{1}+\mathit{Pe}}+C_{\mathrm{1}}.\end{array}$$

Using Eq. (D19), C₁=0. Then using $q_z\left(z\right)=-K_z(z/h)c\left(z\right)\mathrm{d}{m}_{\mathrm{1}}\left(z\right)/\mathrm{d}z$ together with Eq. (D19) we obtain

$$\begin{array}{}\text{(D22)}& {\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}=-{\displaystyle \frac{Sh}{K_z(\mathrm{1}+\mathit{Pe})}}.\end{array}$$

Integrating and using the boundary condition that m₁(h)=0 then yields

$$\begin{array}{}\text{(D23)}& m_{\mathrm{1}}\left(z\right)={\displaystyle \frac{S{h}^{\mathrm{2}}}{K_z(\mathrm{1}+\mathit{Pe})}}\left(\mathrm{1}-{\displaystyle \frac{z}{h}}\right).\end{array}$$

With dimensionless height $\widehat{z}=z/h$, dimensionless OSL age ${\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)=(W/h)m_{\mathrm{1}}\left(z\right)$ and Péclet number $\mathit{Pe}=Wh/K_z$, Eq. (D23) becomes

$$\begin{array}{}\text{(D24)}& {\widehat{m}}_{\mathrm{1}}\left(\widehat{z}\right)={\displaystyle \frac{S\mathit{Pe}}{\mathrm{1}+\mathit{Pe}}}(\mathrm{1}-\widehat{z}),\end{array}$$

which is Eq. (44) in the text.

Appendix E: Variance of OSL ages

For a set of particles possessing finite OSL ages within any interval dz, their rate of “aging” is fixed, independent of age. This means that the average OSL age increases at this fixed rate, whereas the second and higher moments do not change. Thus, direct production of the variance m₂ of OSL ages is zero.

For steady conditions we start with Eq. (40) in the text, namely,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[-{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\displaystyle \frac{\mathrm{d}{m}_{\mathrm{2}}\left(z\right)}{\mathrm{d}z}}\right]\\ \text{(E1)}& {\displaystyle }& {\displaystyle }-\mathrm{2}{\mathit{\kappa }}_z\left(z\right)c\left(z\right){\left[{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{1}}\left(z\right)}{\mathrm{d}z}}\right]}^{\mathrm{2}}=\mathrm{0}.\end{array}$$

In the following we go directly to nondimensional forms of this.

E1 Uniform mixing

With κ_z=K_z and $\widehat{q}=[h/K_{z}c(h\left)\right(h/W{)}^{\mathrm{2}}]q$ we start with

$$\begin{array}{}\text{(E2)}& {\displaystyle \frac{\mathrm{d}\widehat{q}}{\mathrm{d}\widehat{z}}}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}\widehat{z}}}\left(-\widehat{c}{\displaystyle \frac{\mathrm{d}{\widehat{m}}_{\mathrm{2}}}{\mathrm{d}\widehat{z}}}\right)=\mathrm{2}\widehat{c}{\left({\displaystyle \frac{\mathrm{d}{\widehat{m}}_{\mathrm{1}}}{\mathrm{d}\widehat{z}}}\right)}^{\mathrm{2}}.\end{array}$$

Taking the derivative of Eq. (D17), squaring the result and using $\widehat{c}=e^{-\mathit{Pe}}{e}^{\mathit{Pe}\widehat{z}}$ then leads to

$$\begin{array}{}\text{(E3)}& {\displaystyle \frac{\mathrm{d}\widehat{q}}{\mathrm{d}\widehat{z}}}=\mathrm{2}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}{e}^{\mathit{Pe}\widehat{z}}-\mathrm{4}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}+\mathrm{2}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}{e}^{-\mathit{Pe}\widehat{z}}.\end{array}$$

Integrating this with respect to $\widehat{z}$ from $\widehat{z}=\mathrm{0}$ to $\widehat{z}=\mathrm{1}$ and noting that $\widehat{q}\left(\mathrm{0}\right)=\mathrm{0}$,

$$\begin{array}{}\text{(E4)}& \widehat{q}\left(\mathrm{1}\right)={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}}{\mathit{Pe}}}-\mathrm{4}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}-{\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{e}^{-\mathrm{2}\mathit{Pe}}}{\mathit{Pe}}}.\end{array}$$

More generally,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{q}\left(\widehat{z}\right)={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}}{\mathit{Pe}}}{e}^{\mathit{Pe}\widehat{z}}-\mathrm{4}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}\widehat{z}\\ \text{(E5)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{e}^{-\mathit{Pe}}}{\mathit{Pe}}}{e}^{-\mathit{Pe}\widehat{z}}+C_{\mathrm{1}}.\end{array}$$

Using Eq. (E5) gives C₁=0. With $\widehat{q}=-\widehat{c}\mathrm{d}{\widehat{m}}_{\mathrm{2}}/\mathrm{d}\widehat{z}$, and again using $\widehat{c}=e^{-\mathit{Pe}}{e}^{\mathit{Pe}\widehat{z}}$,

$$\begin{array}{}\text{(E6)}& {\displaystyle \frac{\mathrm{d}{\widehat{m}}_{\mathrm{2}}}{\mathrm{d}\widehat{z}}}=\mathrm{4}{S}^{\mathrm{2}}\widehat{z}{e}^{-\mathit{Pe}\widehat{z}}-{\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}}{\mathit{Pe}}}+{\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}}{\mathit{Pe}}}{e}^{-\mathrm{2}\mathit{Pe}\widehat{z}}.\end{array}$$

Integrating with respect to $\widehat{z}$ and evaluating the constant of integration with the condition that ${\widehat{m}}_{\mathrm{2}}\left(\mathrm{1}\right)=\mathrm{0}$ then yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}}{\mathit{Pe}}}(\mathrm{1}-\widehat{z})+{\displaystyle \frac{\mathrm{4}{S}^{\mathrm{2}}}{{\mathit{Pe}}^{\mathrm{2}}}}\left[\right(\mathrm{1}+\mathit{Pe})e^{-\mathit{Pe}}\\ \text{(E7)}& {\displaystyle }& {\displaystyle }-(\mathrm{1}+\mathit{Pe}\widehat{z})e^{-\mathit{Pe}\widehat{z}}]+{\displaystyle \frac{S^{\mathrm{2}}}{{\mathit{Pe}}^{\mathrm{2}}}}\left(e^{-\mathrm{2}\mathit{Pe}}-e^{-\mathrm{2}\mathit{Pe}\widehat{z}}\right),\end{array}$$

which is Eq. (45) in the text.

E2 Nonuniform mixing

With $\widehat{q}=[h/K_{z}c(h\left)\right(h/W{)}^{\mathrm{2}}]q$ we start with

$$\begin{array}{}\text{(E8)}& {\displaystyle \frac{\mathrm{d}\widehat{q}}{\mathrm{d}\widehat{z}}}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}\widehat{z}}}\left(-\widehat{z}\widehat{c}{\displaystyle \frac{\mathrm{d}{\widehat{m}}_{\mathrm{2}}}{\mathrm{d}\widehat{z}}}\right)=\mathrm{2}\widehat{z}\widehat{c}{\left({\displaystyle \frac{\mathrm{d}{\widehat{m}}_{\mathrm{1}}}{\mathrm{d}\widehat{z}}}\right)}^{\mathrm{2}},\end{array}$$

with $\widehat{c}={\widehat{z}}^{\mathit{Pe}}$ then $\widehat{z}\widehat{c}={\widehat{z}}^{\mathrm{1}+\mathit{Pe}}$. Using this and taking the derivative of Eq. (D24) with respect to $\widehat{z}$ and squaring the result leads to

$$\begin{array}{}\text{(E9)}& {\displaystyle \frac{\mathrm{d}\widehat{q}}{\mathrm{d}\widehat{z}}}={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{\mathit{Pe}}^{\mathrm{2}}}{(\mathrm{1}+\mathit{Pe}{)}^{\mathrm{2}}}}{\widehat{z}}^{\mathrm{1}+\mathit{Pe}}.\end{array}$$

Integrating this with respect to $\widehat{z}$ from $\widehat{z}=\mathrm{0}$ to $\widehat{z}=\mathrm{1}$ and noting that $\widehat{q}\left(\mathrm{0}\right)=\mathrm{0}$,

$$\begin{array}{}\text{(E10)}& \widehat{q}\left(\mathrm{1}\right)={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{\mathit{Pe}}^{\mathrm{2}}}{(\mathrm{2}+\mathit{Pe})(\mathrm{1}+\mathit{Pe}{)}^{\mathrm{2}}}}.\end{array}$$

More generally,

$$\begin{array}{}\text{(E11)}& \widehat{q}\left(z\right)={\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{\mathit{Pe}}^{\mathrm{2}}}{(\mathrm{2}+\mathit{Pe})(\mathrm{1}+\mathit{Pe}{)}^{\mathrm{2}}}}{\widehat{z}}^{\mathrm{2}+\mathit{Pe}}+C_{\mathrm{1}}.\end{array}$$

Using Eq. (E10) gives C₁=0. With $\widehat{q}=-\widehat{z}\widehat{c}\mathrm{d}{\widehat{m}}_{\mathrm{2}}/\mathrm{d}\widehat{z}$, and again using $\widehat{z}\widehat{c}={\widehat{z}}^{\mathrm{1}+\mathit{Pe}}$,

$$\begin{array}{}\text{(E12)}& {\displaystyle \frac{\mathrm{d}{\widehat{m}}_{\mathrm{2}}}{\mathrm{d}\widehat{z}}}=-{\displaystyle \frac{\mathrm{2}{S}^{\mathrm{2}}{\mathit{Pe}}^{\mathrm{2}}}{(\mathrm{2}+\mathit{Pe})(\mathrm{1}+\mathit{Pe}{)}^{\mathrm{2}}}}\widehat{z}.\end{array}$$

Integrating with respect to $\widehat{z}$ and evaluating the constant of integration with the condition that ${\widehat{m}}_{\mathrm{2}}\left(\mathrm{1}\right)=\mathrm{0}$ then yields

$$\begin{array}{}\text{(E13)}& {\widehat{m}}_{\mathrm{2}}\left(\widehat{z}\right)={\displaystyle \frac{S^{\mathrm{2}}{\mathit{Pe}}^{\mathrm{2}}}{(\mathrm{2}+\mathit{Pe})(\mathrm{1}+\mathit{Pe}{)}^{\mathrm{2}}}}(\mathrm{1}-{\widehat{z}}^{\mathrm{2}}),\end{array}$$

which is Eq. (46) in the text.