Transcript Slide 1

Simulation Tool for Variably Saturated Flow with
Comprehensive Geochemical Reactions in 2D and 3D
Coupling COMSOL
Laboratoire de Technologie Ecologique
Station 2, CH-1015 Lausanne, Switzerland
[email protected]
http://ecol.epfl.ch; http://people.epfl.ch/laurin.wissmeier
®
Multiphysics
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.
MODEL CAPABILITIES
2D and 3D variably saturated flow and multi-component
solute transport in COMSOL
MODEL VERIFICATION
Example 16.5. from COMSOL model library (Pesticide Transport and Reactions in Soil, modified); Chemistry
modelled by PHREEQC
•
•
•
•
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
Localized infiltration of Aldicarb into initially dry, layered soil.
Kinetic decomposition into 5 daughter products according to first order rate expressions.
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 COMobject until the following reaction step.
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
Simulation domain:
Irrigation trench
Geochemical Reactions:
•
Kinetic Calcite dissolution/precipitation rate r (Plummer et al. 1978):
Symmetry axis
Convenience of COMSOL to define irregular, heterogeneous
flow domains
van Genuchten/Mualem or user-defined hydraulic property
models
User-defined storativity
•
and PHREEQC
•
Fig. 2: Radial symmetric simulation domain
(from COMSOL documentation).
KCalcite =
Water table
a = 10; b = 0.6; θ: Moisture content; T: Temperature °C; [ ]: Ion activities
Cation Exchange:
2×10-3 mol/L soil negatively charged exchange sites
10-8.48;
Fig. 6: 2D-simulation domain with FEM discretisation.
Fig. 3: Reaction pathways of Aldicarb degradation
(from COMSOL documentation).
Scenario 1: Exchange
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
COMSOL_in.m
Fig. 7: Simulation results with cation exchange 20 h after beginning of infiltration. Contour lines: equal moisture content; black arrows: liquid flow velocity field.
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.
Liquid Phase Flow
COMSOL
Phase flow
velocities
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.
Scenario 2: Exchange + Kinetic Calcite Precipitation/Dissolution
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.
Solute Transport
COMSOL_out.m
PHREEQC
PHREEQC_in.m
Reactions
intra-phase speciation,
inter-phase speciation,
kinetic reactions,
sources and sinks
Phase element concentrations,
mass of water
etc.
PHREEQC_out.m
Fig. 1: Program flow with program files.
Numerical time step Δt
Change in element amounts
due to advection
diffusion/dispersion
Fig. 5: Comparison of simulation results from COMSOL the coupling of COMSOL and PHREEQC, where chemical reactions are simulated by PHREEQC.
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.
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 calciumdominated 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.