3 Note on implementation of diagenesis equation in Matlab

A simulation is initiated by executing the MATSEDLAB.m script. For easier reference, we divided the script into distinct blocks (Fig. 1). Block 1 in MATSEDLAB.m configures the spatial and temporal domains through the definition of the vectors x and t. By using x = linspace (a,b,xres), the user generates xres equally spaced grid points within the spatial domain [a, b]; similarly, the time span [t0, tf] is divided into tres intervals by defining the vector t = linspace (t0,tf,tres).

In order to solve the early diagenetic equation (Eq. 1), we chose the pdepe function in MATLAB. The pdepe function is designed to solve initial-boundary value problems (IBVPs) consisting of systems of parabolic and elliptic PDEs in one space variable and time. The numerical method is based on a piecewise nonlinear Petrov–Garlekin method with second-order accuracy. The method solves the ordinary differential equations (ODEs) resulting from the spatial discretization of the PDEs, using a built-in MATLAB® ODE solver to obtain approximate solutions at specified times within a defined time interval. An algorithm demonstrating the applicability of the method to a variety of advection–diffusion and diffusion problems has been derived (Skeel and Berzins, 1990). An important attribute of the solver is that the time step is determined dynamically by MATLAB® to insure the stability of the integration. However, the user can impose a minimum time step while spatial discretization is fully user-defined. Comparison of the accuracy of classical numerical schemes, including finite difference, shooting and collocation methods, with those of the pdepe function when solving 1D transport-reaction PDEs under steady-state conditions has been done (Yudianto and Xie, 2010). The results show that the pdepe function is more accurate, especially when using larger grid spacing, and yields shorter computation times when small grid spacing are imposed.

The pdepe solver and its input functions are automatically called via the command sol=pdepe(m,pdefun,icfun,bcfun,xmesh,tspan), where integer m imposes the symmetry of the problem (m can take values of 0, 1, and 2 corresponding to slab, cylindrical, or spherical symmetry, respectively; here we set it to 0 for 1-D early diagenesis problems), and xmesh (a<x<b) and tspan (t0<t<tf) are the depth and time steps that have been defined by the x and t vectors, respectively. pdefun, icfun and bcfun are functions defining the components of the PDEs, the initial conditions and the boundary conditions, respectively, and are detailed below.

pdefun is solved in Block 2 within temporal (t0<t<tf) and spatial (a<x<b) domains, and implemented under the following general form:

c(x,t,u,∂u/∂x) ∂u/∂t=x^(-m) ∂/∂x (x^m f(x,t,u,∂u/∂x))+s(x,t,u,∂u/∂x) (6)

in which vector u contains all the unknown variables (here the solid and solute concentrations). The coefficients of the time derivatives are collected in the diagonal matrix c. The coupling of the partial derivatives with respect to time is restricted to the multiplication by c. The diagonal elements of the c matrix can be either identically zero or positive. On the right hand side of Eq. (6), the functions f and s – the flux and sink/source terms, respectively – are vector functions, which depend on x, t , u and ∂u/∂x. The term f represents the total flux of constituent i; it encompasses advective fluxes (advection rate ϑ), diffusive fluxes by molecular diffusion (diffusion coefficient Dm) and sediment mixing by benthic organisms (bioturbation coefficient Db). The term s accounts for the net production or consumption of the given chemical species by (bio)geochemical kinetic reactions and, for solutes, also bioirrigation.

icfun (Block 3) defines the initial condition at t=t0 and all depths x and takes the form of: u(x,t_0 )=u_0 (x) (7) The function’s u = icfun(x) returns initial values of all the chemical species at depth x in the column vector u. The user may thus specify any set of initial depth profiles, including spatially heterogeneous ones. bcfun (Block 4) defines boundary conditions at x=a and x=b using the equation: p(x,t,u)+q(x,t)f(x,t,u,∂u⁄∂x)=0 (8) where f is the flux-vector from Eq. (6). bcfun has the form [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t), where ul is the numerical solution at the upper boundary, xl = a, and ur at the lower boundary, xr = b. In order to formulate the boundary conditions (Eqs. 3 and 4) using the pdepe formulation (Eq. 6) the column vectors pl and ql (at xl), and pr and qr (at xr) need to be defined. The commonly used formulations for boundary conditions, such as Dirichlet, Neumann and Cauchy/Robin are imbedded in Eq. (8) and for each type, the coefficients p(x,t,u) and q(x,t) take different values. As can be seen from Eq. (8), the boundary conditions can depend on t, which means that MATSEDLAB evaluates the boundary values at each time step. This gives the user the ability to impose transient boundary conditions. This is particularly useful when simulating the fate of compounds whose inputs are changing due to, for example, anthropogenic activity, or when dealing with systems where the bottom water chemistry varies over time. The MATSEDLAB script therefore allows non-steady state boundary conditions for the concentration of dissolved species and for the fluxes of solid species at the SWI.

The transformation of the early diagenetic Eqs. 1 to 5 into the format of Eqs. 6 to 8, through which the values of the arguments of the pdepe functionn are replaced by physical and biogeochemical terms, yields the matrices shown in Table 2 and detailed in their MATLAB® form in the user manual. The matrices c, f and s described above account for the total number of species L, that is, S pore water solutes and (L-S) solid-bound species.