Numerical Simulation of Non-Linear Models of Reaction—Diffusion for a DGT Sensor

: In this work, we present a novel strategy for the numerical solution of a coupled system of partial differential equations that describe reaction–diffusion processes of a mixture of metals and ligands that can be absorbed by a sensor or a microorganism, in an aqueous medium. The novelty introduced in this work consisted of an adequate database management in conjunction with a direct iterative schema, which allowed the construction of simple, fast and efﬁcient algorithms. Except in really adverse conditions, the calculation is converging and satisfactory solutions were reached. Computing times showed to be better than those obtained with some commercial programs. Although we concentrate on the solution for a particular system (Diffusive Gradients in Thin Films [DGT] sensors), the proposed algorithm does not require major modiﬁcations to consider new theoretical or experimental conﬁgurations. Since the quality of numerical simulations of reaction–diffusion problems often faces some drawbacks as the values of reaction rate constants increase, some additional effort has been invested in obtaining proper solutions in those cases.


Introduction
The high risk of contamination due to the presence of heavy metals in aqueous systems requires the development of reliable analytical techniques, capable of measuring the metal flux that reaches microorganisms, algae, plants or other live organisms present in the media [1].
Generally, in aquatic ecosystems, the metal flux toward an interface where it is consumed (for instance, by an organism or an analytical sensor), is the result of diffusion and kinetics of interconversion between metal and the various species present in the system (see scheme outlined in Figure 1).
A quantitative assessment of the metal flux that reaches a consuming interface requires the rigorous solution of the resulting system of transport and reaction equations, taking into account the geometry and the temporal and spatial scales. The mathematical models that describe these reaction-diffusion processes consist of systems of partial differential equations (PDE) with polynomial reaction terms [2].
the change of ML concentration (c ML ) with time will be given by [3] ∂c ML ∂t = D ML ∂ 2 c ML ∂x 2 − k a,1 c ML c M + k d,1 c ML − k a,2 c ML c L + k d,2 c ML 2 , where c i is the concentration of specie i, t stands for the time, x represents the spatial dimension, D ML is the diffusion coefficient of the ML species, and k a,i , k d,i (i = 1, 2) represent the association and the dissociation rate constants given in reactions (1) and (2). In general, for any species i in the media (i = M, L, ML, ML 2 , . . . ), the variation of the concentration with time will depend on the diffusive term plus the reaction terms, as follows: Thus, the system will contain as many equations as species, coupled by the reaction terms. However, in most practical cases, the resulting system of PDE cannot be analytically solved. In these cases, numerical approximations will be the only way to reach solutions for these systems.
The most common form of solving these systems consists of, as a first step, the discretization of equations, which transforms the differential system in an algebraic system different for each increase of time. Next, when the resulting system is non-linear, it could be solved by Newton-type techniques [4].
In this second step, two difficulties can be mentioned. The lower one is the lack of robustness of the method, which requires the introduction of Armijo or norm minimization strategies [5], to improve the convergence of the iterative process. However, the greatest obstacle appears as a consequence of the dimension of the linear system (n). In these problems, the number of unknown variables, n, is equal to the product of the number of species by the number of nodes in the spatial discretization. Considering that n may be large and that resolution time of a linear equations system is of the order of n 3 /3 s [5], the computational cost can be unaffordable in a general-purpose computer.
In this work, a set of algorithms that allow simulating systems with a large number of chemical species, and with a computational cost that supports the simulation of processes of some complexity in a reasonably short time, are presented. For this, Section 2 introduces the theoretical framework and the mathematical model. Section 3 discusses the iterative resolution of the system. Section 4 shows the results and analysis. The Appendix A contains the general resolution algorithm.

The Model
The proposed algorithm was developed to solve the reaction-diffusion equations that model the operation of a Diffusive Gradients in Thin Films (DGT) sensor [6], submerged in an aqueous solution containing different metal species. Nevertheless, the algorithm does not require major modifications to consider new theoretical or experimental configurations.

