get_al_be_erosion.m

This is the main control function that carries out the erosion rate calculation. al_be_erosion_one and
al_be_erosion_many call it.

The argument sample is a structure containing sample information. The fields are as follows:

sample.sample_name	Sample name	string
sample.lat	Latitude	double
sample.long	Longitude	double
sample.elv	elevation in meters	double
sample.pressure	pressure in hPa (optional if sample.elv is set)	double
sample.aa	Flag that indicates how to interpret the elevation value.	string
sample.thick	Sample thickness in cm	double
sample.rho	Sample density, g $\cdot$ cm $^{-3}$	double
sample.othercorr	Shielding correction.	double
sample.N10	$^{10}$ Be concentration, atoms $\cdot$ g $^{-1} \cdot$ yr $^{-1}$	double
sample.delN10	standard error of $^{10}$ Be concentration	double
sample.N26	$^{26}$ Al concentration, atoms $\cdot$ g $^{-1} \cdot$ yr $^{-1}$	double
sample.delN26	standard error of $^{26}$ Al concentration	double

The argument consts is a structure containing the constants. It is typically the structure created by
make_al_be_consts_v12.m, although only a subset of the fields in that structure are actually used by this function.

The argument nuclide tells the function which nuclide is being used; allowed values are 10 or 26. This is a numerical value, not a string.

This function returns a structure called results that contains the following fields:

Non-scaling-scheme-specific information:
results.flags	Non-fatal error messages, mostly to do with nuclide concentrations above saturation	string
results.main_version	Version number for this function	string
results.obj_version	Version of objective function called internally	string
results.mu_version	Version of P_mu_total called internally	string
results.Pmu0	Surface production rate due to muons (atoms $\cdot$ g $^{-1} \cdot$ yr $^{-1}$ )	double
Scaling-scheme-specific information:
results.P_St	Surface production rate due to spallation (atoms $\cdot$ g $^{-1} \cdot$ yr $^{-1}$ ) following non-time-dependent Lal (1991) - Stone(2000) scaling scheme	double
results.Egcm2yr	1x5 vector containing erosion rates (g $\cdot$ cm $^{2}\cdot$ yr $^{-1}$ ) according to the St, De, Du, Li, and Lm scaling schemes respectively	1 x 5 double
results.EmMyr	1x5 vector containing erosion rates (m $\cdot$ Myr $^{-1}$ ) according to the St, De, Du, Li, and Lm scaling schemes respectively	1 x 5 double
results.delE_int	1x5 vector containing internal uncertainties (m $\cdot$ Myr $^{-1}$ ) according to the 5 scaling schemes. Note that these should all be essentially equal.	1 x 5 double
results.delE_ext	1x5 vector containing external uncertainties (m $\cdot$ Myr $^{-1}$ ) according to the 5 scaling schemes.	1 x 5 double
Diagnostic information:
results.fzero_status	1x5 vector containing output flags returned by rootfinding function `fzero` at solving for the erosion rates from the 5 scaling schemes.	1 x 5 double
results.fval	1x5 vector containing objective function values (atoms $\cdot$ g $^{-1}$ ) at the erosion rates from the 5 scaling schemes.	1 x 5 double
results.time_fzero	1x5 vector containing the time (s) required for each solution.	1 x 5 double
results.time_mu_precalc	Time (s) required to calculate $P_{\mu}(z)$ .	double

for the erosion rate $\epsilon$ (here in g $\cdot$ cm $^{-2} \cdot$ yr $^{-1}$ ), where $\Lambda_{sp}$ is the effective attenuation length for spallogenic production (g $\cdot$ cm $^{-2}$ ), $\lambda$ is the decay constant for the nuclide in question (yr $^{-1}$ ), $N_{m}$ is the measured nuclide concentration (atoms $\cdot$ g $^{-1}$ ), $P_{sp,0}(t)$ is the thickness-averaged surface nuclide production rate due to spallation as a function of time (atoms $\cdot$ g $^{-1}\cdot$ yr $^{-1}$ ), and $P_{\mu}(z)$ is the thickness-averaged production rate due to muons (atoms $\cdot$ g $^{-1}\cdot$ yr $^{-1}$ ) as a function of the depth

