MATLAB code for computing steady-state mineral concentrations and soil thickness in regolith-mantled hillslopes, with application to soil organic carbon storage. Supplemental material for Roering et al. (2023).
Roering, J. J., Hunter, B. D., Ferrier, K. L., Chadwick, O. A., Yoo, K., Wackett, A. A., Almond, P. C., Silva, L., and Jelline, A. M. (2023). Quantifying erosion rates and weathering pathways that maximize soil organic carbon storage. Biogeochemistry. https://doi.org/10.1007/s10533-023-01054-7
These files are used to solve the equations for steady-state mineral concentration and soil thickness (equations 7 and 8). The files are modified versions of those developed by Ferrier and West (2017). The modifications include accounting for the alteration of mineral phases and the production of secondary mineral phases. The numerical implementation otherwise follows that presented by Ferrier and West (2017).
lag_parameters.m is used to specify the model parameters described in the manuscript and follows the conventions of Ferrier and West (2017). steadyq_nd_eqns.m sets up the soil thickness and mineral concentration equations. Notably, the mineral concentration equations are mineral-specific: the allophane, plagioclase, and kaolinite equations account for gains and/or losses as described in the text. steadyq_nd.m invokes the MATLAB equation solver and outputs soil thickness and mineral concentrations, as well as additional weathering metrics used by Ferrier and West (2017) that are not invoked here. steadyEplot.m sets up a recursive loop to solve for mineral concentrations and soil thickness for an array of steady-state erosion rate values specified by the user, and generates plots of erosion rate vs. allophane concentration and vs. all mineral concentrations. pcm_2D.m calculates the erosion rate associated with peak allophane concentration and the peak allophane concentration itself for a range of rate parameter values.
- MATLAB (developed and tested with standard MATLAB; no additional toolboxes required beyond the Optimization Toolbox for
fsolve)
To reproduce Figure 2 from the paper:
E = [0.01:0.01:0.99]; % vector of non-dimensional erosion rate values
p = lag_parameters; % load parameter values
[Ec, conc, h] = steadyEplot(p, E); % solve and plotTo compute peak allophane concentration over a parameter grid:
p = lag_parameters;
E = [0.01:0.01:0.99];
kplag = logspace(-3, 0, 20);
kallo = logspace(-3, 0, 20);
[E_nd_max, allo_max] = pcm_2D(E, p, kplag, kallo);Ferrier, K. L., and West, N. (2017). Responses of chemical erosion rates to transient perturbations in physical erosion rates, and implications for relationships between chemical and physical erosion rates in regolith-mantled hillslopes. Earth and Planetary Science Letters, 474, 447–456. https://doi.org/10.1016/j.epsl.2017.07.002
MIT — see LICENSE for details.