Diffusive Gradients in Thin Films (DGT)
The analytical technique known as DGT consists of simple devices developed for the in-situ measurement of bioavailability of heavy metals in soils and aquatic environments [6,7].
The DGT sensors are made up of three layers (see Figure 2): a hydrogel layer impregnated with a complexing material (resin layer), a hydrogel used as a diffusion layer, and a membrane filter. These three layers are placed in a plastic holder device [7].  DGT sensors were designed to emulate the metal absorption process in aqueous media. When a DGT is deployed in a sampling solution, where all species are in equilibrium, the metal species diffuse through the diffusive gel (where they might associate or dissociate), toward the resin layer where they accumulate as metal resin complexes.
As long as the resin does not get saturated, a concentration gradient is maintained in the diffusion layer. Very soon, (quasi) steady-state conditions are reached, so the initial transient can be neglected to obtain easy approximate analytical expressions. Fick's first law of diffusion describes the flux (J i ) for each species i, where D i is the diffusion coefficient, c i is the concentration and dc i /dx is the gradient of spatial concentration of the i species.
After the DGT has been exposed to the solution, it must be opened and the analytes eluted from the resin to be quantified.
When only metal is present in the solution, the accumulated mass in DGT is related to the solution concentration through a simple equation that neglects the effects of dynamic speciation [6]: where n stands for the accumulated moles, D represents an average diffusion coefficient, c * M is the metal concentration at the solution, A is the DGT area, t is the deployment time, and g is the diffusion layer thickness.
In the presence of metal complexes, c * M is replaced with c * DGT , which indicates the effective or apparent metal concentration leading to this accumulation. c * DGT is usually called the labile metal concentration of the solution.
The rigorous interpretation of c * DGT in terms of the real species in solutions requires taking into account all reaction-diffusion processes involved in the system, which can only be achieved via numerical simulation.

Mathematical Model
Consider the complexation of different metals and ligands in solution, that react according to the elementary schema represented in the Equation (7): where A and B can represent a metal, a ligand, or any combination of them, and k a and k d represent the kinetic constants of association and dissociation, respectively. The equilibrium conditions for this reaction are expressed in the form where K is the equilibrium constant, and c * i indicates the concentration of species i in the solution (bulk concentration).
When a DGT sensor is deployed in solution, the species diffuse to the resin (located between 0 < x < r), through the gel layer between r < x < r + g (see Figure 3).
x [mm] resin diffusive gel solution c r 0 g Inside the DGT, the concentration c(x, t) of each species is modified by diffusion, reaction with other species present in the solution, and by the possible reaction with the resin.
Before presenting the mathematical formulation, we will make some assumptions.
• Since the resin discs are made of complexing material embedded in the same type of gel that composes the diffusive layer, it is reasonable to consider that species that diffuse through the gel can penetrate the resin [8].

•
Species that diffuse through the gel can bind the resin after. They do it by the described schema in Equation (9): where R denotes the free-sites in the interior of the resin, and k a,R , k d,R and K R = k a,R /k d,R represent the corresponding kinetic and stability constants.

•
When the material used for the resin is Chelex, the grains of resin are located mainly in the layer adjacent to the diffusion gel [9]. In a first approach, we consider that binding sites in the resin are uniformly distributed.

•
Because diffusion coefficients in the gel, the filter and diffusive boundary layer (DBL) are similar [10], the diffusive layer width g in Figure 3 includes the size of all these layers.

•
Migration forces or any electrostatic effects due to the charge of resin are not considered, i.e., the background salt is enough to screen these interactions.
Under these conditions, diffusion is the most relevant transport mechanism in the system and, therefore, the transport equations for each species i can be written as where D R i and D i are diffusion coefficients of the species i in the resin and gel, respectively, and the so-called «terms of reaction» refer to the kinetic terms corresponding to the reaction of the specie i with other substances of the solution.
Technically, the expressions shown in Equation (10) are second-order parabolic partial differential equations that can be analytically solved only under very specific conditions.
In general, in the cases of practical interest, these equations do not have analytical solutions. Thus, numerical solutions may be the best way to get the concentration values for each point of space in every instant of time.
As described in the introduction, the most widely used approach to numerically solve these systems consists of discretizing the equations and, subsequently, iteratively solve the resulting system of algebraic equations. The following sections describe a novel strategy to numerically solve the systems of equations of the form (10).

