Uncertainties in modelling diagenetic self-organisation in limestone-marl sequences

Rhythmic variations in the properties of sediments are commonly used as archives of paleoclimate changes driven by variation in insolation caused by the changes in the Earth's orbit and the tilt of its axis. But rhythmicity can also arise from diagenetic self-organization. Distinguishing between these two drivers requires simulating self-organization. We started with an attempt to reproduce the main results from a paper by Ivan L'Heureux (2018)¹ - who proposed a mathematical model of a nonlinear dynamical system, in which self-organized oscillations arise from homogenous initial sediment and result in sediment layers with different compositions. The model consists of five stiff differential equations, for the composition of calcite and aragonite, two mineral polymorphs of CaCO3, of which aragonite is metastable, for the concentrations of calcium and carbonate ions in the pore water and for the porosity, as functions of depth and time. The self-organized patterns are in this model the result of two processes happening at different temporal scales: rapid dissolution of aragonite and slow sediment compression in response to increased porosity as aragonite is removed from the solid phase. Reproducing the steady-state distributions along depth required a major effort, mostly with regard to understanding what triggers numerical instabilities, but was finally successful.  
Currently, we have not yet succeeded in reproducing oscillations, that L'Heureux predicted, without requiring an external force, for high initial and boundary sediment porosity.  It is essential that we are able to determine for which inital and boundary conditions oscillations should occur, beyond the uncertainties introduced by numerical algorithms for solving partial differential equations, e.g. for many sets of parameters the integrations over time can easily "derail". We have formulated two questions that we want to share with the audience in order to seek help. These are our questions: 

1) Do the five differential equations describe the underlying physics adequately?  
2) Our current software implementation of the five differential equations does not yield any oscillations, is that a flaw on our side, or does this agree with mathematical insights?

• "Diagenetic Self-Organization and Stochastic Resonance in a Model of Limestone-Marl Sequences" by Ivan L'Heureux (2018). https://doi.org/10.1155/2018/4968315