Computational Efficiency of Decoupling Approach in Solving Reactive Transport Model : A Case Study of Pyrite Oxidative Dissolution

Pyrite existed widely in nature and its oxidative dissolution might lead groundwater to become acidic, which was harmful to the environment and indeed to artificial building materials.The reactive transport model was a useful tool to predict the extent of such pollution. However, the chemical species were coupled together in the form of a reaction term, which might lead the equations to be nonlinear and thus difficult to solve. A decoupling approach was presented: linear algebraic manipulations of the stoichiometric coefficients of the chemical reactions for the purpose of reducing the number of equation variables and simplifying the reactive source were used. Then the original and decoupled models were solved separately, by both a direct solver and an iterative solver. By comparing the solution times of two models, it was shown that the decoupling approach could enhance the computational efficiency, especially in situations using denser meshes. Using a direct solver, more solution time was saved than when using an iterative version.


Introduction
Pyrite is a common, naturally occurring mineral.In the open atmosphere pyrite oxidative dissolution occurs under the action of groundwater.On one hand the resulting acid water may cause environmental problems, such as contamination of surface and ground waters directed to urban and agricultural supply [1][2][3].Some toxic elements especially, such as arsenic, are closely associated with pyrite.The kinetic oxidative dissolution of As-bearing pyrite due to dissolved oxygen in the ambient groundwater is an important mechanism for arsenic release in groundwater under both natural conditions and engineering applications [4,5].On the other hand the formation of acid water also has some impacts on artificial building materials because of sulfate attack and acid attack [6,7].All the above lead the management of potentially acid generating waste rock to be very important [8].To study the extent and scope of acidic water pollution, some hydrogeochemical models and transport ones are developed to simulate such a system [9][10][11].
In recent years reactive transport model is widely used to simulate the contaminant transport, water-rock interaction, and other processes in earth science fields [12,13].To improve computational efficiency of the model, Friedly and Rubin [14] present a general, concise formulation (decoupling approach) by means of linear algebraic manipulations of the stoichiometric coefficients of the chemical reactions, which can reduce the number of unknown variables and simplify the reaction source/sink terms.Based on this, De Simoni et al. [15] and Molins and Mayer [16] build up the decoupling matrix according to the equilibrium and kinetic reactions.And Huo et al. [17] extend its applications to heterogeneous media.Some efficiency tests are done by Kräutle and Knabner [18] and Hoffmann et al. [19] to study the resulting improvement.In recent years, the decoupling approach is widely used in both engineering applications and laboratory experiments.Saaltink et al. [20] apply the approach into the modelling of multiphase flow for CO 2 injection and storage in deep saline aquifers.And the approach is also used in the simulation of two-phase multicomponent flow with reactive transport in porous media [21,22].And in identifying geochemical processes using end member mixing analysis, Pelizardi et al. [23] uses the decoupling approach to help in the identification of both end members and such reactions, so as to improve mixing ratio calculations.In laboratory experiments and its simulation, the approach is applied in a laboratory experiment where a sand column saturated with a MgSO 4 solution is subject to evaporation [24].And some programs and models are built up based on the decoupling approach for hydrogeochemical calculations, such as CHEPROO++ [25] and MRWM [26].
In this paper, a reactive transport model of pyrite oxidative dissolution is built up in COMSOL Multiphysics, a finite element software platform for the simulation of physicsbased problems.COMSOL is a multiphysics modelling tool that solves various coupled physical problems based on Finite Element Analysis and Partial Differential Equations.It provides a user-friendly interface for mesh generation, equations configuration, and results visualization.And it is widely applied in earth science field.For example, Shao et al. [27] uses it to couple a dual-permeability model with a soil mechanics model for landslide stability evaluation on a hillslope scale.Azad et al. [28] build up an interface between COMSOL and GEMS, a chemical modelling platform, for the reactive transport modelling in variably saturated porous media, while Nardi et al. [29] and Jara et al. [30] couple two standalone simulation programs, COMSOL and PHREEQC, for the reactive transport modelling.
Although some studies have been done on computational efficiency, they are carried out from different research directions.Hoffmann et al. [19] mainly study the impact from theory view, while Kräutle and Knabner' work [18] is based on a transient model to study computational efficiency in different time steps of two approaches.In this paper, we focused on the number of meshes and different solvers.Based on a brief introduction to the theories and mathematical methods behind the decoupling approach, a reactive transport model of pyrite oxidative dissolution is solved by a traditional method and a decoupling approach separately to compare their computational efficiencies.In both 1D and 2D models, the study area is meshed to different grid refinements in each situation.The original and decoupled models are solved and their solution times are compared.Meanwhile, the 2D models are solved by both a direct and an iterative solver to study the effect on the computational efficiency of the decoupling approach compared to different solvers.It is aimed at providing a more convenient and efficient method of calculation to solve the reactive transport model of pyrite oxidative dissolution.