The Problem
Let us consider a set of elementary reactions, in a number equal to numreac, of the form Reaction 1 , where k a,j , k d,j and K j = k a,j /k d,j , for j = 1, 2, . . . , numreac, represent the kinetic and stability constants, respectively. It is worth mentioning that in the set of reactions (11) some species could be repeated, e.g., X L 11 might be equal to X L 21 The system is encoded according to the indexing used in Equation (11), without any restrictions on the number of metals or ligands. As it has already been said, the only restriction is that all chemical reactions should be written in the form A + B C. Given that starting data are the total concentrations of the pure substances (i.e., those that are not the result of a reaction), it is first necessary to solve the equilibrium problem, that is to say, it is needed to find the concentration value in solution (or bulk concentration) for each of the species in the system.
To do this we build mreac and kreac, arrays that allow introducing the information related to the chemical reactions that must be considered in x = r + g. The first of these arrays, mreac, is the index array that contains three columns and a number of rows equal to the number of numreac reactions, The mreac array is constructed in such a way that, if in the j-th reaction the substance 5 is combined with the 7 to produce 8, the indices of the row j will be 5 7 8; that is, the first two columns contain the components and the third indicates the formed compound.
The kreac matrix contains the same number of rows as the mreac matrix, but only two columns. The first column indicates the value of the association constant (k a,j ), and the second one the value of the dissociation constant (k d,j ) of each j reaction, with j = 1, 2, . . . , numreac, The strategy to solve the problem is to express the non-pure of species (i.e., species that are listed in the third column of mreac, because they are the result of a reaction), as a combination of the total concentrations of the pure substances, which are known for the system. The result is a system of nonlinear equations, which is solved using a variation of the Newton-Raphson method.
Once the equilibrium problem is solved, that is to say, once the bulk concentration values are obtained for each species at t = 0 and x = r + g, it is necessary to pose the problem of transport for an arbitrary number of chemical species.
Let numeq be the number of species of the problem, which, of course, coincides with the number of equations. In the first place two vectors of coefficients are considered, one for the region (0; r), hereinafter referred to as: , which contains the diffusion coefficient for each species inside the resin phase, and the other denoted by D = (D 1 , D 2 , . . . , D numeq ), which corresponds to diffusion coefficients inside the gel (r < x < r + g). Diffusion coefficients of the resin and its complexes are zero since they are considered to be fixed in the gel matrix of the resin disc.
The mreac and kreac matrices, that were used for solving the equilibrium problem, must be extended incorporating the reactions corresponding to the resin and its complexes.
These two new matrices contain the information to express the general problem. For each species i, where i = 1, . . . , numeq, the corresponding reaction-diffusion equation will be as: where D R i and D i act as a diffusion coefficients in (0, r) and (r, r + g), respectively. In Equation (14), for a given value of index i (which runs across chemical species), the coefficients A i,j and B i,j can be calculated using the following algorithm: if in the reaction j does not participate the species i, A i,j = B i,j = 0. Otherwise, the absolute value of A i,j is equal to k a,j and the absolute value of B i,j is equal to k d,j . Finally, it must be taken into account that the coefficient (A i,j or B i,j ) of the term that contains the substance i-th must be negative.

Initial and Boundary Conditions
When the DGT is deployed in the solution, at time t = 0 the concentration of each of the species i inside the spatial domain (0, r + g) is zero, that is c i (x, 0) = 0, except for the free places inside the resin, whose initial concentration is c R (x, 0) = c TR (the total concentration of resin sites).
Resin sites are considered immobile, but for species with D i = 0, it must be considered one boundary condition at x = r + g, i.e., at the interface between the diffusive gel and the solution: where c * i stands for the bulk concentration of species i. At the DGT bottom, i.e., at x = 0, there should be no flow of any species, Finally, at the resin-gel interface (x = r), there must be continuity in the concentration profiles and flows of each species,

Dimensionless Formulation
To adapt scales to the order of magnitude of data, the variable changes that lead to the nondimensional problem, are described.

•
For all species in solution, that is, for all the species that have a bulk solution concentration, obtained from the solution of the equilibrium problem, we define standard concentrations, • For the rest, that is to say, the resin and its complexes, • For the spatial coordinate, instead of nondimensionalisation, a change of scale is used, where D max is the highest diffusion coefficient. It makes it necessary to define r * = r/ √ D max ; Given that the aim is to simulate experiments involving different time scales, the time axis is not changed.
With these new variables, the reaction-diffusion equations described in Equation (14) are rewritten as follows: (Â i,j q mreac(j,1) q mreac(j,2) +B i,j q mreac(j,3) ) , where d R i acts as a diffusion coefficient in(0, r * ) and d i in (r * , r * + g * ), and where the new coefficientsÂ i,j andB i,j are calculated asÂ For the concentrations of the resin and its complexes, c * i must be replaced by c TR . The initial and boundary conditions for all substances that do not contain resin, that is to say, for all those whose diffusion coefficient is different from zero, are now and The remaining substances only need an initial condition: for the resin For complexes of the resin q i (x, 0) = 0 in (0, r * + g * ).

