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

Authors: Jixiang Huo; Fuheng Ma; Hanzhou Song
Journal: Geofluids
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 CO2 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 MgSO4 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 physics-based 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.


## 2. 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:(1)∂m∂t=MLc+ST·r=MLc+SeTre+SkTrk,where vector m contains the mass of species per unit volume of porous medium, and it can be split into two parts, mc and ma, 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=(c1c2)T, while c1,  c2 are primary and secondary species: the number of secondary species is equal to the number of reaction equilibrium, and the linear operator L in (1) is defined as L(c)=-∇·(q·c)+∇·(θ·D∇c), where q is the water flux andD is dispersion coefficient; S is a matrix containing the stoichiometric coefficients of reactions involving reactants and product(s) and S=Se/Sk, where Se and Sk represent the matrices of equilibrium and kinetic reactions such that Se=S1∣S2 due to the primary and secondary species. S1 is stoichiometric coefficients matrix of primary species and S2 is stoichiometric coefficients matrix of secondary species. Vector r contains the reaction rate and is also divided into two parts: re and rk.A full rank matrix,U, can be established, orthogonal to Se, which satisfies U·SeT=0. The component matrix U can be calculated by means of Gauss-Jordan elimination which leads to the following expression [31]:(2)U=INs-Ne∣-S1T·S2T-1,where INs-Ne is a diagonal matrix of dimension Ns-Ne, with all diagonal elements equal to one; Ne and Ns are the number of reactions and species. Now a component vector of u is defined as u=U·c and its number, Nu, can be calculated by Nu=Ns-Ne. 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.Table 1
Types of reactive transport system.
Paradigm
Homogeneous
Heterogeneous
Fast
Slow
Fast
Slow
(Equilibrium)
(Kinetic)
(Equilibrium)
(Kinetic)
Tank
√
Canal
√
√
River
√
√
√
Aquifer
√
√
√
√It can be seen that four types of reactive transport systems are classified. The characteristics and calculation of component matrixU 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,Utank, can be calculated by the following equation.(3)Utank=Ue=ISe1t,where Se1t 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,UCanal, can be calculated by the following equation.(4)UCanal=Uk·Ue,where Uk=ISk1t0I and Sk1t 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,URiver, can be calculated by the following equation.(5)Uriver=F∙Ucanal,whereF is a factor matrix which is multiplied by UCanal to eliminate the immobile kinetic species. More detailed solution steps ofF 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,UAquifer, can be calculated by the following equation.(6)UAquifer=E·URiver,whereE is a factor matrix which is multiplied by URiver 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].


## 3. Decoupling Approach in Pyrite Oxidative Dissolution

### 3.1. The Chemical 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 assumed to be formed in deionized water, the main reactions occurring in the open system are as follows.(7)O2g⟺O2aq(8)H2O⟺H++OH-(9)FeS2+3.75O2aq+0.5H2O⟹Fe3++2SO42-+H+where there are equilibrium reactions in both (7) and (8), with reaction rates of R1 and R2 and (9) reflects the process of pyrite oxidative dissolution, which depends on the concentration of H+ and O2(aq) in solution. According to these, the stoichiometric coefficient matrix of the systemS could be written as(10)S=O2aqH+OH-SO42-Fe3+FeS2O2gH2O100000-100110000-1-3.751021-10-0.5.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(11)U=0-11010000001-200000001000.And a new vector of components, u, was defined as(12)u=U·c=-c2+c3+c5c4-2·c5c5=u1u2u3,where the vector c in this system comprised eight species: O2(aq), H+, OH-, SO42-, Fe3+, FeS2, O2(g), and H2O in order. The calculated aqueous components wereu1,u2,u3. 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 Ns, with Ns=8 and the number of equilibrium reactions Ne is 2, while the number of secondary species with fixed activity is defined as N0, which includes pyrite, H2O, and O2(g). In this system N0=3, so the number of components, Nc, could be calculated as(13)Nc=Ns-Ne-N0=8-2-3=3.Meanwhile, the reaction terms of species in the transport equations were as follows:(14)ST·R=R1-3.75·R3,R2+R3,R2,2·R3,R3,-R3,-R1,-R2-0.5·R3T.Multiplying by the decoupling matrix,U, this term could be expressed as(15)U·ST·R=000000001·R1R2R3=00R3.It could be seen that the reaction term of the original model was more complicated and contained the expressions of the equilibrium reaction ratesR1 andR2, 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,u1 andu2, were 0, and the one of componentu3 involved onlyR3. 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).