Mathematical Description
The chemical reactions involved in aqueous species are divided into two kinds, equilibrium reaction and kinetic one.Reaction rates of the former are fast in comparison to transport, so that local chemical equilibrium can be assumed at every point within the system.Kinetic laws are applied to represent the processes of latter one, which is not sufficiently fast enough.So without considering the influence of activity, the mass balance of each species can be written in concise vector notation as follows: where vector m contains the mass of species per unit volume of porous medium, and it can be split into two parts, m  and m  , respectively, related to the constant activity species (such as minerals in solid phase and gases) and to the remaining species.Matrix M is diagonal and its diagonal terms are unity when a given species is mobile and zero otherwise.c contains species concentrations in mol/mass of liquid (m =  ⋅ c for mobile species,  is porosity) and c = (c 1 c 2 )  , while c 1 , c 2 are primary and secondary species: the number of secondary species is equal to the number of reaction equilibrium, and the linear operator  in ( 1) is defined as where q is the water flux and D is dispersion coefficient; S is a matrix containing the stoichiometric coefficients of reactions involving reactants and product(s) and S = (S  /S  ), where S  and S  represent the matrices of equilibrium and kinetic reactions such that S  = (S 1 | S 2 ) due to the primary and secondary species.S 1 is stoichiometric coefficients matrix of primary species and S 2 is stoichiometric coefficients matrix of secondary species.
Vector r contains the reaction rate and is also divided into two parts: r  and r  .
A full rank matrix, U, can be established, orthogonal to S  , which satisfies U ⋅ S   = 0.The component matrix U can be calculated by means of Gauss-Jordan elimination which leads to the following expression [31]: where I   −  is a diagonal matrix of dimension   −   , with all diagonal elements equal to one;   and   are the number of reactions and species.Now a component vector of u is defined as u = U ⋅ c and its number,   , can be calculated by   =   −   .Writing the transport equations in terms of u is helpful because the source/sink term becomes simple.Species concentration c can also be solved with the equilibrium reaction constants.According to Molins et al. [32], four types of reactive transport systems are classified by the types of reactions, as shown in Table 1.
It can be seen that four types of reactive transport systems are classified.The characteristics and calculation of component matrix U of each system are shown as follows.
(1) The first is tank system, in which all reactions take place in equilibrium in the aqueous phase, which means a large aqueous reservoir with residence times long enough for aqueous species to reach equilibrium, and no interaction with other solid or gas phases assumed.The component matrix of this system, U tank , can be calculated by the following equation.
where S  1 is stoichiometric coefficients matrix of equilibrium reactions corresponding to primary species.
(2) The second is canal system, in which all reactions are homogeneous, but some may be slow (kinetic).The component matrix of this system, U Canal , can be calculated by the following equation. where 0 I ] and S  1 is stoichiometric coefficients matrix of kinetic reactions corresponding to primary species.
(3) The third is river system, in which heterogeneous reactions also take place, but they are slow relative to flow.The component matrix of this system, U River , can be calculated by the following equation.
where F is a factor matrix which is multiplied by U Canal to eliminate the immobile kinetic species.More detailed solution steps of F can be seen in Molins et al. [32].
(4) The fourth is aquifer system, where some heterogeneous reactions are fast enough to be considered in equilibrium.Some fixed activity species (e.g., minerals and H2O) can be found among the equilibrium reactions.These species can be eliminated from the equations by reducing the components to be solved.The component matrix of this system, U Aquifer , can be calculated by the following equation.
where E is a factor matrix which is multiplied by U River to eliminate constant activity species and reduce the number of components.More detailed solution steps of E can be seen in Molins et al. [32].

