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

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 roperties 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 n vadose zone modelling.

Hydrogeochemical Modelling of Variably Saturated Soils

Laurin Wissmeier

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.

EPFL2010

Implementation of variably saturated flow into PHREEQC for the simulation of biogeochemical reactions in the vadose zone

David Andrew Barry, Laurin Wissmeier

A software tool for the simulation of one-dimensional unsaturated flow and solute transport together with biogeochemical reactions in the vadose zone was developed by integrating a numerical solution of Richards' equation into the geochemical modelling framework PHREEQC. Unlike existing software, the new simulation tool is fully based on existing capabilities of PHREEQC without source code modifications or coupling to external software packages. Because of the direct integration, the full set of PHREEQC's geochemical reactions with connected databases for reaction constants is immediately accessible. Thus, unsaturated flow and solute transport can be simulated together with complex solution speciation, equilibrium and kinetic mineral reactions, redox reactions, ion exchange reactions and surface adsorption including diffuse double layer calculations. Liquid phase flow is solved as a result of element advection, where the Darcy flux is calculated according to Richards' equation. For accurate representation of advection-dominated transport, a total- variation-diminishing scheme was implemented. Dispersion was simulated with a standard approach; however, PHREEQC's multicomponent transport capabilities can be used to simulate species-dependent diffusion. Since liquid phase saturation is recalculated after every reaction step, biogeochemical processes that modify liquid phase saturation, such as dissolution or precipitation of hydrated minerals are considered a priori. Implementation of the numerical schemes for flow and solute transport have been described, along with examples of the extensive code verification, before illustrating more advanced application examples. These include (i) the simulation of 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 in varying geochemical conditions.

Elsevier2010

Carbon, oxygen and nitrogen dynamics in a soil profile: Model development and application

David Andrew Barry, Jordi Batlle-Aguilar, Alessandro Brovelli

In order to meet demands for crops, pasture and firewood, the rate of land use change from forested to agricultural uses steadily increased over several decades, resulting in an increased release of nutrients towards groundwater and surface water bodies. In parallel, the degradation of riparian zones has diminished their capacity to provide critical ecosystem functions, such as the ability to control and buffer nutrient cycles. In recent years, however, the key environmental importance of natural, healthy ecosystems has been progressively recognized and restoration of degraded lands towards their former natural state has become an area of active research worldwide. Land use changes and restoration practices are known to affect both soil nutrient dynamics and their transport to neighboring areas. To this end, in order to interpret field experiments and elucidate the different mechanisms taking place, numerical tools are beneficial. Microbial decomposition is the main driving force on biogeochemical transformations of soil organic matter and soil nutrients. The activity of the soil biota is primarily controlled by water availability and by the pore-solution oxygen concentrations, which ultimately depend to a large extent on meteorological conditions, e.g., precipitation. In this work a model is presented that simulates carbon and nitrogen turnover and transport in a 1D profile under variably-saturated conditions. The model is based on the mechanistic batch model of Porporato et al. (Adv. Water Res., 26: 45-58, 2003), but extends its capabilities to simulate the vertical transport of the mobile components. Furthermore, oxygen dynamics are included such that the pore-water concentration is dependent on microbial degradation rates and soil moisture level. The model was applied to simulate the effect of land use change from forested to agricultural soils on pedo-fauna activity and nutrient distribution and abundance across the vertical profile. External forcing, i.e., precipitation time series, were generated stochastically. Modeling results showed that soil tillage practices, which modify soil structure and thus soil aeration, were responsible for higher decomposition and mineralization rates in agricultural soils, thus for the lower soil carbon and nitrogen concentrations. Furthermore, higher plant uptake rates and leaching of agricultural soils contributed also to lower nitrate concentrations in such soils.

2009