PHREEQC. COMSOL. Simulation Tool for Variably Saturated Flow with Comprehensive Geochemical Reactions in 2D and 3D Coupling COMSOL Multiphysics ® and PHREEQC. Laboratoire de Technologie Ecologique Station 2, CH-1015 Lausanne, Switzerland
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
[email protected]://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.
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).
Numerical time step Δt
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.
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
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
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.
Fig. 1: Program flow with program files.