1 Scientific rationale

From the global ocean to engineered ponds and reservoirs, sediments play a major role in regulating the ecology, biological productivity and environmental quality of surface water bodies (Hedges, 1992;Jones and Holmes, 1996;Bratvold and Browdy, 2001;Juracek, 2012;Breithaupt et al., 2012). Chemical transformations and (bio)physical processes occurring near the sediment-water interface (SWI) profoundly affect the recycling and speciation of carbon, nutrients, metals and contaminants (Berner, 1980;Canavan et al., 2006;Burdige, 2006;Couture et al., 2010). The associated chemical exchanges across the SWI directly impact the biogeochemical functioning of the sediments and overlying water. The release of dissolved phosphate (PO4) from bottom sediments, for example, can trigger algal blooms in lakes and reservoirs (Sondergaard et al., 2005), while the mobility and toxicity of pollutants may vary dramatically along the steep geochemical gradients found in aquatic sediments (Couture et al., 2010;Bessinger et al., 2012). Reactive transport models (RTMs) are useful tools to develop a process-based understanding of sediment biogeochemistry (Berner, 1980;Van Cappellen and Gaillard, 1996;Boudreau, 1997, 1999;Meysman et al., 2003;Burdige, 2006). RTMs of aquatic sediments have reached high levels of sophistication in the representation of the complex reaction networks that control element cycling and redox dynamics in sediments (Wang and Van Cappellen, 1996;Luff et al., 2001;Meysman et al., 2003;Katsev et al., 2004;Katsev and Dittrich, 2013;Regnier et al., 2011;Schulz, 2006 ). Through the continuous integration of new observational and theoretical knowledge, sediment RTMs have remained abreast of advances in biogeochemistry and transport theory. In addition to the quantitative interpretation of pore water and solid sediment data, these RTMs can help identify gaps in our conceptual understanding of the functioning of aquatic sediments, assess uncertainties in model structure and parameter values, and predict the response of benthic processes to changes in external forcing.

Despite their successes and proven potential, sediment RTMs are used by only a small group of researchers. Students in fields where sediment biogeochemistry is typically part of the curriculum, for example, geochemistry, aquatic ecology, limnology, oceanography or environmental engineering, rarely receive in-depth training in the application of sediment RTMs. Not surprising then, these models tend to be underutilized by environmental professionals. Even detailed data sets combining, say, pore water profiles, solid sediment geochemistry, benthic fluxes, and sedimentation rate determinations, often do not benefit from a quantitative analysis based on RTM calculations.

One of the main obstacles to a more systematic use of sediment RTMs is that many existing models are problem- and site-specific (e.g., Meysman et al. (2003) and, hence, not necessarily easily transferable to a new user’s application. In addition, most codes have been generated in traditional programming languages, such as FORTRAN or C++. Therefore, new applications require the rewriting and compilation of the source codes, which tends to discourage many potential users with limited programming skills. The widespread availability of high-level mathematical software environments, however, should help to overcome these obstacles. Interactive programming packages, such as MAPLE®, Mathematica®, MATLAB® and open source R, provide flexible platforms in which an individual user can easily develop, adapt, test and manage models. By avoiding the hurdles of traditional programming languages, the users can focus on the definition, performance and application of their model, rather than on the technicalities of writing code. Because the new programming environments provide many visualization capabilities, the graphical analysis of the results is also greatly facilitated and enhanced.

An example of the efforts to make sediment RTMs more flexible and user-friendly is the BRNS (Biogeochemical Reaction Network Simulator) developed by Regnier et al. (2003). It is based on the symbolic programming language MAPLE® and offers a modelling framework within which a variable number of reaction processes can be quantitatively described and interfaced with one-dimensional transport processes by means of an expandable, web-distributed Knowledge Base (KB) of process equations and parameters. The user builds their own model application in a MAPLE® template. The template is then sent to an automatic code generator (ACG), which in turn generates an executable FORTRAN code on the home server. The server returns the executable file to the users who can run it on their personal computer. That is, the user no longer has to change and compile the source code herself. While the BRNS-KB-ACG offers a general, adaptive platform for early diagenetic applications, it requires sustained institutional support to periodically upgrade the KB, maintain the server and provide trouble-shooting.

Here, we present MATSEDLAB, a MATLAB® script for one-dimensional sediment reactive transport calculations. We describe the theoretical background of the model and introduce our approach to solve early diagenetic problems in MATLAB®, followed by four applications. For each application, we introduce the key processes and explain how they are implemented in the MATSEDLAB script. The applications illustrate how to modify boundary conditions, add new chemical species and add new biogeochemical processes to the model. All the graphical results presented are generated using built-in MATLAB® functions. Model outputs include vertical pore water plus solid sediment concentration profiles, benthic fluxes and integrated rates of biogeochemical processes.