Reactive transport in unsaturated soil: Comprehensive modelling of the dynamic spatial and temporal mass balance of water and chemical components

Abstract

By implementing the moisture-based form of Richards’ equation into the geochemical modelling framework PHREEQC, a generic tool for the simulation of one-dimensional flow and solute transport in the vadose zone undergoing complex geochemical reactions was developed. A second-order, cell-centred, explicit finite difference scheme was employed for the numerical solution of the partial differential equations of flow and transport. In this scheme, the charge-balanced soil solution is treated as an assembly of elements, where changes in water and solute contents result from fluxes of elements across cell boundaries. Therefore, water flow is considered in terms of oxygen and hydrogen transport. The direct implementation into the geochemical framework provides access to the full set of reactions available in PHREEQC, giving capabilities beyond existing software for unsaturated flow and reaction. Possible reactions include complex aqueous speciation, cation exchange, equilibrium phase dissolution and precipitation, formation of solid solutions, redox reactions, gas phase exchange, surface adsorption considering electrostatics and kinetic reactions with user-defined rate equations, among others. Geochemical reactions were coupled to transport processes by non-iterative sequential operator splitting. The scheme is currently limited to cases where changes in physical fluid properties and hydraulic flow characteristics due to geochemical reactions are negligible. Results from extensive code verification with analytical and accurate numerical solutions as well as HYDRUS-1D show the excellent performance of the scheme for a variety of hydraulic models including the Brooks and Corey model and the van Genuchten model. High accuracy was gained by the use of integrated diffusivities in the finite difference formulation. The integration of complex geochemical reactions was verified with HP1 by simulating the infiltration of a hyperalkaline solution into a clay soil involving aqueous speciation, equilibrium and kinetic phase dissolution/precipitation and cation exchange reactions. The novel capability of the scheme to account for the influence of geochemical reactions on water contents was demonstrated by comparing reactive and non-reactive transport under transient flow conditions. The simulation of surface complexation according to the diffuse double layer model with an explicit calculation of the ion composition in the diffuse layer in transient unsaturated flow conditions, as shown, represents an additional advance in vadose zone modelling.