Introduction of a new platform for parameter estimation of kinetically complex environmental systems.

A modeling framework (ReKinSim - Reaction Kinetics Simulator) is introduced, within which biogeochemical reactions in environmental systems can be described and inversely fitted to experimental data. Three key features of this simulation environment are: (1) a generic mathematical tool for solving sets of unlimited, arbitrary, non-linear ordinary differential equations; (2) no limitation to the number or type of reactions or other influential dynamics (e.g., isotope fractionation or small-scale mass-transfer limitations); (3) an easy to use and flexible module for nonlinear data-fitting. It allows users to easily define any kinetic model by a set of biogeochemical reactions relevant to the experimental application and to obtain the values of the kinetic parameters by fitting of the model to data. By allowing users to include the environmentally related processes and solving them along with the chemical kinetics, ReKinSim helps the user to elucidate the extent that these processes are controlled by factors other than kinetics. The novelty of the presented program primary lays in its unique combination of flexibility, computational efficiency and user-friendliness. ReKinSim's usability is showcased by four case studies of varying complexity, and compared against a set of currently available modeling tools.

Introduction

Over the past few decades significant improvements have been made in understanding the processes controlling the biogeochemical cycling of elements in environmental systems. Recent developments in the field of computational technology and algorithms have further helped scientists to simulate complicated biochemical processes and to enhance the exchange of information between recent scientific questions and models (e.g., Cirpka and Valocchi, 2007, Eckert et al., 2013, Hunkeler et al., 2009, König et al., 2017, Rosenzweig et al., 2014, Thullner et al., 2008). For a more quantitative analysis of such processes in natural or man-made environmental systems, it is important to acquire a fundamental understanding of the kinetics of degradation mechanisms and their dependency on environmental factors (Or et al., 2007, Thullner et al., 2005, Thullner et al., 2007, Van Loosdrecht et al., 1990). Such mechanisms and their interplay with other biogeochemical processes are described by a sequence of reactions/transformations which are then translated into mathematical formulations. Novel conclusions can be drawn when hypothetical descriptions of such mechanisms (usually in the form of models) are reconciled with experimental evidence. To this end, computational tools play a key role in stating the hypothetical functional relationships in the form of mathematical models and minimizing the error between suggested formulations and the observed data. For such a purpose, a tool for kinetic modeling and data fitting must (i) carry sophisticated algorithms equipped with the most recent updates and improvements, (ii) not be computationally elaborate, and (iii) be sufficiently versatile for implementation.

The current advances in kinetics modeling have often led to stand-alone platforms which are not modular, with outcomes that cannot be easily exported and linked to other codes (e.g., AQUASIM, Reichert, 1994). Furthermore, some of the programs suffer from the lack of an easily-operated interface and cannot be handled conveniently by scholars who are not strong modelers/informaticians. In special cases, the integration of numerous add-ons and often superfluous options has made some of the programs difficult to operate. Versatility is another issue in many of programs (e.g., ASM3, Iacopozzi et al. (2007), Global KinTek Explorer, Johnson et al. (2009), and Rayleigh Tool, Höhener and Atteia (2014) and Morrill et al. (2006)) that are confined to the solution of a few specific reaction mechanisms and thus are difficult to employ for a wider range of environmentally associated problems with more complicated dynamics. Addressing the above limitations while keeping the advantages of previously presented programs, we felt the need for a new program which better tackles environmentally related research questions. To that end, a newly developed computational tool (ReKinSim - Reaction Kinetics Simulator) is brought forward in this paper that is accurate, adjustable, user-friendly, and reasonably quick. The presented tool (ReKinSim) has been designed with the purpose of offering an easy-to-operate while flexible platform, free from extra mathematical complexities that are often confusing and sometimes misleading. Flexibility is another key attribute of ReKinSim, so there are no limitations in defining the kinetic mechanisms of interest or the actual number of chemicals in a given system. ReKinSim is coded in MATLAB language and the use of MATLAB numerical modules contributes greatly to its performance, extensibility, and competence. The potential relevance of the presented tool for biochemical/environmental systems is demonstrated and verified through several examples. Some case studies have previously been presented and served as challenging test cases using other optimization tools, including AQUASIM and Global KinTek Explorer. We revisited these case studies to compare the results from ReKinSim with those reported. The rest of the case studies were designed to test ReKinSim and are, in comparison, more complicated and challenging.