Decoupling Approach in Pyrite Oxidative Dissolution
where there are equilibrium reactions in both ( 7) and ( 8), with reaction rates of  1 and  2 and (9) reflects the process of pyrite oxidative dissolution, which depends on the concentration of H + and O 2 (aq) in solution.According to these, the stoichiometric coefficient matrix of the system S could be written as Since the reactions involved both aqueous and solid phases, it satisfied the aquifer system in Section 2. So the component matrix, U, could be calculated as And a new vector of components, u, was defined as where the vector c in this system comprised eight species: O 2 (aq), H + , OH − , SO 4 2− , Fe 3+ , FeS 2 , O 2 (g), and H 2 O in order.The calculated aqueous components were u 1 , u 2 , u 3 .As (12) shows (1) the component vector u is a linear combination of species, which is readily calculated; (2) the number of unknowns to be solved in the equations is reduced, from eight species to three components.In this system, the number of all species is defined as   , with   = 8 and the number of equilibrium reactions   is 2, while the number of secondary species with fixed activity is defined as  0 , which includes pyrite, H 2 O, and O 2 (g).In this system  0 = 3, so the number of components,   , could be calculated as Meanwhile, the reaction terms of species in the transport equations were as follows: Multiplying by the decoupling matrix, U, this term could be expressed as It could be seen that the reaction term of the original model was more complicated and contained the expressions of the equilibrium reaction rates R 1 and R 2 , both of which were difficult to obtain explicitly which introduced some difficulties when solving the model.However, it was expressed in the form of extremely simple items by means of the decoupling approach.As shown in (15), the reaction terms of components, u 1 and u 2 , were 0, and the one of component u 3 involved only R 3 .Then the transport equations of component u could be solved.Once system components had been evaluated, the original species, c, was obtained from the nonlinear algebraic system of ( 12) and corresponding equilibrium constants of ( 7) and (8).

Verification.
In order to verify the influence of the decoupling approach on the calculation accuracy, firstly a batch reactor system of pyrite oxidation dissolution was taken as an example.In this system the pyrite was completely immersed in deionized water in a stirred vessel, which meant there was no need to consider the transport problem.The simulation results by both original approach and decoupling one were compared.The chemical parameters are shown in Table 2.  sp1 and  sp2 are the equilibrium constant of ( 7) and ( 8), while , , and  are reaction parameters of pyrite oxidative dissolution.Its reaction rate could be calculated by the following equation.
where  solid is solid phase surface area and  water is water volume;  1 and  2 are concentration of O 2 (aq) and H + .The ratio of solid phase surface area to water volume was set as 3 dm −1 and simulated time was 10 days.The variation and error of pH value and Fe 3+ concentration in two models are shown in Figure 1.
It can be seen from Figure 1 that the results of two models are basically the same.Maximum relative error of pH value is −0.04725%, while the one of Fe 3+ is 0.59542%, which means that decoupling approach has little effect on the accuracy of the calculation.

