Numerical validation of a Eulerian hydrochemical code using a 1D multisolute mass transport system involving heterogeneous kinetically controlled reactions

Abstract

It is demonstrated that at steady state, the 1D thermo-kinetic hydrochemical Eulerian mass balance equations in pure advective mode are indeed identical to the governing mass balance equations of a single reaction path (or geochemical) code in open system mode. Thus, both calculated reaction paths should be theoretically identical whatever the chemical complexity of the water–rock system (i.e., multicomponent, multireaction zones kinetically and equilibrium-controlled). We propose to use this property to numerically test the thermo-kinetic hydrochemical Eulerian codes and we employ it to verify the algorithm of the 1D finite difference code KIRMAT. Compared to the other methods to perform such numerical tests (i.e., comparisons with analytical, semi-analytical solutions, between two Eulerian hydrochemical codes), the advantage of this new method is the absence of constraints on the chemical complexity of the modelled water–rock systems. Moreover, the same thermo-kinetic databases and geochemical functions can be easily and mechanically used in both calculations, when the numerical reference comes from the Eulerian code with no transport terms ( u and D=0) and modify to be consistent with the definition of the open system mode in geochemical modelling. The ability of KIRMAT to treat multicomponent pure advective transport, subjected to several kinetically equilibrium-controlled dissolution and precipitation reactions, and to track their boundaries has been successfully verified with the property of interest. The required numerical validation of the reference calculations is bypassed in developing the Eulerian code from an already checked single reaction path code. A forward time-upstream weighting scheme (a mixing cell scheme) is used in this study. An appropriate choice of grid spacing allows to calculate within the grid size uncertainty the correct mineral reaction zone boundaries, despite the presence of numerical dispersion. Its correction enables us to improve the convergence and to extend the numerical test to mixed advective–dispersive mass transport. However, the skewness factor involves numerical oscillations that prevent to compute different grid spacing. The use of a different chemically controlled time step constraint in both calculations induces some inconsistencies into the validation tests. This numerical validation method may be applied as well as to check a thermo-kinetic hydrochemical finite element based code, from a 1D heterogeneous systems, and 2D–3D systems provided that they are designed so as to be 1D equivalent. A one-step algorithm and the use of a numerical reference coming from the Eulerian code to be tested ensure the potential success (accuracy) of the numerical validation method.