The overall aim of this contribution is to introduce a comprehensive tool for modeling of complex reaction mechanisms, enabling environmental engineers and biochemists to understand and quantify reaction kinetics/dynamics associated with their systems of interest. In this respect, ReKinSim can in principle be employed as an exploratory tool (e.g. using trial and error methods) to unravel the kinetics of the processes that are suspected to take place in a system. However, its main usage is to estimate the values of kinetic parameters by fitting a reaction mechanism to a set of experimental data. ReKinSim is developed with simplicity of usage in mind and unlike some other similarly oriented approaches (e.g., Ferrai et al., 2010, Fang et al., 2009, Kneis et al., 2017, Iacopozzi et al., 2007, Morrill et al., 2006, Regnier et al., 2002) it provides a simpler framework with a more approachable and less complicated set of options that lets scientists focus more on experimental and theoretical features of their research, and less on mathematical and numerical aspects. It is the goal of the program ReKinSim to allow users to specify models as freely as possible, to keep the complexities to a minimum, and to estimate with acceptable precision the values of kinetic parameters.

Fundamentals of using ReKinSim

The definition of reaction kinetics in ReKinSim is straightforward. The conceptual kinetic model, i.e. the definition of the influential processes (e.g., biogeochemical reactions) and the kinetic description of them in form of ordinary differential equations (ODEs), leads to a set of equations that need to be supplied to ReKinSim by user. For predictive reaction modeling only the relevant set of ODEs is required, which then can be solved by ReKinSim provided the values of all kinetic/dynamic parameters are already declared. For parameter estimation, in addition to a clear description of the mechanisms in form of ODEs, a reasonable number of fitting parameters must be specified for which the model equations are fitted to experimental data. It is worth mentioning that ReKinSim is a kinetic modeling tool and does not provide solutions to dynamics in a spatially resolved domain, meaning that processes like hydrodynamic flow and solute transport are exempted. Nevertheless, for inclusion of such transport mechanisms ReKinSim code-based module can be easily coupled with any transport code via operator splitting techniques (Valocchi and Malmstead, 1992).

ReKinSim is delivered in two modes: code-based and with a graphical user interface (GUI). The GUI mode facilitates ReKinSim application for the majority of users with little knowledge of programming. More advanced users can benefit further from the code-based mode that allows user to couple or integrate it into other codes.

Computer algorithms

The ReKinSim program is written in MATLAB. The code-based mode needs the MATLAB environment and the Optimization Toolbox to run. The GUI mode is stand-alone and for running only requires the proper MATLAB libraries which are patched with the program and are alternatively available for download on MathWorks website. The internal numerical solver is a multi-variable-step, variable-order solver based on the numerical differentiation formulas of orders 1 to 5, and therefore can handle relatively stiff problems (Shampine and Reichelt, 1997). The nonlinear regression method for fitting a model to data is based on least-squares minimization algorithms and can be chosen between trust-region-reflective or Levenberg-Marquardt. The trust-region-reflective algorithm is based on the interior-reflective Newton method described in Coleman and Li (1996). For this algorithm, the nonlinear system of equations cannot be under-determined, that is, the number of equations must be at least as many as the number of fitting parameters. On the other hand, the Levenberg-Marquardt algorithm cannot handle bound constraints, as described in Levenberg, 1944, Marquardt, 1963. As a consequence, under-determined problems with parameters that are constrained to a certain range cannot be solved with ReKinSim.

Definition of mechanisms in ReKinSim