(g $\cdot$ cm $^{-2}$ ) . As this cannot be solved analytically, we use the MATLAB rootfinding algorithm fzero to find the zero of an objective function, ET_objective.m, that computes the left-hand side of this equation, that is, the predicted nuclide concentration for a particular erosion rate less the measured nuclide concentration. The actual integration is described in the documentation for ET_objective.

The objective function involves several numerical integrations, so finding its zero can be slow. The easiest way to speed it up is to provide it with an initial guess for the solution that is close to the actual solution. We do this by first estimating the erosion rate using the simple equation of Lal (1991), which disregards production by muons:

We then assemble the functions $P_{sp,0}(t)$ (for each scaling scheme) and $P_{\mu}(z)$ and pass them to the objective function, then find its zero to obtain the corresponding erosion rate for each scaling scheme.

The actual assembly of $P_{sp,0}(t)$ for the four time-dependent scaling schemes is the same as is used in get_al_be_age.m, and is described above in section 1.3.

$P_{\mu}(z)$ is generated on a 100-point log-spaced vector from 0 to $2 \times 10^{5}$ g $\cdot$ cm $^{-2}$ using the function P_mu_total.m. The spacing of this vector affects the accuracy of the solution; this is discussed in more detail in the documentation for ET_objective.m.

The fact that the full erosion rate equation cannot be solved analytically also makes error propagation difficult. We propagate errors linearly via the standard error propagation formula, but we cannot compute the derivatives of the derived erosion rate with respect to the uncertain input parameters, that are required in the formula, analytically, so two additional solutions of the full equation are required to estimate the partial derivative with respect to each uncertain parameter. It is time-consuming to do this using the full forward model, so we use a simplified erosion rate - nuclide concentration relationship for the uncertainty analysis:

where $\lambda$ is the decay constant for the nuclide of interest, $P_{eff,sp}$ is the effective surface production rate due to spallation, $P_{\mu}$ is the surface production rate due to muons, $\Lambda_{sp}$ is the effective attenuation length for production by spallation, and $\Lambda_{\mu,eff}$ is the effective attenuation length for production by muons. Both $\Lambda_{\mu,eff}$ and $P_{eff,sp}$ depend on the erosion rate, which is why this simplified equation cannot be used to accurately calculate the erosion rate in the first place. However, as we have already calculated $N_{\mu}$ , the nuclide concentration attributable to production by muons, in the course of solving the full equation above, we can calculate the value of $\Lambda_{\mu,eff}$ for which the simplified equation would yield the true erosion rate by:

where $\epsilon_{true}$ is the `true' erosion rate calculated by solving the full equation.

In like manner, we have already calculated $N_{sp}$ , the nuclide concentration attributable to production by spallation, so we can also calculate the value of $P_{eff,sp}$ that would yield the correct erosion rate by:

Finally, the simplified erosion rate equation cannot be solved directly either, so it is coded as an internal function that can be used as an argument to fzero. We can then calculate the partial derivatives of the erosion rate with respect to the uncertain parameters $\partial \epsilon / \partial N$ , $\partial \epsilon / \partial P_{ref,Xx}$ , etc. numerically by perturbing the input parameters and repeatedly solving the simple erosion rate equation. We use a first-order centered difference scheme.

Finally, we use the resulting partial derivatives and input parameter uncertainties to evaluate the standard error propagation formula for internal and external uncertainties. Note that the internal uncertainty, that is, the uncertainty due to measurement error alone, is for all practical purposes the same for all the different scaling schemes. On the other hand, the external uncertainties, which reflect uncertainties in production rates as well, are different for each scaling scheme.