Discretization
The spatial discretization of a PDE model can be performed using different techniques. The most frequently used are those of finite differences, finite element (with all its variants) or finite volumes [4,11]. Each of these techniques has its advantages and disadvantages. For instance, discretizing by the finite difference method leads to simple computer programming, as long as fixed spatial steps are used. However, this latter could imply the use of a large number of points and, therefore, longer calculation times.
On the other hand, the finite element method allows working with arbitrary point distributions, using smaller spatial steps in regions where it is assumed that function could rapidly changes, so error in the iteration could be large, and larger widths where it is known that the error will be smaller [12]. Although this flexibility results in a computational cost lower than in the case of the finite differences, it involves an algorithmic and programming effort that can be significantly high.
In this work, the finite differences method was used, taking fixed spatial steps of length ∆z, and temporal intervals equal to ∆t (see Figure 4). Figure 4. Scheme of space and time discretization.
Let F(z, t) be a function that can contain spatial derivatives (of different orders) of the function f (z, t). If Equation (19) can be discretized as a linear combination of the form where, 0 ≤ λ ≤ 1.
In general, for the PDE of parabolic type, the backward Euler method for discretization (i.e., taking λ = 1 in Equation (20)), allows choosing arbitrary increases in z and t, without risking its stability or introducing spurious oscillations in the iterative process [4]. Next, for every substance i, the differential Equation (15) will by discretized using the backward Euler scheme to obtain a system of algebraic equations.
For species with diffusion coefficient different from zero, the Equation (15) is transformed into In this case, the described boundary conditions in Equations (17) and (18) will be q i (∆z, t + ∆t) − q i (0, t + ∆t) = 0 , and There is no need for any expression for the solution continuity at z = r * . Simply, from now on, it will be considered for such a position, that each substance has a unique concentration value (a simple way to state that it is continuous, due to the absence of electrostatic effects).
Due to the fact that the resin and its components are considered motionless, the differential equations for these substances do not contain spatial derivatives. In this way, the equations described in (15) for species with zero coefficient diffusion are re-written as follows: for (0, r * ), and q i (z, t + ∆t) = 0, ∀ z ∈ (r * , r * + g * ) .

Iterative Resolution of the Resulting System of Equations
In Equations (21) and (24), the terms derived from reactions couple the resulting algebraic equation system, by involving second-order terms as products of concentrations from different species. As mentioned in the introduction, these types of non-linear systems are generally solved by using iterative-type Newton techniques. The novelty introduced in this work consisted in convenient database management operated in conjunction with a direct iterative schema, which allowed us the construction of simple, fast and efficient algorithms for solving reaction-diffusion problems.
The basic idea of these algorithms consists on the assumption that the solution for all species is known at time t; then, the solution at time t + ∆t is obtained according to the following scheme: Initialize all concentrations at t + ∆t with the solution at time t.