A MATLAB function file allows the direct input of the kinetic mechanisms into ReKinSim in the form of a system of ODEs. The file thus contains the differential equations that are needed for fitting. The ReKinSim method of input is structurally different and more straightforward than other previously mentioned programs. For example, KinTek Explorer (Johnson et al., 2009) uses a simple text file with a specific syntax to input the mechanisms, from which the program then automatically derives the differential equations. As a result, KinTek Explorer is limited to the elementary reactions only and cannot cover the more complex environmentally related dynamics. Input of kinetic mechanisms into AQUASIM (Reichert, 1994) is more versatile, though it requires the error-prone procedure of determining the stoichiometry matrix from the system of differential equations.

Apart from trivial reaction mechanisms such as reversible and irreversible first, second or higher order reactions, ReKinSim's versatility in solving different, arbitrary types of ordinary equations allows the inclusion of more complex dynamics into the kinetic models, such as Bioavailability limitations (Best, 1955), isotope fractionation and enrichment (Elsner, 2010), absorption/desorption (e.g., Hill-Langmuir equation, Hill, 1910, Langmuir, 1918), biomass growth (e.g, logistic function, McKendrick and Pai, 1912), equations of chemostats/retentostats, etc. Therefore, its applicability goes beyond the examples shown later in Section 3, and ReKinSim can in principle handle any mechanism of interest as long as its description can be written in the form of a first-order differential equation,where are the species concentrations, are the kinetic parameters, and the function is a closed form mathematical expression linking them together.

Estimated standard errors

ReKinSim reports estimated parameter values with their respective standard deviation. The square roots of the diagonal elements of the covariance matrix, which is calculated from the Jacobian matrix, are used as estimators for the standard errors. The reported standard deviations can serve as a measure for the sensitivity of model results to the parameters and moreover can be used for uncertainty analysis. Sensitivity and uncertainty analyses are not included in ReKinSim.

Kinetic modeling scenarios

A few case studies were designed to test ReKinSim's accuracy and applicability. We first examined cases for which an analytical solution exists. Thereafter, scenario complexities were increased to an extent that they could only be solved numerically. For some of the scenarios, ReKinSim was used to fit the model to noisy data generated from a solution of the model for a set of parameter values, and then the estimated parameter values were compared with the actual ones. Some of the scenarios were chosen from literature and were previously presented as challenging tasks of kinetic modeling.

Enzymatic reaction

Enzymatic reactions are described by Michaelis–Menten kinetics (Michaelis and Menten, 1913)where is the maximum volumetric conversion rate and is the half-saturation constant. Eq. (1) is usually solved numerically, however, an analytical solution was recently found (Corless et al., 1996, Tang and Xiao, 2007) based on the Lambert function (Lambert, 1758, Lambert, 1772),where W is the Lambert function and is the initial concentration (at initial time ). Considering an enzymatic reaction with kinetic parameters and (taken from Harms and Zehnder (1994)), the solution of Eq. (1) with the initial concentration at the initial time is given as in Fig. 1. A noisy dataset was thereafter generated from this solution with the standard deviation equal to 15% of the actual concentrations at 20 time points. Subsequently, ReKinSim was used to estimate the kinetic parameters and by fitting Eq. (1) to this dataset. In a similar manner, the analytical solution Eq. (2) was also fitted to this dataset and the results were then compared to the fit with ReKinSim. The values of and , respectively, were estimated as and by the analytical fit, and as and by ReKinSim. The results show that the parameter values were estimated closer to the actual values by ReKinSim.

Fig. 1

Enzymatic reaction. Estimating the enzymatic kinetic parameters ( and ) by fitting both the Michaelis-Menten analytical solution Eq. (2) and the numerical form (Eq. (1) using ReKinSim) to the data. Both fits agreeably overlap with the actual solution from which the noisy data was generated.

Sequential first order decay

The analytical solution of a sequential degradation chain of finite intermediate products following first order irreversible decay (with equal stoichiometric yields) is given by (Cho, 1971):

