Thermodynamic equilibrium solutions through a modified Newton Raphson method

Get Complete Project Material File(s) Now! »

Numerical experiments

A first simplified numerical example is proposed in order to highlight the effects of matrix conditioning and round-off errors on chemical equilibria solutions. The simplified chemical system is composed of three species and two components and is detailed in Table 3.1. The system is solved by the Newton Raphson method (NR), NR with Row-Column scaling (NR-RC), the positive continuous fraction (PCF) method, and PCF coupled with NR. For the initial conditions that are defined in Table 3.1, the associated matrix for the NR method has a condition number of 1.4 103 and the scaled matrix (NR-RC) has a condition number of 9.3 102. These condition numbers are defined in equation (17) by using the infinity norm. As shown by equation (16), the accuracy of the solution depends on the condition number and the round-off errors. Figure 3.1 shows the effects of the round-off errors on the path from the initial guesses to the solution of the above system.
Figure 3.1 – Effect of the round-off error 10-d on the iterative solution. The picture shows the path throughout the solution of four variants of the same simplified problem. Figures (a), (b), (c) and (d) show the results for different values of d, the exponent of the round-off error (d=3, 4, 5 and 6 for, respectively). Since the problem has only two unknowns (components 1 and 2) the path from the common initial guesses ( 1 =6.0 and 2 =3.0) to the solution ( 1 =0.3 and 2 = -2.9) is easily
represented on a 2D graph. Within a single variant of the problem, different paths occur with different algorithms. When round-off error is higher (a) the implementation of scaling makes the difference between reaching the solution or not.
This example shows that the NR method is sensitive to the computational accuracy. In the worst case, the method does not converge (Figure 3.1a) or converges after some amount of iterations, even if the system is not accurately solved during the first iterations (Figure 3.1b and c). Once an acceptable accuracy is defined, NR and NR+RC lead to very similar solutions (Figure 3.1d). NR+RC is much less sensitive to the accuracy because the matrix is better conditioned, even if the contrast in the condition numbers is not very high for this example. The path to the solution does not change for the studied round-off errors. PCF is a first-order method that does not require a system to be solved. The method is provided here to illustrate its advantages and drawbacks; specifically, the method is robust far from the solution but inefficient in the neighborhood of the solution (the computation is stopped after  iterations for this example). Finally, the association of PCF and NR provides an interesting alternative. The algorithm starts with a prescribed number of PCFs (3 in this example) and switches to NR to reach the solution. This association is not sensitive to the round-off errors in this example. Finally, NR may reach the solution in fewer iterations than NR+RC even if the system matrix is not solved properly (Figure 3.1c). The wrong direction of descent for the first iteration is more efficient than the correct one. This process may occur a very few times and is not reliable.
Real or realistic experiments involve a larger number of species and components than those involved in the previous experiment. A visual representation is therefore impossible, but the issues remain the same as before. The efficiencies of the different algorithms on more complex experiments are studied with 7 test cases of increasing complexity. Each test case is solved through the technique clarified in the previous sections: the concentrations of all the chemical species Ci in the system are computed based on a set of fixed components Xj, assuming their total amount T j ¬ in the system is known. The majority of the test cases are available in the literature and considered to be numerically challenging in reason of large range of stoichiometric coefficients and equilibrium constants .The number of species and components for each test case are summarized in Table 3.2.
The first and simplest test case is Gallic Acid, a system proposed by Brassard and Bodurtha (2000) as an example of the onset of problems in numerical methods. The system was originally studied in relation to Al(III) speciation in natural waters (Öhman 1983). This first test case is characterized by the presence of 17 chemical species that can be described through the combination 3 components. Since no solid phase is taken into consideration here, the Jacobian matrix has size 3×3, the smallest of the whole set of test cases. Also the range of variation of equilibrium constants is the smallest of those examined. Increasing complexity is found in MoMaS Easy test case, a synthetic benchmark designed to evaluate the performances of computational codes and published in a special issue of Computational Geosciences (Carrayrou et al. 2009). This test case is characterized by the presence of only 12 species and the number of base components required to describe the system is 5 (2 more than those needed for Gallic acid). At the same time the interval of variation of equilibrium constants is higher (47 against 35 orders of magnitude) than in the previous case. Complexity increases significantly approaching Pyrite31 test cases. This example describes the environment for the potential precipitation of Pyrite ( FeS2 ) . In a first variant of the test case precipitation of a solid phase is denied and the system is composed of 40 dissolved species described through 4 components. The number of components is limited but the differences between equilibrium constants are huge: the lowest is -520 and the highest 19. A second variant of the test case, Pyrite Mineral, is examined. This example is a copy of the previous enriched with the test for the formation of 3 minerals ( Fe,FeSO4,FeS2 ). The size of the Jacobian matrix becomes 7×7. Continuing in the presentation of test cases, it becomes harder to precisely assess the order of complexity. MoMaS Hard test case comes from the same set of synthetic examples of MoMaS Easy (Carrayrou et al. 2009). It s characterized by the presence of 12 chemical species and described through 6 components. Two minerals are tested for precipitation, making the Jacobian matrix 8×8 while the range of equilibrium constants remains the same of MoMaS Easy. Two test cases involving iron and chrome are also studied. Fe Cr test case is the simplest and describes a chemical system of 40 species and 7 components with no test for precipitates. In this case the number of components is higher than in the Pyrite test cases but the range of variation of equilibrium constants is considerably smaller. Fe Cr Min is a variant of the previous case in which 3 minerals are tested for precipitation, making the Jacobian matrix size 10×10. The whole sets of equilibrium constants and a complete description of the chemical species and components are presented in Annex II.
Because the motivation of this work is to evaluate the behavior of the Newton Raphson method and associated algorithms, we only compute the first solution of the chemical system, i.e., if precipitation occurs and no precipitation was assumed, the computation is not repeated. Therefore, we also assumed that the activity coefficients were constant and equal to one. Because the efficiency of the Newton Raphson method is very sensitive to the initial guesses, we searched the solutions of chemical equilibria for 30000 initial component concentrations X j , except for Gallic Acid (10000), which appeared to be the easiest test case. We assumed that 30000 (10000 for Gallic Acid) simulations would adequately represent the behaviors of the convergence rates for each example because the percentage of failure (actually the mean number of iterations that is required to converge) remains unchanged after 25000 runs at most.
Usually, initial guesses for the Newton Raphson method are chosen with great care (i.e., from the concentrations at the previous time step in transient computations). However, initial guesses are not always known, especially for the first time step or for sharp fronts, where the concentrations show abrupt changes from cell to cell in the reactive transport code. Because the aim of this work is to test the robustness of the method no matter the initial conditions, these conditions are chosen randomly in a reasonable interval, as described and explained in Annex III. The robustness of the procedure is also enhanced by prescribing boundaries to the NR increments (Annex III).
We remark that the 30000 solutions are intended for a coupling test case/method of resolution. This means that each test case has been solved 30000 times with each method of resolution. To sum up, different resolutions methods are: regular Newton Raphson method, Newton Raphson method implemented with the four scaling techniques and with the matrix equilibration procedure, regular Newton Raphson coupled with Positive Continuous fractions (implemented at different degrees of ill-conditioning).
The Newton Raphson iterations end when the convergence criterion is satisfied, i.e., when the highest residual is lower than a given threshold ( tol 10 12 ). The Newton iterations are also stopped when the maximum number of iterations exceeds 2000. When possible (with symmetric matrices), the evaluation of the condition number is performed with two approaches: (i) the absolute value of the ratio between the highest and lowest eigenvalues of the Jacobian matrix (3.18) and (ii) the product of the norm one of matrix J and its inverse J 1 (3.17). The subroutine DE3LRG in the IMSL library is used to evaluate the eigenvalues. The linear systems (3.13) and/or (3.27) within the Newton Raphson procedure are solved through LU factorization with quadruple precision. The matrix equilibration algorithm performs scaling on the Jacobian matrix until the infinity norm of each row and column equals one. We set a maximum number of iterations nmaxR 5 to avoid excessive slowdowns in the computation. The algorithm starts with positive continuous fractions when implemented and depending on the condition number, i.e., the algorithm is activated if the condition number is higher than a fixed threshold. Only ten iterations of this method are performed because the aim is to use the algorithm as a type of preconditioner and not to reach the solution (also in light of the results in Figure 3.1).

