Model Verification

1. PHREEQC COMSOL Simulation Tool for Variably Saturated Flow withComprehensive Geochemical Reactions in 2D and 3DCoupling COMSOL Multiphysics® and PHREEQC Laboratoire de Technologie Ecologique Station 2, CH-1015 Lausanne, Switzerland laurin.wissmeier@epfl.chhttp://ecol.epfl.ch; http://people.epfl.ch/laurin.wissmeier L. Wissmeier, D.A. Barry We present a coupling between the multipurpose finite element solver COMSOL Multiphysics® and the geochemical modelling framework PHREEQC for simulations of flow and multispecies solute transport in 2 and 3-dimensional domains in combination with complex intra-phase and inter-phase geochemistry. To demonstrate COMSOL’s capabilities to simulate environmental flow phenomena, we solve Richards’ equation for variably saturated porous media. As a major advantage over COMSOL’s built-in reaction capabilities, the coupling to PHREEQC accounts for ion activities, and implements major state-of-the-art adsorption models such as ion exchange and surface adsorption with diffuse double layer calculations. In addition, PHREEQC provides a framework to integrate user-defined kinetic reactions with possible dependencies on solution speciation (e.g., pH). Extensive compilations of geochemical reactions and their parameterization are accessible through associated databases. The coupling procedure employs a simple sequential split operator to separate the computation of spatial derivatives for flow and transport from reaction calculations. PHREEQC, compiled as a Windows COM-object (D.L. Parkhust, USGS, in development), updates the geochemical status of the liquid and solid phases after changes in the solution composition due to solute transport calculations in COMSOL. Matlab® is conveniently used as the controlling application for the overall simulation process, and specifically to direct output from COMSOL to PHREEQC and vice versa. The PHREEQC COM-version facilitates a ‘close coupling’ approach, where information between the transport and reaction module is rapidly exchanged via the memory stack at runtime as opposed to slow reading and writing of data files. At the end of a reaction step, COMSOL is reinitialized with concentrations of each element in the mobile liquid phase (including oxygen and hydrogen) and the liquid phase saturation. After the following transport step, the element concentrations are converted to total moles of elements using the updated liquid phase saturation and the elemental composition of the solution in PHREEQC is adjusted. Subsequently, solution speciation, equilibration with mineral phases and adsorbing surfaces, and kinetic reactions are carried out. As major advantage of this procedure, only mobile phase elements and their oxidation states must be passed between the flow and reaction module. The status of immobile adsorbents (surfaces and exchange sites) and mineral compositions is stored by the PHREEQC COM-object until the following reaction step. • Model Capabilities • 2D and 3D variably saturated flow and multi-component solute transport in COMSOL • Convenience of COMSOL to define irregular, heterogeneous flow domains • van Genuchten/Mualem or user-defined hydraulic property models • User-defined storativity • Comprehensive geochemical reactions in PHREEQC • Complex solution speciation including activity corrections • Equilibrium mineral dissolution/precipitation • Permanent charge adsorption (ion exchange) • Variable charge surface complexation including diffuse layer calculations • Redox reactions • Gas phase exchange reactions • Kinetic reactions with user-defined rate equations • Direct access to extensive reaction databases distributed with PHREEQC • Software Requirements • COMSOL Multiphysics® (version 3.5 or higher) • Matlab® (version 7.0 or higher) • PHREEQC (COM version) • The coupling files are available for scientific purposes upon request • Model Verification • Example 16.5. from COMSOL model library (Pesticide Transport and Reactions in Soil, modified); Chemistry modelled by PHREEQC • Localized infiltration of Aldicarb into initially dry, layered soil. • Kinetic decomposition into 5 daughter products according to first order rate expressions. • Cause of deviations between COMSOL and PHREEQC+COMSOL • Split-Operator Error • The separate treatment of the equations of flow, transport (PDE’s) and reaction equations (coupled system of ODE’s) leads to a splitting error. The magnitude of splitting error depends on the time increment between flow/transport (COMSOL) and reaction calculations (PHREEQC). The coupling scheme allows for the flexible definition of the coupling time steps with a possible dependency on flow properties (e.g., maximum solute velocity). • Negative Solute Concentrations in COMSOL’s Transport Computation • The numerical treatment of transport in COMSOL leads to oscillatory solute concentrations near sharp concentration fronts. This may lead to negative concentrations locally. Negative concentrations are taken as zero in PHREEQC reaction calculations, which is why a small mass balance error may be introduced. Negative concentrations can be avoided by using logarithmic transformation of concentrations in the finite element scheme. However, this may introduce larger absolute concentration oscillations. • Application Example • Infiltration of Sodium-Carbonate Solution From Irrigation Trench into Calcite Containing Loam • Description: • Infiltration of sodium-carbonate solution (Na2CO3 0.0116 M; pH 9.44) in equilibrium with atmospheric O2 (partial pressure 10-0.7) and CO2 (pp 10-3.5) from irrigation trench (const. hydraulic head: 0 m). • Initially wetted soil (hydraulic head: -1m, van Genuchten/Mualem hydraulic properties of loam). • Initial soil solution (pH 6.97) in equilibrium with Calcite (CaCO3) and root zone CO2 (pp 10-1.5). • Bottom boundary at water table (const. hydraulic head: -1m); advective flux boundary for solutes • Geochemical Reactions: • Kinetic Calcite dissolution/precipitation rate r (Plummer et al. 1978): • KCalcite= 10-8.48; a = 10; b = 0.6; θ: Moisture content; T: Temperature °C; [ ]: Ion activities • Cation Exchange:2×10-3 mol/L soil negatively charged exchange sites • Scenario 1: Exchange Simulation domain: Irrigation trench Symmetry axis Water table COMSOL_in.m Fig. 6: 2D-simulation domain with FEM discretisation. Liquid Phase Flow Phase flow velocities Change in element amounts due to advection diffusion/dispersion Fig. 2: Radial symmetric simulation domain (from COMSOL documentation). Solute Transport Numerical time step Δt COMSOL_out.m PHREEQC_in.m Fig. 4: Simulation results for Aldicarb and concentrations of its toxic daughter products (color), pressure head (white contour lines) and liquid flow velocity field (black arrows) 5 days after beginning of infiltration event. Reactions intra-phase speciation, inter-phase speciation, kinetic reactions, sources and sinks Fig. 7: Simulation results with cation exchange 20 h after beginning of infiltration. Contour lines: equal moisture content; black arrows: liquid flow velocity field. Phase element concentrations, mass of water etc. The infiltrating moisture front has reached the phreatic groundwater surface. Flow velocities are highest directly below the irrigation trench. Calcium concentration shows a distinct snow-plough at the sodium front. Due to continuous desorption from the exchanger, concentrations ahead of the sodium front exceed their original value by a factor of 3. The exchanger composition reflects the replacement of calcium by infiltrating sodium. Due to the preferences of the exchanger for divalent calcium, the exchanger composition shows a smooth transition from calcium dominated to sodium controlled. Scenario 2: Exchange + Kinetic Calcite Precipitation/Dissolution PHREEQC_out.m Fig. 3: Reaction pathways of Aldicarb degradation (from COMSOL documentation). Fig. 5: Comparison of simulation results from COMSOL the coupling of COMSOL and PHREEQC, where chemical reactions are simulated by PHREEQC. Fig. 8: Simulation results with cation exchange and kinetic calcite dissolution/precipitation 20 h after beginning of infiltration. Calcite precipitates as a kinetic phase. Due to elevated carbonate concentrations in the infiltration water, calcite remains relatively stable in the once precipitated. The snow-plough as a result of exchange is limited by kinetic calcite precipitation, and calcium ahead of the sodium front only slightly exceeds its original concentration. The exchanger composition abruptly changes from calcium-dominated to sodium-dominated, since behind the concentration front calcium concentrations are reduced through precipitation with the high carbonate concentrations in the infiltration water. • Concluding Remarks • The coupling of COMSOL and PHREEQC provides novel capabilities for the detailed simulation of geochemical reactions in variably saturated flow conditions on 2 and 3D domains. • Extensions to other environmental flow phenomena (2-phase flow, density dependent reactive flow) are envisaged. • Drawbacks of the coupling to COMSOL are slow computation times for flow and transport after re-initialization of element concentrations and oscillations of element concentrations near sharp concentration front, which may lead to fluctuations of the redox potential and small mass balance errors. Fig. 1: Program flow with program files.