The above Eq. (3) was extended by Höhener and Atteia (2014) to describe isotope ratios of intermediate products in the finite chains of degrading chemicals following a first-order kinetics. The proposed analytical formulas were then fitted as an example to the experimental data (containing both concentration and isotopic profiles) of a mixed initial contamination of PCE() and TCE () decaying to less chlorinated ethenes (DCE (), VC(), and ethene ()). Assuming that a specific chemical has two light and heavy isotopes in its structure (e.g., and ), the following ODEs can be written for light and heavy fractionswhere and are specific rate constants for light and heavy fractions, respectively (Elsner, 2010). α is the isotopic fractionation factor which links the change in concentration of isotopologues with the current concentrations according to . The ODE system of sequential first order decay with five components (Fig. 2) for both isotopologues of all components in respect to Eqs. (4a), (4b) is therefore written as below:

Fig. 2

Sequential first order decay. Reaction scheme of PCE (perchloroethylene) degradation to the less chlorinated ethenes.

ReKinSim was used to fit Eqs. (5a)–(5j) to the same experimental data as used in Höhener and Atteia (2014). Fig. 3 shows the solution to this example by ReKinSim, which is virtually identical to that from the analytical approach of Höhener and Atteia (2014). For this specific case, AQUASIM (Reichert, 1994, Reichert, 1998) was also used to calculate the values of rate constants and isotope fractionation factors. The values of kinetic parameters obtained from all three approaches were similar and are listed in Table 1.

Fig. 3

Sequential first order decay. The data points show concentration profiles and isotopic data of an example of PCE dechlorination (reaction scheme in Fig. 2) presented in Höhener and Atteia (2014). The lines show the ReKinSim fit to the data from which the parameter values in Table 1 are obtained.

Table 1

The estimated values of first order rate constants () and isotope fractionation factors () of the PCE decay chain (Fig. 2) by analytical solution of Höhener and Atteia (2014), ReKinSim (Fig. 3) and AQUASIM (Reichert, 1995).

	k_1	k_2	k_3	k_4	α_1	α_2	α_3	α_4
Höhener and Atteia (2014)	0.041	0.0695	0.0975	0.1917	0.9960	0.9903	0.9853	0.9810
ReKinSim (Fig. 3)	0.040	0.0698	0.0993	0.1934	0.9958	0.9901	0.9852	0.9808
AQUASIM (Reichert, 1995)	0.041	0.0696	0.0976	0.1921	0.9953	0.9889	0.9863	0.9792

The example shows the suitability of ReKinSim for datasets that contain not only concentration profiles but also isotopic data. This is particularly important since stable isotope data improve our understanding about the underlaying processes in an ecosystem and reduce parameter uncertainty (Adiyanti et al., 2016, Eckert et al., 2013, Qiu et al., 2013, Thullner et al., 2012, Van Engeland et al., 2012). With few exceptions (e.g., AQUASIM), most of the currently available, similarly-oriented programs have difficulties in providing a platform in which isotope data analysis can contribute as another compelling measure to data fitting and parameter estimation.

Thermal isomerization of α-pinene

α-Pinene, a primary component of turpentine, is the most widely encountered terpenoid in nature and is well-known for its natural insect repellent characteristics. The thermal isomerization of α-pinene includes two irreversible reactions from which Dipentene and Alloocimene are formed. Subsequently, Alloocimene transforms irreversibly to Pyronene or forms a dimer which can reversely break into Alloocimene again (Fuguitt and Hawkins, 1945). The diagram of the above reactions is shown in Fig. 4. In previous nonlinear regression analyses (e.g., Bates and Watts, 1988, Box et al., 1973), this was considered a rather difficult test case for fitting a complex kinetic scheme to data using analytical methods and was revisited as a test case for the program Global KinTek Explorer (Johnson et al., 2009). As discussed in Johnson et al. (2009), in previous studies the forward reaction of Alloocimene leading to formation of a dimer was repeatedly, and incorrectly, assumed to be unimolecular. Since two Alloocimene molecules form a dimer (Fuguitt and Hawkins, 1945), the forward reaction must be considered bimolecular. On the contrary, the backward reaction remains unimolecular. The system of ODEs for this reaction scheme (Fig. 4) is written as:

Fig. 4

Thermal isomerization of α-pinene. The reaction scheme.

Through fitting the data to differential equations of the corrected reaction scheme (Eqs. (6a), (6b), (6c), (6d), (6e)) using ReKinSim, the values of rate constants were estimated close to those obtained from KinTek Explorer program (Johnson et al., 2009) and are compared in Table 2. The values for and are slightly different between the two programs. This might be due to the differences between the two codes such as the internal numerical schemes, relative tolerances, initial starting points, or the definition of upper and lower bounds for each parameters. This might also be due to the insensitivity of the results of this specific kinetic model (Eqs. (6a), (6b), (6c), (6d), (6e)) to the above parameters, which therefore makes it possible to still obtain a good fit even with a slightly bigger variation of these values. The absolute-relative sensitivity of the model results to all the kinetic parameters was calculated for this specific test case. As shown in the sensitivity analysis results (Fig. 6 (2–6)), the concentrations of α-pinene and Dipentene were only sensitive to the values of and . Alloocimene concentrations were negligibly sensitive to and Pyronene concentrations showed a low sensitivity to . The concentrations of dimer were however insensitive to both and . This led to a lower sensitivity of the model results as a whole to and values, which is well mirrored by a larger uncertainty ranges of these parameters (shown as standard deviations in Table 2).

Table 2

Rate constants values of thermal isomerization of -pinene (reaction diagram is shown in Fig. 4) estimated by the KinTek Explorer program (Johnson et al., 2009) and ReKinSim (Fig. 5). It is worth noting the dimension of being as per minute per concentration of Alloocimene in the mixture (as percentage) due to the bimolecular transformation of two Alloocimene into a dimer.

×10^5	k_1[min^−1]	k_2[min^−1]	k_3[min^−1]	k_4[min^−1%Allo−1]	k_5[min^−1]
KinTek Explorer	6.0	3.0	1.6	2.5	5
ReKinSim	5.9±0.005	3.0±0.005	2.1±0.029	2.6±0.04	3.6±0.134

Fig. 6

Thermal isomerization of α-pinene: sensitivity analysis of the model (Eqs. (6a), (6b), (6c), (6d), (6e)) results. (1) Shows Fig. 5 with the estimated error bounds calculated from uncorrelated linearized error propagation of change of parameter values equal to 1% of their respective standard deviation. Note that the y-axis is scaled logarithmically to better illustrate the approximate standard deviation of the model results. (2–6) Show absolute-relative sensitivity functions of the calculated concentrations with respect to five model parameters.

Fig. 5

Thermal isomerization of α-pinene. The data are obtained from Fuguitt and Hawkins (1945) and the lines show the fit of the reaction scheme in Fig. 4 (Eqs. (6a), (6b), (6c), (6d), (6e)) by ReKinSim. The figure is similar to the Fig. 8 of the Johnson et al. (2009) where a fit by the program KinTek Explorer is demonstrated. The estimated rate constant values by both programs are listed in Table 2.

The data for this example are experimentally driven, therefore the absolute values of these parameters are not precisely known and it is difficult to express which method resulted in a better match. That said, the norm of residuals for our fit with ReKinSim is about 4.29, which was smaller than the norm of residuals of the fit by Kinetic Explorer program (4.67). A lower value of the norm of residuals indicates that a better fit, and thus a more precise set of parameter values, was attained by ReKinSim. Fig. 5 shows the ReKinSim fit of the model to the isomerization data. In a different attempt to fit the data to the erroneous model for which dimerization was assumed as a first-order reaction, ReKinSim confirmed the formerly reported rate constants (data not shown).

A hypothetical test case