### 3.2. 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 Table2.Table 2
Chemical parameters of pyrite oxidation dissolution.
K
s
p
1
K
s
p
2
k
m
n
10−2.8983
10−14
10−10.19
0.5
−0.11K
s
p
1 and Ksp2 are the equilibrium constant of (7) and (8), while k, m, and n are reaction parameters of pyrite oxidative dissolution. Its reaction rate could be calculated by the following equation.(16)R3=AsolidVwater·k·c1m·c2n,where Asolid is solid phase surface area and Vwater is water volume; c1 and c2 are concentration of O2(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 Fe3+ concentration in two models are shown in Figure 1.Figure 1
Variation and error of pH value and Fe3+ concentration.It can be seen from Figure1 that the results of two models are basically the same. Maximum relative error of pH value is −0.04725%, while the one of Fe3+ is 0.59542%, which means that decoupling approach has little effect on the accuracy of the calculation.


### 3.3. 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 Table3.Table 3
Parameters for the 1D and 2D models.
L/m
W/m
v/m·s−1
α
L/m
α
T/m
D
m/m2·s−1
1D
5
—
0.001
0.05
—
1.34
e
-
9
2D
5
1
0.001
0.05
0.02
1.34
e
-
9Initially the fracture was deemed to have been full of deionized water, with a constant flow velocityv 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 memory-efficient 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.Table 4
Comparison of solution time in 1D model using direct solver, UMFPACK.
Maximum element size
Number of elements
Original solution time,t1/s
Decoupled solution time,t2/s
t
2
/
t
1 %
0.01
500
1.079
0.297
27.53
0.02
250
0.625
0.203
32.48
0.04
125
0.406
0.141
34.73
0.05
100
0.359
0.125
34.82
0.08
63
0.313
0.109
34.82
0.10
50
0.282
0.109
38.65
0.20
25
0.234
0.094
35.47Table 5
Comparison of solution time in 2D model using direct solver, UMFPACK.
Maximum element size
Number of elements
Original solution time,t1/s
Decoupled solution time,t2/s
t
2
/
t
1 %
0.04
7800
-
21.033
-
0.05
5036
246.141
11.767
4.78
0.08
2038
43.644
4.141
9.49
0.10
1310
18.095
2.344
12.95
0.12
944
11.032
1.719
15.58
0.15
616
6.062
1.141
18.82
0.20
308
2.672
0.626
20.08
0.25
206
1.750
0.485
23.43
0.30
148
1.344
0.406
27.71
- represents the absence of a result because of an out of memory error during LU factorisation.It can be seen from Tables4 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, it has a better computational efficiency costing only 4.78% to 27.71% of the original.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 Table6.Table 6
Comparison of solution time in 2D model using iterative solver, GMRES.
Maximum element size
Number of elements
Original solution time,t1/s
Decoupled solution time,t2/s
t
2
/
t
1 %
0.04
7800
89.644
46.923
52.34
0.05
5036
50.955
26.313
51.64
0.08
2038
18.251
8.235
45.12
0.10
1310
11.016
5.062
45.95
0.12
944
7.328
3.454
47.13
0.15
616
4.531
0.812
17.92
0.20
308
2.172
0.484
22.28
0.25
206
1.547
0.406
26.24
0.30
148
1.219
0.360
29.53The following can be seen from Table6.(1) Like the results in Table3, 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 Table3, t2/t1 does not reduce with increased numbers of elements: at 5,036 and 7,800 elements, t2/t1 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 andt2/t1 (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 Tables5 and 6, solution times for original and decoupled models with direct and iterative solvers are shown in Figure 2.Figure 2
Comparison of solution times.Figure2 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.


## 3.1. The Chemical 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 assumed to be formed in deionized water, the main reactions occurring in the open system are as follows.(7)O2g⟺O2aq(8)H2O⟺H++OH-(9)FeS2+3.75O2aq+0.5H2O⟹Fe3++2SO42-+H+where there are equilibrium reactions in both (7) and (8), with reaction rates of R1 and R2 and (9) reflects the process of pyrite oxidative dissolution, which depends on the concentration of H+ and O2(aq) in solution. According to these, the stoichiometric coefficient matrix of the systemS could be written as(10)S=O2aqH+OH-SO42-Fe3+FeS2O2gH2O100000-100110000-1-3.751021-10-0.5.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(11)U=0-11010000001-200000001000.And a new vector of components, u, was defined as(12)u=U·c=-c2+c3+c5c4-2·c5c5=u1u2u3,where the vector c in this system comprised eight species: O2(aq), H+, OH-, SO42-, Fe3+, FeS2, O2(g), and H2O in order. The calculated aqueous components wereu1,u2,u3. 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 Ns, with Ns=8 and the number of equilibrium reactions Ne is 2, while the number of secondary species with fixed activity is defined as N0, which includes pyrite, H2O, and O2(g). In this system N0=3, so the number of components, Nc, could be calculated as(13)Nc=Ns-Ne-N0=8-2-3=3.Meanwhile, the reaction terms of species in the transport equations were as follows:(14)ST·R=R1-3.75·R3,R2+R3,R2,2·R3,R3,-R3,-R1,-R2-0.5·R3T.Multiplying by the decoupling matrix,U, this term could be expressed as(15)U·ST·R=000000001·R1R2R3=00R3.It could be seen that the reaction term of the original model was more complicated and contained the expressions of the equilibrium reaction ratesR1 andR2, 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,u1 andu2, were 0, and the one of componentu3 involved onlyR3. 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).


## 3.2. 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 Table2.Table 2
Chemical parameters of pyrite oxidation dissolution.
K
s
p
1
K
s
p
2
k
m
n
10−2.8983
10−14
10−10.19
0.5
−0.11K
s
p
1 and Ksp2 are the equilibrium constant of (7) and (8), while k, m, and n are reaction parameters of pyrite oxidative dissolution. Its reaction rate could be calculated by the following equation.(16)R3=AsolidVwater·k·c1m·c2n,where Asolid is solid phase surface area and Vwater is water volume; c1 and c2 are concentration of O2(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 Fe3+ concentration in two models are shown in Figure 1.Figure 1
Variation and error of pH value and Fe3+ concentration.It can be seen from Figure1 that the results of two models are basically the same. Maximum relative error of pH value is −0.04725%, while the one of Fe3+ is 0.59542%, which means that decoupling approach has little effect on the accuracy of the calculation.


## 3.3. 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 Table3.Table 3
Parameters for the 1D and 2D models.
L/m
W/m
v/m·s−1
α
L/m
α
T/m
D
m/m2·s−1
1D
5
—
0.001
0.05
—
1.34
e
-
9
2D
5
1
0.001
0.05
0.02
1.34
e
-
9Initially the fracture was deemed to have been full of deionized water, with a constant flow velocityv 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 memory-efficient 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.Table 4
Comparison of solution time in 1D model using direct solver, UMFPACK.
Maximum element size
Number of elements
Original solution time,t1/s
Decoupled solution time,t2/s
t
2
/
t
1 %
0.01
500
1.079
0.297
27.53
0.02
250
0.625
0.203
32.48
0.04
125
0.406
0.141
34.73
0.05
100
0.359
0.125
34.82
0.08
63
0.313
0.109
34.82
0.10
50
0.282
0.109
38.65
0.20
25
0.234
0.094
35.47Table 5
Comparison of solution time in 2D model using direct solver, UMFPACK.
Maximum element size
Number of elements
Original solution time,t1/s
Decoupled solution time,t2/s
t
2
/
t
1 %
0.04
7800
-
21.033
-
0.05
5036
246.141
11.767
4.78
0.08
2038
43.644
4.141
9.49
0.10
1310
18.095
2.344
12.95
0.12
944
11.032
1.719
15.58
0.15
616
6.062
1.141
18.82
0.20
308
2.672
0.626
20.08
0.25
206
1.750
0.485
23.43
0.30
148
1.344
0.406
27.71
- represents the absence of a result because of an out of memory error during LU factorisation.It can be seen from Tables4 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, it has a better computational efficiency costing only 4.78% to 27.71% of the original.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 Table6.Table 6
Comparison of solution time in 2D model using iterative solver, GMRES.
Maximum element size
Number of elements
Original solution time,t1/s
Decoupled solution time,t2/s
t
2
/
t
1 %
0.04
7800
89.644
46.923
52.34
0.05
5036
50.955
26.313
51.64
0.08
2038
18.251
8.235
45.12
0.10
1310
11.016
5.062
45.95
0.12
944
7.328
3.454
47.13
0.15
616
4.531
0.812
17.92
0.20
308
2.172
0.484
22.28
0.25
206
1.547
0.406
26.24
0.30
148
1.219
0.360
29.53The following can be seen from Table6.(1) Like the results in Table3, 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 Table3, t2/t1 does not reduce with increased numbers of elements: at 5,036 and 7,800 elements, t2/t1 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 andt2/t1 (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 Tables5 and 6, solution times for original and decoupled models with direct and iterative solvers are shown in Figure 2.Figure 2
Comparison of solution times.Figure2 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.


## 4. 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 matrixU was calculated. When multiplied byU, 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.(2) The more mesh grids the domain generates, the more efficiently the decoupled model finds a solution by direct solver, whether in 1D or 2D.(3) Although the iterative solver takes less time than the direct solver for the original 2D model, it is more efficient to use a direct solver in solving a decoupled problem.(4) The solution times in ascending order are the decoupled model solved by direct solver, a decoupled model solved by an iterative solver, the original model solved by an iterative solver, and the original model solved by a direct solver.As a conclusion, the decoupling approach is of assistance when solving reactive transport of pyrite oxidative dissolution problems, especially over a large domain with more mesh elements. Its applicability is thus demonstrated.



---

*Source: 4670103-2017-11-26.xml*