Comparison of Computational
Efficiency.The decoupling approach not only simplified the reaction term but also reduced the number of unknown variables in the transport equations.As a result, the new transport model of each component should now be solvable with improved computational efficiency.Deionized water flows through a single smooth fracture of pyrite can be simplified to either a 1D or 2D parallel plate model with model parameters as shown in Table 3.
Initially the fracture was deemed to have been full of deionized water, with a constant flow velocity V through the fracture.Without considering the change in the aperture size caused by dissolution, the distribution of aqueous species reached dynamic equilibrium, which can be regarded as a steady state.The two models were thus simulated.One of the two models involved transport of species c and the original model was designated: the other involved transport of decoupled component u and the decoupled model was designated.
Then the two models were separately established in COMSOL 3.5a, a software platform for the simulation of physics-based problems.The central processing unit (CPU) of the computer was an Intel Core Duo P8400 with a clock speed of 2.26 GHz and the motherboard had 3 GB of random access memory.There were two main categories of solver in the software: direct and iterative.The former included UMFPACK, SPOOLES, PARDISO, and TAUCS Cholesky, which solved a linear system by Gaussian elimination.The iterative solvers, GMRES, FGMRES, conjugate gradients, BiCGStab, and geometric multigrid, were more memoryefficient to deal with models with many degrees of freedom.
When the model was solved in COMSOL, the mesh generator partitioned the study domains into mesh elements: the number of elements depended on the maximum element size when they were uniformly subdivided.Then the 1D and 2D models (original and decoupled types) were solved separately.First the direct solver, UMFPACK, was chosen and the model was solved three times in each case.To solve the nonlinear equations in both original model and decoupled one, Damped Newton Method (DNM) was adopted.Relative tolerance was set as 1.0 × 10 −6 and maximum number of iterations was set as 25.The average solution times are shown in Tables 4 and 5.
It can be seen from Tables 4 and 5 that: (1) Solution by direct solver, UMFPACK, costs much more time in the original model than when adopting a decoupling approach in both 1D and 2D models.
(2) The solution time in both 1D models increases with the number of elements, but it saves more time by using a decoupling approach when the number of elements becomes large.When there are only 25 elements in the model, it costs 0.234 s and 0.094 s to solve each model.With the increase in the number of elements, the solution times reach 1.079 s and 0.297 s for 500 elements (some 4.61 and 3.16 times the requirements at 25 elements).
(3) As in 1D, the computing time in 2D also increases with the number of elements in both of the two models and it saves more time when using a decoupling approach for large numbers of elements.At 7,800 elements, the original model cannot be solved due to an out of memory error during LU factorisation.However, by using a decoupling approach it only costs 21.033 s.Compared to the 1D model, has a better   Then the iterative solver GMRES was chosen to solve the 2D model.The solution set of nonlinear equations was the same as the direct solver, UMFPACK.Three replicates were run and the average solution times are shown in Table 6.
The following can be seen from Table 6.
(1) Like the results in Table 3, the decoupling approach also enhances the computational efficiency when solving the two models by use of the iterative solver.The decoupled solution time is only 17.92% to 52.34% of that needed for the original model and the solution time increases with the number of elements no matter whether in the original or decoupled model.
(2) Unlike the situation in Table 3,  2 / 1 does not reduce with increased numbers of elements: at 5,036 and 7,800 elements,  2 / 1 is only 51.54% and 52.34% of that required originally.It shows that the decoupling approach does not have as significant an effect as expected when used iteratively on a large model.
(3) Solving the original model with a direct solver costs much more time than when using an iterative version.However, when dealing with a decoupled model, the direct solver is faster and  2 / 1 (Table 3) ranges from 4.78% to 27.71%, while it is 17.92% to 52.34% in Table 4.This means that solving the decoupled model by using an iterative solver does not save as much time as the direct solver does.According to Tables 5 and 6, solution times for original and decoupled models with direct and iterative solvers are shown in Figure 2.
Figure 2 shows that solving an original model takes more time than a decoupled one, no matter whether by direct solver or iterative solver.In general sorted by time taken: the decoupled model by direct solver < decoupled model by iterative solver < original model by iterative solver < original model by direct solver.

Conclusions
The present work described the basic theory and mathematical methods of the decoupling approach and then took pyrite oxidative dissolution as an example.Based on the analysis of its chemical reaction system, the decoupling matrix U was calculated.When multiplied by U, the concentration vector c was converted to the component vector u, which had fewer variables and simpler reaction terms.Then the original and decoupled models were established in COMSOL Multiphysics 3.5a.Then the study domain was meshed at different degrees of refinement.In each case it was solved by direct and iterative solvers.The results show the following.
(1) Decoupling enhances the computational efficiency in both 1D and 2D models while saving more time for 2D models than 1D models.

Figure 1 :
Figure 1: Variation and error of pH value and Fe 3+ concentration.

Figure 2 :
Figure 2: Comparison of solution times.

Table 1 :
Types of reactive transport system.
assumed to be formed in deionized water, the main reactions occurring in the open system are as follows.
3.1.TheChemical System and Its Decoupling Matrix.It is important to build up a chemical reaction system for pyrite oxidative dissolution reactive transport.When the initial solution is

Table 2 :
Chemical parameters of pyrite oxidation dissolution.

Table 3 :
Parameters for the 1D and 2D models.

Table 4 :
Comparison of solution time in 1D model using direct solver, UMFPACK.

Table 5 :
Comparison of solution time in 2D model using direct solver, UMFPACK.

Table 6 :
Comparison of solution time in 2D model using iterative solver, GMRES.