results = get_al_be_erosion(sample,consts,nuclide)
This is the main control function that carries out the erosion rate calculation.
al_be_erosion_many call it.
sample is a structure containing sample information. The fields are as follows:
|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 cm||double|
|sample.N10||Be concentration, atoms g yr||double|
|sample.delN10||standard error of Be concentration||double|
|sample.N26||Al concentration, atoms g yr||double|
|sample.delN26||standard error of Al concentration||double|
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.
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:
|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 g yr)||double|
|results.P_St||Surface production rate due to spallation (atoms g yr) following non-time-dependent Lal (1991) - Stone(2000) scaling scheme||double|
|results.Egcm2yr||1x5 vector containing erosion rates (g cm yr) according to the St, De, Du, Li, and Lm scaling schemes respectively||1 x 5 double|
|results.EmMyr||1x5 vector containing erosion rates (m Myr) 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 Myr) 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 Myr) according to the 5 scaling schemes.||1 x 5 double|
|results.fzero_status||1x5 vector containing output flags returned by rootfinding function
||1 x 5 double|
|results.fval||1x5 vector containing objective function values (atoms g) 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 .||double|
1. Basic idea.
This function solves the equation:
for the erosion rate (here in g cm
is the effective attenuation length for spallogenic production (g cm), is the decay constant for the nuclide in question (yr), is the measured nuclide concentration (atoms g),
is the thickness-averaged surface nuclide production rate due to spallation as a function of time (atoms g
is the thickness-averaged production rate due to muons (atoms g
yr) as a function of the depth (g cm) . 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
2. Initial guess.
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:
Here we calculate using the non-time-dependent St scaling scheme.
This can be solved directly for to obtain the initial guess.
3. Calculate the production rates; pass to objective function; solve.
We then assemble the functions (for each scaling scheme) and 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
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.
is generated on a 100-point log-spaced vector from 0 to
g cm 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
4. Error propagation.
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 is the decay constant for the nuclide of interest, is the effective surface production rate due to spallation, is the surface production rate due to muons, is the effective attenuation length for production by spallation, and is the effective attenuation length for production by muons. Both and 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 , the nuclide concentration attributable to production by muons, in the course of solving the full equation above, we can calculate the value of for which the simplified equation would yield the true erosion rate by:
where is the `true' erosion rate calculated by solving the full equation.
In like manner, we have already calculated , the nuclide concentration attributable to production by spallation, so we can also calculate the value of 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
, 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.