Hydrogeochemical Modelling of Variably Saturated Soils

Abstract

This thesis is about the simulation of flow, transport and geochemical reactions in variably saturated soils. Motivated by the problematic disposal and revegetation of residue from aluminium refining, novel modelling tools for general geochemical reactions in variably saturated flow conditions were developed. In addition, a new model for the effect of mineral reactions on unsaturated hydraulic properties was established. The most sophisticated of the presented series of models was applied in a comprehensive simulation of the geochemical evolution of residue from aluminium refining in response to rainfall- and evaporation-induced leaching. Model development commenced with the direct implementation of the one-dimensional moisture-based form of Richards' equation into the geochemical modelling framework PHREEQC for the simulation of biogeochemical reactions in the vadose zone. The implementation gives access to the full suite of reactions in PHREEQC with activity corrected solution speciation, advanced adsorption models and the flexible definition of kinetic mineral reactions. The numerical approach is based on the conceptual treatment of liquid phase flow and solute transport as time- and space-dependent reactions. The accuracy of the applied finite difference scheme profits from the Kirchhoff transformation on the diffusivity term. In case of the van Genuchten hydraulic model, the integration of soil moisture diffusivity led to analytical series solutions that necessitate numerical approximation. The model's novel capabilities for the simulation of unsaturated flow and ion complexation with variable charge surface sites as well as water content modifying mineral reactions, such as dissolution of hydrated minerals were demonstrated in examples. In a second step, the applicability of this simulation tool was extended to soil profiles with heterogeneous hydraulic properties and variably saturated conditions where parts of the domain may be fully saturated. For this, the head-based form of Richards' equation was implemented using a mass balancing, implicit finite difference scheme with Picard iteration. The use of a total variation-diminishing scheme for solute transport yields second order accurate yet oscillation free solutions for the advective transport of sharp concentration fronts. Applications of the scheme were illustrated for (i) infiltration with saturation-dependent cation exchange capacity, (ii) changes in hydraulic properties due to mineral reactions and (iii) transport of mobile organic substances with complexation of heavy metals. Using this extended simulation tool, an advanced model for the effect of mineral reactions on unsaturated hydraulic properties was developed. The selective radius shift model accounts for the fact that mineral reactions in unsaturated conditions only affect the water filled pore size fraction. In addition, the volume change in each pore due to dissolution/precipitation is proportional to the pore volume. Changes in mineral volume according to geochemical speciation are translated to the pore-size distribution considering the actual moisture content. The resulting discontinuous pore-size distribution is related to the water retention curve by the pore-bundle concept. From the updated discontinuous water retention curve, a smooth water retention function is deduced by means of parameter optimization. Applications of the selective radius shift model were presented for kinetic halite dissolution and calcite precipitation as a consequence of cation exchange. In a further development, the simulation of geochemical reactions in variably saturated soils was extended to two and three-dimensional domains by coupling the general finite element solver COMSOL Multiphysics® to IPhreeqc. IPhreeqc is an object-oriented version of PHREEQC that is specifically designed for couplings to flow and solute transport simulators. The resulting model combines COMSOL's capabilities for the computation of variably saturated flow and solute transport in domains of arbitrary complexity with the full suite of geochemical reaction calculations in PHREEQC. Code verification was performed for point source infiltration and degradation of the pesticide Aldicarb via a kinetic reaction chain. An application example of the new software features subsurface irrigation and fertigation, with cation exchange, proton adsorption and pressure- and space dependent root-water uptake. Using the coupled software, the coarse fraction of the leftover from aluminium refining using Bayer's process was investigated as a covering substrate for residue disposal areas and substrate for plant growth. The solid phase mineral assembly in the residue sand was derived from measured acid neutralization curves. Ion adsorption was modelled and parameterized as an exchange process. A column experiment under saturated flow conditions partially validated the model. The parameterized model was used to predict the geochemical evolution of fertigated and gypsum amended residue sand as response to rainfall- and evaporation-induced leaching in the climatic conditions of southwest Western Australia. Scenarios were calculated for a vertical flow lysimeter and a sloping cover layer for a disposal facility with dimensions 3 m × 50 m. Results show poor long-term retention of fertilizers and the stabilization of the pH at around 10 for at least 5 years of leaching. Soil moisture conditions in the cover layer could be ameliorated by appropriate design of the layer thickness, slope and distance between lateral drains.