For the hypothetical test case, we assumed a reaction network of the following:with six components and three different reaction types. , and are the rate constants for the Reactions (7b), (7c), (7d). The unimolecular rate constant () has the dimension while the rate constants of bimolecular reactions ( and ) have the dimension . For the Michaelis-Menten reaction (7a) and are respectively half-saturation constant and maximum conversion rate. The ODE system for the above reaction network (7) is as follows:and was solved with the following initial concentrations: , , , , , , for the actual values of the kinetic parameters: , , , , . The goal was to back-estimate the kinetic parameters from noisy observations which are artificially created from the model solution by adding normally distributed pseudorandom errors. The added random noise was independently generated from a normal distribution with specific mean equal to the model solution at 40 time points (for each component) and standard deviation of 10% of the mean value. Therefore, for each component 40 randomly generated data points were obtained, which in total led to 240 data points for six components. Using ReKinSim, the kinetic model Eqs. (8a)–(8f) were then fitted to these data points to back-estimate the rate constant values. The actual solution together with the noisy data points and the fitted solution by ReKinSim are illustrated in Fig. 7. Despite the seemingly high complexity of this particular model which included three different types of reaction, the estimated solution by ReKinSim correlated well with the actual solution from which the noisy data points were generated. ReKinSim estimated rate constants and their respective standard deviation as following: , , , , and .

Fig. 7

The hypothetical scenario. The tick dashed lines show the actual solution from which the data points were generated by applying a normally distributed random noise. The solid lines show the fit of the kinetic model Eqs. (8a)–(8f) to the data using ReKinSim. Note that only 20 points out of actual 40 data points for each component are shown.

Discussion

Comparison with some currently available kinetic modeling tools

The ReKinSim program is easy to operate and in comparison to other similar programs provides users with a number of advantages, such as easier management of the system of ODEs, better control on the minimization process and minimizing function, flexible data import/export from/to other programming environments, and the chance for post-processing analyses in the popular MATLAB environment. Furthermore, the use of MATLAB's sophisticated functions equips ReKinSim with capabilities such as superior handling of numerical instabilities, easier detection of errors, and the options to choose the algorithm of minimization. In comparison to some of its counterparts (e.g., Dynafit (Kuzmič, 1996), COPASI (Hoops et al., 2006), and KinTek Explorer (Johnson et al., 2009)), ReKinSim is more of a general-purpose program that can take over a wider range of applications and is not limited to few predefined sets of elementary or enzymatic reactions. It therefore gives the user the freedom of solving any system of ODE that might be relevant for an environmental system. For instance, environmentally related processes and observables such as bioavailability limitations, isotopic data, and more complicated mechanisms than elementary reactions (e.g., Michaelis-Menten and inhibition) can be easily included in ReKinSim. Almost all the ReKinSim alternatives are written as standalone programs which are difficult to connect to other environments and sources. Since ReKinSim is written in MATLAB, it can be easily linked to other programming environments (e.g. C, C++, Fortran, Python) through the MATLAB application programming interface (API). ReKinSim has the additional advantage of being non-platform-dependent. The code-based mode can run on any operating system (Windows, Linux or OSX) where MATLAB is installed. The GUI mode is also available for any of the mentioned operating systems. Apart from that, ReKinSim is not a commercial program and is available for free, unlike most of its equivalents (e.g., Dynafit and KinTek Explorer).

In terms of application and versatility, ReKinSim is similar to the program AQUASIM (Reichert, 1994, Reichert, 1995). Several advantages of AQUASIM over ReKinSim are: no dependency on other environments or libraries and thus no requirement of installing extra packages; extra modules for parameter sensitivity analysis and uncertainty analysis; and the solution for a simple linear reactive transport system. As for the drawbacks, because AQUASIM is a standalone program, it causes difficulties and inconsistencies when AQUASIM is required to communicate with or be coupled to other programs (e.g., external transport codes). In addition, we found the procedure of entering a kinetic model into AQUASIM, which involves calculating the stoichiometric matrix and inserting the right coefficients at the right input windows, particularly difficult. Although a matter of preference, we believe characterizing the kinetic model via the system of ODEs and addressing it directly into ReKinSim is more straightforward and less error prone. Since ReKinSim is written in MATLAB and the changes are directly applied by modifying its source code, it gives more flexibility for modifications than AQUASIM. It is also easier for ReKinSim to communicate with other environments, for example to import large datasets into ReKinSim from other sources (e.g., Excel or ASCII files), or to export the results back to them. Furthermore, it is relatively easy and computationally inexpensive to couple ReKinSim code-based mode with other programs, in particular to the recent scientific codes which are mostly written in MATLAB (e.g., Gharasoo et al., 2012, Gharasoo et al., 2014) or in languages that are easy to link with MATLAB. The direct linking of programs considerably reduces computational time, amongst many other benefits such as improving code consistency, increasing code readability, and simplifying error handling.

