# Numerical Simulation of Non-Linear Models of Reaction--Diffusion for a DGT Sensor
## 1 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]. For instance, if the ML species schematized in Figure 1 reacts as follows:  
\[\mathrm{M}+\mathrm{L}\ \xrightleftharpoons[k_{\mathrm{a,1}}]{k_{\mathrm{d,1}}}\ \ \mathrm{ML} \tag{1}\]  
and  
\[\mathrm{ML}+\mathrm{L}\ \xrightleftharpoons[k_{\mathrm{a,2}}]{k_{\mathrm{d,2}}}\ \ \mathrm{ML}_{2}, \tag{2}\]  
the change of ML concentration (\(c_{\mathrm{ML}}\)) with time will be given by [3]  
\[\frac{\partial c_{\mathrm{ML}}}{\partial t}=D_{\mathrm{ML}}\frac{\partial^{2} c_{\mathrm{ML}}}{\partial x^{2}}-k_{\mathrm{a,1}}c_{\mathrm{ML}}c_{\mathrm{M}}+k_{ \mathrm{d,1}}c_{\mathrm{ML}}-k_{\mathrm{a,2}}c_{\mathrm{ML}}c_{\mathrm{L}}+k_{ \mathrm{d,2}}c_{\mathrm{ML}_{2}}, \tag{3}\]  
where \(c_{i}\) is the concentration of specie \(i\), \(t\) stands for the time, \(x\) represents the spatial dimension, \(D_{\mathrm{ML}}\) is the diffusion coefficient of the ML species, and \(k_{\mathrm{a,i}}\), \(k_{\mathrm{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=\mathrm{M}\), \(\mathrm{L}\), \(\mathrm{ML}\), \(\mathrm{ML}_{2}\),...), the variation of the concentration with time will depend on the diffusive term plus the reaction terms, as follows:  
\[\frac{\partial c_{i}}{\partial t}=D_{i}\frac{\partial^{2}c_{i}}{\partial x^{2 }}+\text{*reaction terms}^{\gg}. \tag{4}\]  
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  
Figure 1: Schematic representation of coupled reaction–diffusion processes of \(\mathrm{ML}_{i}\) complexes, when the metal is adsorbed at an interface. Here, \(D_{j}=\) is the diffusion coefficient of the species \(j\) (\(j=\mathrm{M}\), \(\mathrm{L}\), \(\mathrm{ML}\), \(\mathrm{ML}_{2}\),...), and \(K_{i}\), \(k_{\mathrm{a,i}}\), \(k_{\mathrm{d,i}}\) represent the equilibrium, the association and dissociation constants of \(\mathrm{ML}_{i}\), respectively. 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.