2.
for i = 1, numeq To solve discrete i-th equation (for the species i-th) considering all other j-th species (i = j) constant in this stage and whose value at t + ∆t being the corresponding to its most recent calculation. endfor 3.
Repeat step 2 until the difference between concentrations in two consecutive iterations be less than a small fixed value.
The described iterative structure has the advantage that resolution, for each substance i with diffusion coefficient different from zero, is reduced to a tridiagonal system in which the diagonal is strongly dominant, due to the negative sign of the terms that contain c i . Even better, for the resin and its complexes, species for which the diffusion coefficient is zero, the system is diagonal.
All this results in a great increase on calculating speed, since it is better to solve a dozen or hundred times these systems (each one of them solved in a time proportional to n) than to solve the complete non-linear system in Equations (21) and (24). In the latter case, even if the system is solved with a small number of iterations, each iteration involves a calculation time of n 3 order. Bearing in mind that n equals to the number of species times the number of intervals in the spatial domain, any attempt to directly solve the system can be prohibitive.
Implicitly, when describing the discretizations, it has been considered that the complete domain (0, r * + g * ) is discretized in intervals of equal length ∆z in both subdomains. Then, for each species, different vectors that represent the solution at each position of the mesh, at each time, and for a given k state of the iterative process, will be used.
If the concentrations vector for a species is represented by where nint indicates the number of intervals in which the domain is divided (see Figure 5), the vectors used in the algorithms will result from the juxtaposition of vectors for each species at different times and iterative levels. Essentially, these will be used: q aux ≡( q k 1 (t + ∆t), q k 2 (t + ∆t), . . . , q k numeq (t + ∆t)) q nex ≡( q k+1 1 (t + ∆t), q k+1 2 (t + ∆t), . . . , q k+1 numeq (t + ∆t)) .
Given that solution at z = r * + g * is known, its value is not included within the resolution. Then, the index point l is on the left end from the subinterval l, that is to say, it corresponds to z = (l − 1)∆z position. With index nr the corresponding last interval of the resin subdomain is indicated, so z = r * corresponds to the nr + 1 index. Figure 5. Scheme of space discretization.
The length of the vectors q pre , q aux and q nex is equal to ntot = numeq × nint. All this implies that to have access to the j-th spatial position of a vector corresponding to the i-th substance, it must be calculated l = (i − 1) * nint + j. This relationship is used several times within the general resolution algorithm detailed in the Appendix A.
It is worth noting that although the complete domain (0, r * + g * ) has been divided into equal-length intervals ∆z in every subdomain, this discretization is not essential. In fact, if required, the problem can be easily split into several subdomains, each with its own length.
This provides two major advantages: on the one hand, it allows choosing small intervals where functions quickly vary (usually near the interface) and large intervals in the rest of the domain, increasing calculation speed even more. On the other hand, it allows solving configurations different from those of the conventional DGT sensor (as the case of the voltammetry detectors (r → 0) [13], or DGT sensors with several resins [14,15]). That is to say, with the condition that the reactions may be expressed in the basic form A + B C, the proposed algorithm can be adapted to solve different arrangements, simply by modifying the initial and boundary conditions.
Generally, the quality of numerical simulations of reaction-diffusion problems often faces some drawbacks as the values of reactions rate constants increase. Nevertheless, if the kinetics rate constants increase, while keeping their quotient constant (i.e., K = k a /k d =const), the concentration profiles tend to be those that would result in the fully labile case, i.e., the case in which the rate constants are infinite but the stability constant, K, remains fixed. Thus, the use of the solution of the fully labile case could provide good initialization values for the problem with large kinetic rate constants. Unfortunately, it is not possible to obtain always the fully labile solution in the general case. For this reason, specific techniques have been developed to initialize in these cases.
The general strategy is to solve the system, at each time, for a set of kinetic constants lower than desired but that provide a solution. With this result, the problem is initialized for a larger value of the constants of reaction (keeping the same K value). The procedure is repeated until the values of the require constants are reached. The described algorithm ceases to operate when the constants are enormous, but even so, it allows to explore a range of sufficiently large values to obtain conclusions about the behavior of the systems of experimental interest.

Results and Analysis
Initially, the algorithm was implemented in FORTRAN and used to calculate the concentration profiles in a DGT, when it was introduced into a solution containing a metal (M) and a ligand (L) that react according to the scheme (28) Figure 6 shows the normalized concentration profiles for metal (q M ) and complex (q ML ) in the DGT numerically computed after 4-h of deployment. Numerical results are in good agreement with the experimental accumulation of Cd by DGT sensors in Cd-NTA systems [3]. Also, our results are in full agreement with the data obtained by simulating the same physicochemical system using commercial COMSOL Multiphysics R software (COMSOL, Inc., Los Altos, CA, USA).
Besides, systems containing a mixture of metal and various ligands were successfully simulated [16], showing a higher computer speed than commercial software. Table 1 shows a comparison in computing times between our program and the commercial COMSOL Multiphysics R software, for the case of a DGT that was submerged during 4-h in solution containing a metal (M) and one or more ligands (L).  Although the finite differences method was used, the method adopted for spatial discretization is not essential in terms of the work shown. As it was already explained, the suggested algorithmic novelty consists of proper database management, parallel to the iterative resolution of the non-linear system, solving species one by one for every time. Thus, for a standard simulation involving the resolution of hundreds of tri-diagonal systems, the computational cost is proportional to the number of points n. Instead, trying to solve a few systems of complete linear equations, whose dimension equals to n times the number of species (numeq), the computational cost will be proportional to the cube, not of n, but of n × numeq. The gain in runtime is clear; it is enough to consider n = 1000, whatever the number of involved chemical species, to realize that. Certainly, the finite elements method would increase a little more the speed of the simulations, so the possibility of its use should not be rejected. However, the proposed iterative scheme has proved to be efficient even in the worst possible scenarios, as is the use of finite difference method.
Using this algorithm, analysis of rigorous digital simulation of the diffusion-reaction processes in experimental accumulation by DGT sensors in several systems led to important theoretical developments related to DGT functioning [3,[17][18][19].

Conclusions
To be realistic, and bearing in mind the more physicochemical than mathematician feature of this work, it should be noted that a rigorous study of the convergence of the iterative method was not performed. Despite this, it was ensured that, except in really adverse conditions, the calculation is converging and, in fact, satisfactory solutions were reached. In general, these adverse conditions tend to be the result of poor initializations and, as it was as described above, some additional effort has been invested in obtaining proper initializations.