AQUASIM gives users the choice between two minimization algorithms, the downhill simplex method (Nelder and Mead, 1965) and the secant method (Press et al., 1992, Ralston and Jennrich, 1978). Both of these methods are slightly inferior in comparison to more advanced algorithms used by ReKinSim - trust-region-reflective (Coleman and Li, 1996) and Marquardt-Levenberg (Levenberg, 1944, Marquardt, 1963) - mainly due to the fact that the latter algorithms use derivative information to converge and are equipped into MATLAB with the latest up-to-date improvements.

ReKinSim usability as an external reactive module

ReKinSim has two numerical modules, kinetic simulation and parameter estimation. When the rate constants and the initial concentrations are known, a kinetic model can be solved by the kinetic simulation module within a finite time interval similar to AQUASIM (Reichert, 1994, Reichert, 1995) or Biological Reaction Network Simulator - BRNS (Regnier et al., 2002). The results then show the change of concentrations in a zero dimensional system over time. ReKinSim's kinetic simulation module can be easily coupled through operator splitting techniques with a simple linear finite difference transport code to give a similar functionality as the Reactor Network in AQUASIM, or the reactive 1D longitudinal transport in BRNS. By coupling this module to more advanced transport codes such as those published by Gharasoo et al., 2012, Harbaugh et al., 2000 and Zheng (2006), it is possible to simulate reactive transport of contaminants inside more complex environments. Given that ReKinSim can solve a variety of ODEs and is not limited to elementary reactions, it can be easily adapted to solve environmentally related equations such as mass transfer limitations (Best, 1955, Thullner et al., 2008), inhibition (Gharasoo, 2014, Gharasoo et al., 2015), and adsorption/desorption (Goldberg et al., 2007). In comparison to other linkable kinetic simulation codes (e.g., Farley et al., 2008, Kneis et al., 2017, Morrill et al., 2006, Regnier et al., 2002), coupling of the ReKinSim kinetic simulator proved more flexible and technically less challenging. For instance, the coupling of BRNS to a pore network model (Gharasoo et al., 2012), or to OGS (formerly known as GeoSys, Centler et al., 2010), required compiling the code to create a dynamic-link library file (.dll) which then needed to be loaded in the new environment and called through a complicated set of Windows API commands. This is not only an elaborate process, but also a common source of errors.

Computational runtimes

ReKinSim runtimes highly depended on the model complexity, number of components and mechanisms, parameters initial guess, absolute and relative tolerances, amount of data points, and required accuracy of fit. For the presented test case scenarios, on a Linux system with Intel Core i5-4590 CPU at 3.30 GHz and 8 GB RAM, the fits were achieved in fractions of a minute suggesting that for cases of similar complexity the computational times are quite reasonable. For the second scenario (Section 3.2), ReKinSim runtime was of the same order of magnitude as was required by AQUASIM (a few seconds).

Summary and conclusions

The ReKinSim program can smoothly fit experimental data to various environmentally-related mechanisms. Though powerful, ReKinSim is easy to operate, allowing users to concentrate more on data interpretation than on the otherwise complicated mathematical procedures of data fitting. ReKinSim is also a very effective tool for learning and understanding the dynamics of chemical reaction processes. It provides insights to the mechanistic understanding of reaction mechanisms inside an environmental system, and lets the user predict the behavior of an investigated system and interpret the changes of concentration profiles under different circumstances. It thus offers an invaluable simulation platform for planning ahead before performing experiments and quantitatively explaining experimental data.