Numerical simulations: discussion

The purpose of scaling is to reduce the condition number of a linear system. We evaluated the condition number of the Jacobian matrix of the nonlinear system before and after scaling. We presented the results for problems with condition numbers at different orders of magnitude and for a round-off error of 10-32.
When the condition number was on the order of 1060 or lower (as in Figure 2a), both methods of computing the condition number (as the product of the norm one of the Jacobian matrix and its inverse (equation (3.17)) and as the ratio of the eigenvalues (equation (3.18)) were effective and provided the same results.
Figure 3.2 – Variation in condition numbers throughout the minimization for the Pyrite test case with favorable initial guesses (a) and poor initial guesses (b) of the component concentrations.
Under this circumstance, the condition number was generally reduced by RC scaling, even if the
amplitude of this reduction varied widely. When the Jacobian matrix had a condition number on the
order of 10200 or more (Figure 3.2b), the estimation of the condition number was no longer reliable: the inversion of the Jacobian matrix was heavily imprecise, and the DE3LRG subroutine from IMSL
failed to estimate the eigenvalues for condition numbers that were greater than 1050. RC scaling was effective except for condition numbers that were approximately 10120-10150, which shows that scaling may be useful but not universally so. Figure 3.2 shows that the condition number decreased with the distance from the solution. When a large number of iterations were required to converge, the initial values of the condition number were extreme. When a modest number of iterations were necessary to reach the solution, the condition number had smaller values and showed more regular behavior (Machat and Carrayrou 2016).
The results are presented through the relative number of obtained solutions as a function of the number of iterations for each test case (Figures 3.3 to 3.9). When the solution is obtained for all the initial conditions (Ntot=30000 for all test cases except Gallic Acid), the relative number equals 1. The distributions of the initial condition numbers are also provided in Figures 3.3 to 3.9. Table 3.3 provides the number of iterations that are required to solve 50%, 70% and 90% of the problems for each test case and each algorithm. The number of failures (solution not reached after 2000 iterations) is listed in Table 3.4. The results vary widely between test cases, particularly for the implementation of scaling techniques.
In the case of Gallic Acid speciation (Figure 3.3), no differences existed among the curves that represented the standard Newton Raphson method, Newton Raphson with scaling techniques or PCF, which were activated when the condition number was greater than 1020, because most of the condition numbers were below 1020 (Figure 3.3a). Because the condition numbers were smaller than the threshold (1032) that complicates the computation of reliable digits in terms of the round-off error (10-32), the system was solved accurately with or without scaling. When PCFs were performed by default, the number of iterations that were required to solve 100% of the problems dropped from 250 to approximately 30. The results that correspond to the activation of PCFs when the condition number was greater than 1010 provide an intermediate result because some computations were run without PCFs (condition number less than 1010).
The distribution of the initial condition numbers reached 1090 for the MoMaS Easy test case (Figure 4a). Newton Raphson when implemented with scaling performed better than standard Newton Raphson except for DI (Figure 3.4b). 90% of the problems were solved within 140 iterations with the matrix equilibrium technique where other scaling techniques required more than 200 (Table 3.3). The solution of NR coupled with PCFs constantly outperformed the integration of scaling (Tables 3.3 and 3.4). This test case also demonstrated the main shortcoming of the PCF method: as shown in Figure 3.4, PCFs can slow down the solution process in the neighborhood of the solution. When PCFs were applied by default (i.e., regardless if they were needed), almost no solutions were obtained within 20 iterations and only a few within 40. This result confirms that the NR method needs no strengthening when the initial conditions are favorable.
The difficulty in finding a solution increased for the Pyrite test case, which exhibited numerous condition numbers that were greater than 10200 (Figure 3.5a). Approximately 85% of the problems were solved by NR or NR with scaling with the same efficiency, except the matrix equilibration technique which was able to solve 90% of the problems in about 1800 iterations (Table 3.3) and solved 95% of the problems (Table 3.4). The increase in the number of iterations (from approximately 200 to 600) did not change the probability of success in solving the problem. This probability slightly increased after 1000 iterations, but none of the NR methods could solve all the problems. Scaling by MEq was the most efficient scaling technique for this example. The PCFs appeared to be very efficient, solving all the problems with about 50 iterations.
Most of the condition numbers of the standard NR ranged from 1040 to 10160 for the Fe Cr test case (Figure 3.6a). Scaling did not improve the convergence probability (Figure 3.6b and Table 3.3) for that range of condition numbers. The scaling techniques improved the system s matrix properties but not enough to compute an accurate solution of the system. Therefore, scaling was inefficient in this case. Moreover, no significant differences existed in the activation of PCFs after different thresholds of condition numbers, which means that the large majority of the problems had a condition number that was greater than 1040. Pyrite Mineral, exhibited high contrasts in the condition numbers, with 1/3 of the system matrices having a condition number that was greater than 10300 (Figure 3.7a). At this level of complexity, scaling techniques were not efficient enough to reduce the high condition numbers. Condition numbers that were greater than 10300 were indeed reduced but still remained enormous (10260-10280) (Figure 3.7a). Therefore, NR with scaling was as efficient as NR without scaling at solving this system of equations (Figure 3.7b). Again, the PCFs appeared to be very efficient, and 90% of the problems were solved in fewer than 50 iterations (Table 3.3). The MoMaS Hard test case contemplated the presence of precipitates, and its condition numbers ranged from 1010 to 10200 (Figure 8a). Therefore, some problems were solved with a small number of iterations (less than 50 for more than 50% of the initial conditions Table 3.3), except for the MEq. For more difficult problems, the implementation of scaling techniques (except MEq) improved the situation slightly. MEq reduced the condition numbers significantly (Figure 3.8a) and the number of failures more than other scaling techniques but made the process of convergence slower (Figure 3.8b). The implementation of PCFs coupled with the standard Newton Raphson method improved in a more noticeable way both the robustness and speed of the convergence. Again, the standard PCFs required more iterations than the adapted PCFs to solve the easiest problems.

Table of contents :

Chapter 1 – General introduction
1.1 Reactive transport
1.2 Governing equations
1.3 Coupling techniques: operator splitting and global approach
1.4 A universe of reactive transport codes
1.5 LHyGeS numerical models: a garden worth gardening
1.6 Structure of the work
Chapter 2 Implementing isotopes
2.1 Theoretical background
2.1.1 What is an isotope?
2.1.2 Notation
2.1.3 Isotopic abundance and its variations
2.1.5 Why do we care about isotopes?
2.2 Modeling isotopes
2.2.1 Modeling stable isotopes equilibrium fractionation
2.2.2 Modeling stable isotopes kinetic fractionation
2.2.3 Conclusions about modeling isotopes
Chapter 3 Thermodynamic equilibrium
3.1 Thermodynamic equilibrium solutions through a modified Newton Raphson method
3.1.1 Abstract
3.1.2 Introduction
3.1.3 Thermodynamic equilibrium: governing equations
3.1.4 Newton Raphson algorithm
3.1.5 Condition of the linear system
3.1.6 Working on a logarithmic base
3.1.7 Preconditioning
3.1.8 Scaling procedures in this work
3.1.9 Positive continuous fraction method
3.1.10 Numerical experiments
3.1.11 Numerical simulations: discussion
3.1.12 Conclusions about the strategies to improve Newton Raphson method
3.2 Thermodynamic capabilities of the code
3.2.1 Modeling surface complexation
3.2.2 Modeling ion exchange
Chapter 4 Mixed equilibrium and kinetics
4.1 Theoretical background and generic formulation
4.1.1 Generic formulation I
4.1.2 Generic formulation II
4.1.3 Systems of equations
4.2 Solving the systems of equations
4.2.1 Implicit and explicit, one-step or multistep methods of integration
4.2.1.1 Implicit and explicit methods
4.2.1.2 One-step and multi-step methods
4.2.1.3 Variable stepsize
4.2.2 An implemented explicit method: Richardson extrapolation of QSSA method
4.2.3 An implemented implicit method: BDF in DASPK
4.2.4 Solving systems with DASPK
4.2.4.1 Residual computation for DASPK 1
4.2.4.3 Residual computation for DASPK 2
4.2.4.3 Residual computation for DASPK 3
4.3 Numerical simulations
4.3.1 TST model, verification of results with PHREEQC and KINDIS
4.3.1.1 Description of the problem
4.3.1.2 Numerical simulations: results
4.3.2 Chilakapati test case: verification with publication
4.3.2.1 Description of the problem
4.3.2.2 Numerical simulations: results
4.3.2.3 Numerical simulations: effect of convergence criteria
4.4 Conclusions about mixed equilibrium and kinetics
Chapter 5  Solid solutions
5.1 Introduction and theoretical background
5.1.1 The interest in solid solutions
5.1.2 Theoretical background: Thermodynamics of solid solutions
5.1.3 Modeling solid solutions and their interaction with the aqueous phase
5.1.3.1 Equilibrium models for solid solutions
5.1.3.2 Kinetic models for solid solutions
5.1.3.3 Exploiting solid solutions concept for stable kinetic isotope fractionation
5.2 Numerical simulations of solid solutions
5.2.1 Modeling solid solutions at thermodynamic equilibrium
5.2.1.1 Verification with PHREEQC
5.2.1.2 Fe-Cr redox reaction, a reactive transport example
5.3 Conclusions about solid solutions
Chapter 6  Building SpeCTr, a reactive transport code
6.1 Coupling flow, transport and reaction
6.1.1 Governing equations
6.1.2 Global approach
6.1.3 Operator splitting
6.1.4 War of the approaches
6.1.5 Multicomponent reactive transport
6.2 TRACES
6.2.1 Code capabilities
6.2.2 Numerical schemes
6.2.3 Coupling with reaction module
6.2.4 SpeCTr
6.3 Validation: coupling and implementation of isotopes
6.3.1 Interest of validation
6.3.2 Presentation of the problem
6.3.2.1 Spatial discretization, flow characteristics and ground properties
6.3.2.2 Boundary and initial conditions, transport parameters
6.3.2.3 Reaction network
6.3.2.4 Considerations about Courant number
6.3.3 Results of numerical simulations: SpeCTr
6.3.3.1 Results: 􀀁t(CFL), 􀀂L = 0 m, no Cr fractionation
6.3.3.2 Results: 􀀁t(CFL=1), 􀀂L = 0 m, Cr fractionation
6.3.3.3 Results: reduced time step, 􀀂L = 0 m, Cr fractionation
6.3.3.2 Results: 􀀁t (CFL=1), 􀀂L = 1.0 m, Cr fractionation
6.3.3.3 Results: reduced time step, 􀀂L = 0.54 m, Cr fractionation
6.4 Conclusions about SpeCTr validation
Chapter 7 Application of SpeCTr: modeling Calcite dissolution & precipitation
7.1 From mixed flow reactor to column experiments and modeling: upscaling of calcite dissolution rate
7.1.1 Abstract
7.1.2 Introduction
7.1.3 Materials and experimental methods
7.1.3.1 Sample preparations
7.1.3.2 Aqueous solution preparations
7.1.3.3 Mixed flow reactor experiments
7.1.3.4 Column experiment
7.1.3.5 Aqueous sample analyses and thermodynamic calculations
7.1.3.6 Determination of calcite dissolution rate
7.1.4 Mathematical modeling of flow and reactive transport for the column experiments
7.1.5 Results and discussion
7.1.5.1 Mixed-flow reactor experiments
7.1.5.1.a Etching and etch pits morphology
7.1.5.1.b R-􀀃 and R-􀀄G relationships as determined from VSI measurements
7.1.5.2 Column experiment
7.1.5.2.a Dissolution rates determined from VSI measurements
7.1.5.2.b Etch pit morphology
7.1.5.2c Comparison of the mean dissolution rates retrieved with VSI to those inferred from pit morphology
7.1.5.3 Modeled dissolution rates using 2D reactive transport simulations of the column experiment
7.1.5.4 Modeled dissolution rates: 1D versus 2D simulations.
7.1.5.5 Simulation of the Calcium breakthrough curve.
7.1.5.6 1D and 2D-reactive transport simulations of the column experiment: overview and perspectives
7.1.5.6a Mineralogical considerations
7.1.5.6.b Calcium breakthrough
7.1.5 Conclusion
7.2 Preparing Calcite dissolution rate modeling
7.2.1 Computation of reactive surface area for TST and SWM models
7.2.2 Time and spatial discretization
7.3 Mixing induced CaCO3 precipitation
7.3.1 Presentation of the test case
7.3.1.1 Spatial discretization, flow characteristics and ground properties
7.3.1.2 Boundary and initial conditions, transport parameters
7.3.1.3 Reaction network
7.3.1.4 Algorithms for porosity changes
7.3.2 Results of numerical simulations: SPeCTr
7.3.2.1 Results: constant porosity 􀀂 equilibrium and kinetic CaCO3 precipitation
7.3.2.2 Results: variable porosity – equilibrium CaCO3 precipitation
7.3.2.3 Results: variable porosity – kinetic CaCO3 precipitation
7.3.2.4 Conclusion about Calcite precipitation and porosity changes
7.4 – 3D Calcite dissolution modeling
7.4.1 Presentation of the problem
7.4.2 Results of 3D simulation
7.4.2 Conclusions about 3D simulation
General conclusion

GET THE COMPLETE PROJECT