Coupled THMC Processes in Radionuclide Waste Management
Gour-Tsyh (George) Yeh, CAMRDA
College of Earth Sciences, National Central University
No 300 Jhongda Road, Jhongli District, Taoyuan City 32001, Taiwan

Chia-Hsing (Peter) Tsai
Water Sciences and Management, 253 Skeen Hall, New Mexico State University
1780 East University Ave., Las Cruces, NM 88003

Key words; Multiphase Flow, Thermal Transport, Reactive Transport, Geo-Mechanics

Abstract
Physics-based mathematical-computational models provide an invaluable tool for repository design,
performance/safety assessments, site clean-up, and environmental remediation in nuclear spent fuel and highlevel waste disposal in shallow or deep geological media. The migration of nuclear wastes is controlled bycoupled THMC (Thermal-Hydrology-Mechanics-Chemical) processes. These include multiphase flow orvariably saturated flow, thermal transport, reactive transport, and geo-mechanical deformation. Multiphase orvariably saturated flow in subsurface media is the key driving mechanism for thermal transport, chemicaltransport, and geo-mechanics. Conversely thermal transport, chemical transport, and geo-mechanic deformation
will affect flow directly or indirectly. Temperature changes are induced by fluid injection, associated phasechanges, and chemical reactions. Conversely, temperature change will alter fluid flow, reactive transport, andgeo-mechanical deformation. Faults and fractures will affect fluid pressure, thermal transport, and chemical
migration. Conversely fluid pressure, thermal transport, and chemical migration will induce rock deformationand fault displacement. Reactive transport is controlled by fluid flow, affected by thermal transport, andinfluenced by geo-mechanical deformation. Conversely, reactive transport will have significant feedbacks on
fluid flow, thermal transport, and geo-mechanics. This talk presents the development of a series of computationalmodels that fully or partially couple these processes. Demonstrative applications are presented to illustrate theinterplay of THMC processes and their implications for radionuclide waste managements.

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan1. Introduction
The key issue in groundwater modeling is the ability
to mathematically describe and numerically simulate
coupled or decoupled physical, chemical, andbiological processes, and their interactions with themedia. The media include heterogeneous porous
media and discrete fractures. The processes include
fluid flow, thermal transport, reactive transport andbio-geochemistry, and geo-mechanical deformation.
These are widely known as the THMC (Thermal,
Hydrology or Hydraulics, geo-Mechanics, andreactive bio-geo-Chemistry) processes. Subsurfacehydrology deals with flow in single or multi-fluid
phases, including aqueous, super liquid, NAPL
(non-aqueous phase liquid), gaseous phases, etc.
Thermal transport deals with the evolution oftemperature over the spatial domain. Reactive
transport considers advection, dispersion, anddiffusion of chemical species while they undergo
biogeochemical reactions. Biogeochemical reactions
include aqueous complexation, adsorption-desorption,ion-exchange,
precipitation-dissolution,
redox,
acid-base reactions, microbial-mediated redox, thenutrient cycle, carbon cycle, metal cycle and biota
kinetics. Geo-mechanics deals with the deformation
or displacement of a solid skeleton.
The aforementioned THMC processes may bestrongly or weakly coupled. Fluid flows are the key
driving mechanism for thermal transport andbio-geo-chemical transport. The variations ofpressure and degree of saturation in fluid phases may
be a significant driving force in geo-mechanic
deformation through the pressure force and gravity.
In contrast, thermal transport, bio-geo-chemicaltransport, and geo-mechanic deformation may affectfluid flow. Temperature changes are induced by fluid
injection, associated phase changes and chemical
reactions. Conversely, temperature changes can affectfluid flow, reactive transport, and geo-mechanics.
Reactive transport is controlled by fluid flow, asaffected by media deformation and temperaturechanges. In contrast, reactive transport will havefeedback to fluid flow, thermal transport, andgeo-mechanical deformation. Geo-mechanicallydriven faults and fractures will affect fluid pressure,thermal transport, and chemical migration.
Conversely fluid pressure changes, thermal transport,
and bio-geo-chemical migration may induce rockdeformation and fault displacement.
Table 1 lists some of the well-known andwidely-used models that have included a single
process or partially or fully coupled THMC processes.
Column 1 lists the model name, while column 2
indicates what media these models can be applied to.
Some models are applicable to only porous media(Macro), some are applicable to multi-porosity media(Macrco, Meso, and Micro), while others areapplicable to porous media with embedded or
imbedded fractures. Column 3 indicates if flow
processes are included. When these are included,
some models are for saturated media (1P) only, somefor variably saturated (VF), while some for multiple
phase fluids (2P, 3P, or MP). There are basicallytwo different ways to model multiphase flow:pressured-based (p-) and fractional flow-based (ff-)
approaches. In general, fractional flow based
approaches are more robust for numerical
discretization. Column 4 indicates if thermaltransport is included in a model, while Column 5shows if a model has integrated commerciallyavailable solid mechanic models with other processes,
or if geo-mechanics has been coupled with partial or
full THC processes.
Column 6 shows (1) the number of chemicalcomponents and the number of species in eachcomponent a model can handle, and (2) how
components are determined. Some models cansimulate the transport of Multiple Species (MS)without introducing the concept of components.
Many models can simulate the transport of Multiple
Components and each component includes Multiple
Species (MCMS). A vast majority of models havepre-determined what the possible components are,and only a few can automatically determine
components based on the input of reaction networks.
However, those models that deal with multiple
species might as well be said to deal with multiple
components.
Column 7 indicates how many THMC processes amodel has been developed for, and these include H, C,
C(R), HM, THM, THC, THC(R), THMC, and
THMC(R). An H-model can only simulate flow
processes. A C-model can only simulate solutetransport or solute transport with adsorption of linear
or nonlinear isotherms. A C(R)-model can simulate
reactive chemical transport. An HM- model cansimulate flow and mechanics. A THM-model cansimulate thermal transport, flow, and mechanics. A
THC- or THC(R)-model can simulate flow, thermaltransport, and solute transport or reactive transport. A
THMC- or THMC(R)-model has the capability tosimulate all processes.
Column 8 indicates what type of biogeochemical
reactions a model can include. These reactions coveraqueous complexation (A/C), adsorption-desorption(A/D), ion exchange (I/E), precipitation-dissolution(P/D), and biological or microbe-mediated reactions(Bio). When a type of reaction is included, thiscolumn also shows if the reaction is treated as anequilibrium reaction (EQ), kinetic reaction (KI), or
optional equilibrium or kinetic reaction (EK). For an

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanequilibrium reaction, a User Defined Algebraic
equation (UDA) rather than the mass action equation
may be used to implicitly define its rate. If a modelallows UDA for all types of equilibrium reactions, it
is labeled all (All). Column 8 also shows if a modelallows Users to Define the Rate equation (UDR) for a
kinetic reaction. If only biological reaction rates areallowed to be defined by users, it is labeled “Bio.”
If a model allows users to define rate equations forany type of kinetic reaction, it is labeled “All.”
Column 9 shows what Approach is employed to
Couple Reactions and Transport (ACRT) in thereactive transport process. Basically, there are five
popular approaches to coupling reactions andtransport: the Direct Substitution Approach (DSA) or
Global Implicit Method (GIM), Sequential Iteration
Approach (SIA), Operator Splitting (OS), Predictor
Corrector, and mixed Predictor Corrector andOperator Splitting (PC-OS). The advantages anddisadvantages of these have been widely discussed in
the literature. In general, the DSA or GIM will taketoo much CPU and RAM for large two- andthree-dimensional problems, but they are perhaps themost accurate methods. The OS, PC, or PC-OS willtake much less CPU than DSA or GIM, but are muchless accurate. The SIA approach offers a compromise
between these two previous kinds of approaches.
Column 9 also shows what form of numerical
discretization is adopted to discretize transport
equations. These methods include Finite Element
Methods (FEM), Finite Difference Methods (FDM)with Total Variation Diminishing (FDM/TVD)
schemes for the advection transport, Integrated Finite
Difference Methods (IFDM), Finite Volume Methods(FVM), and hybrid Lagrangian-Eulerian Finite
Element Methods (LE-FEM). The LE-FEM methods
and TVD are perhaps the methods to resolve theoscillation problems caused by advection transport.
Finally, Column 10 indicates if a Mesh Generator(Mesh Gen) and Users Graphical Interface (GUI) arebuilt into a model to facilitate its use with real-world
problems.
2. Governing Equation of THMC models
It is believed that the best way to address the fully
coupling issue of THMC processes is via thegoverning equations that are derived based on theprinciples of the conservation of mass, momentum,
and energy. The explicitly coupled governing
equations are given in the following subsections.
2.1 Flow Equations. The governing equations for
multiphase flow are derived based on theconservation principle of phase mass and simplified
momentum equations. The detailed derivations canbe found elsewhere (Yeh and Tsai, 2015a, 2015b).

Only the governing equations will be summarized inthis paper. The governing equations of multiphaseflow can be formulated using either the pressurebased or fractional flow-based approach.
Pressure-Based Flow Equation The pressure-based
governing equations of multiphase flow can be statedas follows (White and Oostrom, 1995; White et. Al.,
1995; Xu and Pruess, 1998; Xu et al., 1998, 2003)
 S

dSˆ dpˆ dTd ˆ   ˆ S  S dtdt
p dtT dt   SiM  d   u ˆ dCi ˆ S
Ci dt
dt

 ˆ k    r k   p   gz    M   R ,    L ;
 


(2.1)


  s  1; M    M i , R   Ri ; S   ,iM  iM    e  r , e   e , r   r ,     .
 L

 L

 L

where  is the porosity; S is the degree ofsaturation of Phase  ; ˆ is the molar density of

Phase  ; p is the pressure of Phase  ; t is
the time; T is the temperature; M   having M members is the set index of all global species that arein Phase  , where M  is the number of species

in Phase  ; Ci is the molar concentration ofSpecies i in Phase  ; u is the displacement;kr is the relative permeability of Phase  ;  isthe dynamic viscosity of Phase  ; k is theintrinsic permeability tensor;  is the density of
Phase  ; g is the gravity constant; z is thepotential head; M  is the artificial mass rate ofPhase  ; R is the mass production rate of Phase due to reactions;  L is equal to 1, 2, , L inwhich L is the number of fluid phases; s is thevolume fraction of the solid phase; M i is thesource/sink rate of Species i in Phase  ; Ri isthe production rate of Species i in Phase  due toreactions;  is the porosity of Phase  ; e isthe effective porosity of Phase  ; r is theresidual porosity of Phase  ; e is the effective
porosity; and r is the residual porosity.
The first and second equations in Eq. (2,1) consistof two equations involving 11 unknowns, p ,  ;S , kr , ˆ ,  ; T ; Ci , Ri ; s and u . The
temperature T is governed by the thermal transport
equation. The species concentrations Ci can beobtained from reactive transport equations, whileRi can be obtained via the biogeochemical reaction
module. The solid phase volume fraction s anddisplacement u are governed by the mass and6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanmomentum conservation equations of solid phase in
the geo-mechanic module. If we associate p and with the first and second equations in Eq. (2.1),
we would need four more equations for four
unknowns S , kr , ˆ , and  .
The densities ˆ and  are calculated using
equations of states. In a reactive system, the equationof state for a phase should be derived based on thedefinition of phase density and the conservation of
volume, as follows (Chang, 1995; Cheng, 1995; Yeh
and Tsai, 2015a, 2015b)
ˆ 


 maj





M maj


   majM  maj  
1  i
Ci and
M maj i 
iM  ,i  maj   maj 
1   Ci M i ,   Li iM  ,i  maj 





(2.2)



where i pi , T :   L , being function of partial
phase pressure and temperature, is the intrinsic
density of Species i in Phase  . It is noted thati  i  pi , T  is in fact the equation of state forchemical species i in phase  . It is thus seen that

the equation of state for a phase should be obtained
from the equations of state of all chemical species in
the phase in a reactive system. There are many
specific forms of phase-based equations of state inreactive systems in the literature (e.g., Pruess, 1991;
Finsterle et al., 1994; White and Oostrom, 1995;
White et al., 1995).
The degree of saturation and relative permeability,S and kr , are evaluated using K-S-P(Conductivity-Saturation-Pressure) constitutive laws,
which are posed in general terms as follows
S  S  p  and kr  kr  p  ,  M (2.3)
There are many specific forms of K-S-P equations in
the literature (e.g., Burdine, 1953; van Genuchten,
1980; Parker et al, 1987a, 1987b; Lenhard and Parker,
1998; White and Oostrom, 1995; White et al., 1995;
Yeh and Tsai, 2011; Tsai and Yeh, 2012, 2013).
The pressure-based multiphase flow equation is given
in Eq. (2.1). This equation constitutes L mixed
parabolic and hyperbolic partial differential equations.
The dominance of the hyperbolic or parabolic foreach of L equations varies over space and time.
Furthermore, each equation is highly nonlinear. As aresult, numerical solutions of these equations often
encounter convergence difficulties. Furthermore,cumbersome procedures of variable switch have to beadapted to deal with phase disappearances.
Fractional Flow-Based Flow Equation It is widely
known that using a fractional-flow based approach
could alleviate some of these problems (Allen, 1983;
Abriola and Pinder, 1985a, 1985b; Guarnaccia andPinder, 1997; Suck, 2003; Suk and Yeh, 2008),

because such a method will produce one almost linear
parabolic equation and  L  1 hyperbolic-dominantequations. Hyperbolic-dominant equations can beideally dealt with using hybrid Lagrangian-Euleriannumerical discretization (Yeh, 1990; Yeh, et al. 1995;
Cheng et al., 1998; Suk and Yeh, 2007; Zhang et al.,2008). Phase disappearances were automaticallytaken care of (Suk and Yeh, 2007).
The fractional-flow based approach that originatedin the petroleum engineering literature employs water
saturation, total liquid saturation, and a total pressure
as the primary variables for three-phase flow. For an
arbitrary number of L phases, the total pressure and L  1 accumulated effective degree of saturations
can be employed as the primary variables. Detailed
derivations of fractional-flow based governing
equations for L phases can be found elsewhere(Yeh and Tsai, 2015a, 2015b). In this paper, theseequations are given as follows
S PPt
SeP

L
dPt L 1 P dSt
dCi
d
dT
S
 SP
 STP
   SCPidt  1 St dt
dt
dt
dt
 1 iM  d   u    κ   Pt   gz    M L  R L
dt
in which M L   M  and R L   R L

 L

dP L 1
dS
dCi
d
dT
S Pt t   S St t  S  S  ST
dt  1dt
dt
dt iM  Ci dt

(2.4)



Se

d   u dt L 1    κ Pt      κ    DSt St    1   κ   gz   M   R ,    L  1in whichPt 
1 L 1L  11
 p1  p2  ...  pL1  pL  Lpc 1,   d , p    p  p ;   

0

L

c 1,

 1 1

1

  r ,     for   L  1;  rL Lˆ  kr 
St 0  0 and StL  1; κ  
k ;    S  ˆ  ; 1   1
 St   S  , S   1L       ; and  
 1

where

Pt

(2.5)

ˆ kr L ˆ  kr /.
  1  is the total pressure,St is theaccumulated effective saturation up the  -thphase, κ is the total mobility tensor,  is themobility-weighted density, M L is the mass
source/sink to all L fluid phases, R L is the mass
production rate in all L fluid phases due toreactions,  is the coefficient of fractionalmobility for the  -th phase, DSt is the diffusion

coefficient of the accumulated degree of saturation,pc 1, is the capillary pressure between the relative
non-wetting phase   1 and wetting phase  ,S  is the effective degree of saturation of the  -thphase, and 
density.

is the saturation-weighted molar

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, JapanP
In Equation (2.4), S Pt and S Pt are the storageP
coefficient of total pressure; S St and S St are thestorage coefficient of accumulated effectivePsaturation; S and S are the storage coefficientof porosity; STP and ST are the storage coefficientSof temperature;

P
Ci

and



SC 

are the storageicoefficient of species; and SeP and Se the storage
coefficient of dilation. The superscript P in thesestorage coefficients denotes their appearance in the
total pressure equation while the superscript 
denotes their appearance in the  -th massconservation equation.
2.2 Thermal Transport Equation. Assume thermal
equilibrium in the system, the energy per phasevolume can be related to the temperature via thespecific capacity. In the meantime, the dispersivediffusive-conductive heat flux is related to thegradient of temperature. Invoking the principle ofenergy conservation, we can derive the thermaltransport equation as follows (Yeh and Tsai, 2015a,
2015b)STdT
 L
 KT T      V Ch T  dt  1 L     S DT  kT δ  s k Ts   T   S Hout  S Hin  S R
  1






(2.6)

in whichLST    S Ch  s Chs , KT 
 1

ST
 ST   Vs ,tk N given by H    jk  H     jk  H  &jM objjM of H  is the referenced formation enthalpy of the
o
b

j

j ;

product species jk

is the stoichiometriccoefficient of the reactant species j in the k -threaction; and H  is the referenced formationofjenthalpy of the reactant species j . It is noted that,similar to the Van’t Hoff relationship for theequilibrium
constant,
we
have
used
F T   exp H ko T  To  /  RTTo  .
2.3 Reactive Transport. When the unit of species
concentration is in moles per unit phase volume, themass conservation equation for any chemical species
is given as follows (Yeh and Tsai, 2015a, 2015b)
 Ci    V Ci      Vs Ci  t   ˆ   D    am δ   i   M i  Ri , i  M   ,   L ; S D  aT V δ   aL  aT  V V / V ;(2.8)

s Cis   s Vs Cis      ˆ ss ams δ is   Ris , i  M s .
twhere Ci is the molar concentration of Species iin Phase  ; D is the dispersion coefficient tensor
for all species in the  -th phase;   is the
tortuosity for Phase  ; am is the molecular
diffusion coefficient of all species in the  -thphase; δ is the Kronecker delta tensor;  i is themole fraction Species i in Phase  ; aL is theS R    H k rk , H k  H ko F T  with H ko & F T okF T  is an empirical function of temperature;

(2.7)

j

 H ko  1 1  
F T   exp   
 R  To T  

longitudinal diffusivity for Phase  ; aT is thetransverse diffusivity for Phase  ; Cis is the molar
concentration of Species i in the solid phase; ˆ swhere ST is the storage coefficient of temperature;
KT is the effective decay of the temperature; V isis the molar density of the solid phase; ams is themolecular diffusion coefficient of all species in thesolid;  is is the mole fraction Species i in thethe Darcy velocity of the fluid phase  ; Ch is thesolid phase; Ris is the production rate of the solid

heat capacity of the fluid phase  ; DT is theapparent thermal dispersion coefficient tensor of the -th phase; kT is the apparent thermal diffusion
coefficient of the  -th phase; k Ts is the apparent
thermal conductivity tensor of the solid phase; S Houtis the energy rate artificially taken out of the system;
S Hin is the energy rate artificially injected into thesystem; Vs is the velocity of the solid skeleton; S R
is the energy rate produced due to reactions;  H k is
the intrinsic reaction enthalpy of the k -th reaction;
 H ko is the referenced reaction enthalpy of thek -th reaction;  jk is the stoichiometric coefficientof the product species j in the k -th reaction;

species i due to reactions; and M s  having M smember is the set index of all global species that arein the solid phase, where M s is the number of
species in the solid phase.
A reaction-based formulation would result in thefollowing equation for Ri and/or Ris as follows
Ri    ik  ik  rk , i  M  ,   L  1
(2.9)
kN 

in which

rk  f k Ci : i  M  ,   L  1 , k   N K .

    C  /   C   0, k   N .

rk : Fk Ci : i  M  ,   L  1  0, for example,


e
k



iM

i

 ik



iM

ik

i

(2.10)E

Ci is related to Ci by Ci  M i Ci .

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanwhere rk is the rate of the k -th reaction; Ci isthe molality of the k -th species;  N K  having
N K members is the set index of kinetic reactions; ke is the modified equilibrium constant of the k -threaction; and  N E  having N E members is the set
index of independent equilibrium reactions. Note:  L  1 denotes the solid phase; N  N E  N K .
Equations (2.8) through (2.9) constitute a system ofM transport equations and N E nonlinear algebraic
equations for M Ci 's and N E rk 's of unknowns,
after N K rk 's in Eq. (2.10) is substituted into Eq.
(2.8). Simultaneous solutions of this system of M
reactive transport equations and N E nonlinear
algebraic equations over the entire domain of interest
would demand excessive computational time andcomputer storage (Yeh and Tripathi, 1989), even withtoday’s computers. Decoupling of the N E rk 's of
unknowns from the M Ci 's of unknowns can bedone via the Gauss-Jordan reduction of reaction
networks (Chilakapati, et al., 1998; Fang, et al., 2003;
Kräutle and Knabner, 2007). Performing theGauss-Jordan reduction of reaction networks, we
obtaind   u dEi   VEi Ei 
Ei    DEi Eidt
dt
ss
 S Ei  REi , i  M C  M KI 







in which

Ei   aij C j , Ei jM  a  C  ,jM  ij(2.12)

s s
i m

Ei
Ea
D 
δ, S Essi   aij M j , andEi L Ei
jM 

D Ei  

REi  Dik rk j N KD ( k )Dij rjwhere Ei is the i-th transformed variable; VEi
is the advection velocity of

Ei , D Ei

is thess
dispersion-diffusion coefficient tensor of Ei , S Ei
is the artificial source/sink rate contributing to thechange of Ei ; REi is the production rate of Ei

due to reactions,components,

M C 

is the set index ofM KI  is the set index of kinetic

variables, aij is the decomposed unit matrix, Eia
is the content of Ei in the  -th phase, Eis is
the content of Ei in the solid phase, Dik is thei-th row k -th column of the decomposed reaction
matrix, andN

KD ( k )ˆ s dps  dsˆ dTˆ dCis s s  s ss ˆ s  s





p
d
dt
T
dt
Ci dt
iM s ss d   u  ˆ ss
 R s in which R s   Ris
dtiM s (2.13)ands  s

d 2u
0
dt 2



  T     S p       S  s  s  gz  Tg T  L  L  is Cis  iM s (2.14)   C  L iM  iiin whichT  s ps I  s1B  s2B 1 (trD)I  21D  2 (DB  BD)  3 (DB 1  B 1D)
in which B  FFT , F  I  H, H  u,D

(2.15)

1
H  HT , and H  u  Vs
2

where ps is the solid phase pressure; Rs is themass production rate in the solid phase due toreactions; T is the is the Cauchy stress tensor(positive for tension) exerted on the medium; Tg

is the thermal expansion coefficient in thegeo-mechanical module;  is is the expansion
coefficient of the solid species concentration Cis ;j V EE 
VEi     i  ˆ  S D  ˆ  E  L   Ei
Eis ˆ ss ams δ ,ˆ ss Ei
(2.11)

2.4 Geo-Mechanic Equations. The governing
equations for the geo-mechanic module are obtained
based on the conservation of mass and momentum asfollows (Liu, 2002, 2006, 2011; Liu et al., 2010; Yeh
and Tsai, 2015a, 2015b)

 is the set index of kinetic

reactions that are linearly dependent on the
independent kinetic reaction k .

 i is the expansion coefficient of the  -th species
concentration Ci ; s1 , s2 ,  , 1 ,  2 , and 3
are six model parameters, respectively; F  X x is

the gradient of the deformation x  x in thematerial frame X  X ; H   X u is the gradient of
the displacement u in the material frame X ; andH   X u  Vs is the gradient of the solid phasevelocity in the material frame X .

Equations (2.13) through (2.15) consist of threeequations involving 5 unknowns, ps , s , and  s
as well as u and T . If we associate s , u , andT , respectively, with Eqs. (2.13), (2.14), and (2.15),
respectively, we would need two more equations for
two unknowns ˆ and ps . The density ˆ s is
calculated using the equation of state similar to thosefor fluid phases as follows(Chang, 1995; Yeh andTsai, 2015a, 2015b)s
 s  majs maj s
1  s Ci M iiM s ,i  maj i 

(2.16)



where is pis , T , being function of partial phasepressure and temperature, is the intrinsic density ofSpecies i in the solid phase. It is noted that

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanis  is  pis , T  is in fact the equation of state forchemical species i in the solid phase. The
constitutive equation of the solid phase pressure
versus volume fraction of solid may be established as
ps  f s  in general, e.g., ps  pso 1 1  s n
 1  so (2.17)

3. A Demonstrative Example of THMC ModelingThis example is to demonstrate the capability of
HYDROGEOCHEM 6.1 (Yeh and Tsai, 2015a) tosimulate fully coupled THMC processes. Radioactive
wastewater containing a high concentration of
NpO2 is injected into the media filled with athree-fluid phase (water, NAPL, and air) system(Figure 1). The media contains Region B with a zeroadsorption site, TOTSOH  0 , and Region A with arelatively high adsorption site, TOTSOH  104.
m o l/ d3m
Figure 2 Initial and Boundary Conditions andArtificial Source for Flow Simulation

Z=0
Injection well(8.25, -3.25)
Region A
4.125

3.453125
4.125
2.03125
Z=-6.5
X=0

Region B
X=16.5Figure 1 Region of Interest and Location of
Injection Well
For this problem a three-phase flow is considered,and thermal effects are included. Injection effects on
geo-mechanics (porosity change and deformation) arealso considered, while the chemistry includes both
intra (homogeneous) and inter phase (heterogeneous)reactions. The initial and boundary conditions forflow simulation are depicted in Figure 2. The initialdegree of saturation for aqueous, NAPL, and gaseous
phases, respectively, are 0.1, 0.1, and 0.8 at boundary
B4, and vary with a slope of 0.121212 % until theyreach boundary B2. While the degrees of saturation
will change due to a 0.1 dm3 day-1 of water injection
rate, the degrees of saturation on all four boundaries
are kept the same as the initial conditions.
In order to model thermal transport, the initial
temperature is set at 298 oK. The temperatures onfour boundary sides are maintained at 298 oK, whilethat in the injection well is 308 oK (Figure 3).
For reactive chemical transport, a reaction network
consisting of 28 reactions involving 42 species isconceptualized, as shown in Table 2.

Figure 3 Initial and Boundary Conditions andArtificial Source for Thermal Transport
The initial and boundary conditions for thishypothetical simulation of reactive transport are
depicted in Figure 4. The initial and boundary
conditions for three major species, H2O, NAPL, andAir, are identical, and the same is true for the species
O(g), OH(g), CH4(g), and HCO(g). The initial andboundary conditions for two aqueous components
Ca2+ and CO32- are also identical. The major
disturbances to the reactive system are due to theinjection of a high concentration of NpO2+, highinitial concentration of NpO2+ in Region A, as well as
high adsorption site in Region A.
For the geo-mechanical simulation (Figure 5), it isinitially assumed there is no displacement throughout
the region of interest. On boundary BC1, no external
force is applied. On boundaries BC2 and BC4, it isassumed that the displacement in the horizontal
direction and force in the vertical direction are zero,as depicted with rollers on these two boundaries. On

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanboundary BC3, vertical displacement is not allowed
and the horizontal force acting on the surface is zero,
which is represented by rollers attached to theboundary.

Figure 4 IC & BC for Reactive Transport

It is seen that water entering into the simulation
region displaces the two other fluids (Figure 7). It isalso seen that NAPL is displaced by the injection
water and moves downward (Figure 8). On theother hand, air is displaced by the injected water andmoves upward (Figure 9). These interactions among
fluid phases are easily seen in Figure 10, where thestream lines and velocity vectors are shown.

Figure 7 Distribution of Degree of Saturation
at Various Times for the Aqueous Phase

Figure 8 Distribution of Degree of Saturation
at Various Times for the NAPL Phase
Figure 5 IC and BC for Geo-Mechanics
For numerical simulations with HGC 6.0 (Yeh andTsai, 2013; Tsai and Yeh, 2014), a uniform mesh of
45 x 333 =1,485 nodes and 44 x 32 = 1,408 elements
is used (Figure 6). The total simulation time is 2.985
days. The initial time step size is 0.0001 day. The
maximum time step size is 0.001 day, and total
number of time steps is 3,000.

Figure 9 Distribution of Degree of Saturation
at Various Times for the Gaseous Phase

Figure 6 FE Mesh: 1,408 Elements 1,485 Nodes

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

Figure 13 Distribution of Total Dissolved NpO2+
in the NAPL Phase at Various Times

Figure 10 Velocity Field at Time = 2.98 Days
Temperature decreases with the distance from theinjection well and exhibits a non-symmetrical plumetoward the right-upper direction (Figure 11). Thisbehavior is caused by the non-symmetrical streamlines shown in Figure 10.
Figure 14 Distribution of Total Dissolved NpO2+
in the Gaseous Phase at Various Times
NpO2+ is adsorbed when the fronts of the plumes inall three phases reach the adsorption site (Figure 15).

Figure 11 Distribution of Temperature
at Various Times
The total dissolved NpO2+ in all three phases istransported away from the injected well and retarded
when reaching the adsorption region (Figures 12, 13,
and 14).

Figure 15 Distribution of Total Adsorbed NpO2+
at Various Times
The distribution of pH is in the range of 10.50 ~10.71. It is seen that pH is relatively low in Region A
compared to Region B (Figure 16). It is further seenthat pH increases along with the approach of theinjected NpO2+. This is because NpO2+ and H+ arecompeting for the same sites.

Figure 12 Distribution of Total Dissolved NpO2+
in the Aqueous Phase at Various Times

The total concentration of CH4(g) in the gaseous
phase is shown in Figure 17. Although the initial andboundary conditions for CH4(g) are the same, the

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanevolution of total for CH4(g) is due to the complex
inter-phase reactions that occur.

models center around their increasing design
capability with regard to representing these coupledprocesses: from the simplest one-phase groundwaterflow to the most complete of the aforementioned
processes. Many fully or partially coupled THMC
models or single process models have beendeveloped for practical applications. However,current models are still incomplete, sinceelectro-magnetic waves have not yet been coupledwith THMC processes. A hypothetical example was
employed in this study to illustrate simulations offully coupled multiphase flow, thermal transport,
reactive transport, and geo-mechanics. Extensive
test of fully coupled THMC models is expected in the
future (DECOVALEX 2019).

Figure 16 Distribution of pH at Various Times
2.5 References
Abriola, L. M. and G. F. Pinder, 1985a. A multiphaseapproach to the modeling of porous media
contamination by organic compounds, 1. Equation
development.
Water
Resources
Research
1985;21(1):11-18.

Figure 17 Contour of Total CH4(g)at Time = 2.98 Days
The change in porosity results from both theboundary conditions imposed on the boundaries B2,
B3, and B4, and the free displacement on boundary
B1. Internal changes of pressure in all three phases,
which are caused by changes in the degree ofsaturation and total pressure, are the major factors in
the evolution of porosity (Figure 18).

Abriola, L. M. and G. F. Pinder, 1985b. A multiphaseapproach to the modeling of porous media
contamination by organic compounds, 2. Numerical
Simulation,
Water
Resources
Research,
1985;21(1):19-26.
Allen, M.B., 1983. Collocation techniques formodeling compositional flows in porous media. Ph.D.
Dissertation. Dept. Of Civil Eng., Princeton
University, Princeton, New Jersey (1983). Dept. Of
Civil Eng., Princeton University, Princeton, New
Jersey (1983
Burdine, N.T., 1953.
Relative permeability
calculations from pore-size data. Trans. A.I.M.E. 198,
71-77, 1953.
Chang, Jing-Ru C., 1995. Numerical Modeling ofThree-Dimensional Subsurface Flow, Heat Transfer,and Fate and Transport of Chemicals and Microbes.
PhD Dissertation. Department of Civil andEnvironmental Engineering, The Pennsylvania State
University, University Park, PA. 253 pp.

Figure 18 Distribution of Change of Porosity at
Various Times
4. Conclusion and Summary
The coupled groundwater processes include fluid
flow in either variably saturated media or multiple
phases, thermal transport, geo-mechanics, andreactive transport. Advances in the related numerical

Cheng, Hwai-Ping, 1995. Development andApplication of a Three-dimensional Finite Element
Model of Subsurface Flow, Heat Transfer, andReactive Chemical Transport. PhD Dissertation.
Department of Civil and Environmental Engineering,
The Pennsylvania State University, University Park,
PA. 357 pp.
Cheng, H. P., J. R. Cheng, and G. T. Yeh, 1998. A
Lagrangian-Eulerian Method with Adaptively Local

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, JapanZooming and Peak/Valley Capturing Approach toSolve
Three-Dimensional
Advection-Diffusion
Transport Equations.
International Journal for
Numerical Methods in Engineering, 41:587-615.
Chilakapati, A., T. Ginn, and J. Szecsody, 1998. An
Analysis of Complex Reaction Networks in
Groundwater Modeling, Water Resources Research,
Vol. 34, pp. 1767-1780, 1998.
DECOVALEX 2019. DEvelopment of COupledmodels and their VALidation against EXperiments,
2016-2019.
https://www.quintessa.org/latest-news/preparation-for-the-new-decovalex2019-project.html
Fang, F. L., G. T. Yeh, and W. D. Burgos, 2003. A
New
Paradigm
toModel
Reaction-Based
Biogeochemical Processes, Water Resources
Research, Vol. 39, pp. 1083-1108, 2003.
Finsterle, S, G. J. Moridis, and K. Pruess, 1994. A
TOUGH2 Equation-of-State Module for the
Simulation of Two-Phase Flow of Air, Water, and a
Miscible Gelling Liquid. LBL-36086, UC-400. Earth
Sciences Division, Lawrence Berkeley Laboratory,
University of California, Berkeley, CA 94720
Guarnaccia, J. F and G. F. Pinder, 1997. NAPL: A
Mathematical Model for the Study of NAPL
Contamination in Granular Soils, Equation
Development and Simulator Documentation. The
University of Vermont, RCGRD #95-22, 1997.

Computational and Applied Mathematics, 29(3):465–
478
Liu, I-S, 2011. A note on the Mooney-Rivlin material
model. Instituto de Matematica, Universidade Federal
do Rio de Janeiro, 21945-970, Rio de Janeiro, Brasil
Parker, J. C., Lenhard, R. J., and T. Kuppusamy,
1987a: A parametric model for constitutive properties
governing multiphase flow in porous media. Water
Resour. Res., 23(4), 618--624.
Parker, J. C., and R. J. Lenhard, 1987b: A model forhysteretic constitutive relations governing multiphase
flow: 1. Saturation-pressure relations. Water Resour.
Res., 23(12): 2187--2196.
Pruess, K., 1991. EOS7, An equation-of-state module
for the TOUGH2 simulator for two-Phase flow ofsaline water and air, Lawrence Berkeley Laboratory
Report LBL-31114, Berkeley, CA.
Suk, H., 2003. Development of 2-D and 3-D
Simulator for Three-Phase Flow with General Initial
and Boundary Conditions on the Fractional Flow
Approach. PhD Dissertation. The Department ofCivil and Environmental Engineering, Pennsylvania
State University, University Park, PA. pp. 282.
Suk, H. and G. T. Yeh, 2007. 3D, three-phase flow
simulations using the Lagrangian-Eulerian approachwith adaptively zooming and peak/valley capturing
scheme (LEZOOMPC). Hydrologic Engineering,
ASCE, 12(1): 14-32.

Kräutle, S. and P. Knabner, 2007. A reduction
schemeforcoupledmulticomponent
transport-reaction problems in porous media:
Generalization to problems with heterogeneous
equilibrium reactions, Water Resources Research,
Vol. 43, W03429, 2007.

Tsai, C. H. and G. T. Yeh, 2012. Retention
Characteristics for Systems of Multiple Phase-Fluids.
Terrestrial, Atmospheric, and Oceanic Sciences(TAO),
23(4):451-458,
doi:10.3319/
TAO.2012.0.14.01(Hy).

Lenhard, R. J. and J. C. Parker, 1998. Experimental
validation of the theory of extending two phasesaturation-pressure relations to three phase systems
for monotonic saturation paths, Water Resour. Res.,24, 817–830, 1988.

Tsai, C. H. and G. T. Yeh, 2013. An Advanced
Constitutive Law in Multiphase Flow Model for
Simulations in Compressible Media. Recent advancesin Hydrogeology, American Geophysical Union,
Chap 2, pp.27-56., 2013.

Liu, I-S, 2002. Continuum Mechanics, SpringerVerlag, Berlin-Heidelberg

Tsai, C. H. and G. T. Yeh, 2014. Benchmarking a
THMC Process Problem with an Explicitly Coupling
Model in Subsurface Media. Abstract and Program,
SSBench IV. Chateau de Cadarache, France, October
6 - 8, 2014. pp. 40

Liu, I-S, 2006. A Continuum Mechanics Primer On
Constitutive Theories of Materials. Instituto deMatematica, Universidade Federal do Rio de Janeiro,
Caixa Postal 68530
Liu, I-S, R. A. Cipolatti, and M. A. Rincon, 2010.
Successive linear approximation for finite elasticity.

van Genuchten, M. Th., 1980: A closed-form
equation for predicting the hydraulic conductivity ofunsaturated soils, Soil Sci. Soc. Am. J., 44, 892--898.

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

White, M. D. and M. Oostrom, 1995. STOMP
Subsurface Transport Over Multiphases. Theory
Guide. PNL - XXXX. Pacific Northwest National
Laboratory, Richland, WA 99352
White M. D., M. Oostrom, and R. J. Lenhard, 1995.
Modeling fluid flow and transport in variably
saturated porous media with the STOMP simulator. 1.
Nonvolatile three-phase model description. Advances
in Water Resources 1995;18(6):353-364.
Xu, T. and K. Pruess, 1998. Coupled Modeling of
Non-Osothermal Multi-phase Flow, Solute Transportand Reactive Chemistry in Porous and Fractured
Media: 1. Model Development and Validation.
LBNL-42050. Earth Sciences Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720.
Xu, T., K. Pruess, and G. Brimhall, 1998. An
Improved Equilibrium-Kinetics Speciation Algorithm
for Redox Reactions in Variably Saturated
Subsurface Flow Systems. LBNL-41789. Earth
Sciences Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720.
Xu, T., E. Sonnenthal, N. Spycher, and Pruess, 2003.
TOUGHREACT User’s Guide: A Simulation
Program for Non-Isothermal Multiphase Reactive
Geochemical Transport in Variably Saturated
Geologic Media. Berkeley: Lawrence Berkeley
National Laboratory, 2003.

Geomechanics through Multiple Phase Systems
Version 6.0 (FACTM2D): Theoretical Basis andNumerical Approximation. Graduate Institute ofApplied Geology, National Central University, 300
Jhongda Road, Jhongli, Taoyuan 32001
Yeh, G. T. and C. H. Tsai, 2015a.
HYDROGEOCHEM 6.1: A Two-Dimensional
Model of Coupled Fluid Flow, Thermal Transport,
HYDROGEOCHEMical
Transport,
andGeomechanics through Multiple Phase Systems
Version 6.1 (A Two Dimensional THMC Processes
Model)
Theoretical
Basis
and
Numerical
Approximation. Technical Report. Graduate Instituteof Applied Geology, National Central University.
DOI: 10.13140/RG.2.1.1766.2482
Yeh, G. T. and C. H. Tsai, 2015b.
HYDROGEOCHEM 7.1: A Three-Dimensional
Model of Coupled Fluid Flow, Thermal Transport,
HYDROGEOCHEMical
Transport,
andGeomechanics through Multiple Phase Systems
Version 7.1 (A Three Dimensional THMC Processes
Model). Technical Report. Graduate Institute ofApplied Geology, National Central University. DOI:
10.13140/RG.2.1.3339.1125.
Zhang, F., L. Jiang, and G. T. Yeh. 2008. An
adaptive local grid refinement and peak/valley
capture algorithm to solve transport problems withmoving sharp-fronts J Transport in Porous Media.
72:53–69, DOI 10.1007/s11242-007-9135-2.

Yeh, G. T. and V. S. Tripathi, 1989. “A Critical
Evaluation
of
RecentDevelopments
inHydrogeochemical Transport Models of Reactive
Multichemical Components”, Water Resources
Research, vol. 25, pp. 93-108, 1989
Yeh, G. T., 1990. A Lagrangian-Eulerian methodwith zoomable hidden fine-mesh approach to solving
advection-dispersion equations, Water Resources Res.
26(6): 1133-1144 June 1990.
Yeh G. T., J. R. Chang, H. P. Cheng, and C. H. Sung,
1995. An adaptive local grid refinement based on theexact peak capturing and oscillation free scheme tosolve transport equations. Int. J. Comput. Fluids
1995;24:293-329.
Yeh, G. T. and C. H. Tsai, 2011. Proposing a New
Retention Function for Multiple Phase-Fluids.
Abstract and Program, 2011 AGU Fall Meeting, CD.
December 5-9, 2011
Yeh, G. T. and C. H. Tsai, 2013.
HYDROGEOCHEM 6.0 A Two-Dimensional Modelof Coupled Fluid Flow, Thermal Transport,
HYDROGEOCHEMical
Transport,
and6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

Table 1 List of Partially or Fully Coupled THMC Process Models
0

No
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

1

Model
MODFLOW (USGS)
GMS-FEMWATER (PSU, Army Corps)
SEWAT (USGS)FEFLOW (DHI)
AT123D (ORNL, Seview)
PHREEQC (USGS)
OS3D/GIMRT (PNNL)
CrunchFlow (LBNL)NUFT-C (LLNL)
FEHM (LANL)STORM/STOMP/ECKEChem (PNNL)PARSSIM1 (UT Texas)
RT3D (Auburn Univ)
MT3D (Univ Alabama)
HP1/HP2/HP3 (UC Riverside )UTCHEM
MULTI-FLO/PFLOTRAN (LANL)TOUGH2/TOUGHREACT (LBNL)
MIN3P (Canada)
CORE (Spain)
COMSOL Multiphysics-PHREEQC (SKB)
OpenGeoSys (UFZ, Germany)
FRACHem (Switzerland)
HYTEC-CHESS-R2D2 (France)
COUPLYS (PNC, Japan)
THM (Korea)
MURF-MURT (PSU-ORNL)
COWADW123D (PSU)
HBGC123 (PSU - ORNL)MPS (PSU-UCF)
HYDROBIOGEOCHEM (PSU-UCF)
HYDROGEOCHEM 4.0 5.0 (UCF-ORNL)
HGC 4.1 to 4.6 & 5.1 to 5.6 (NCU-TPC)
HGC 6.0 to 6.2 & 7.0 to 7.2 (NCU)

23(H) Flow

4(T)
Thermal5

6f(C ) Chemical Transport

7(M) GeoNo of
ComponentaePhaseb Approachc Transportd
Components Automated? THMCg included?
Multiple Scales Media
Mechanics
Macro
1P
pNone
None
N/A
N/A
H
Macro
VF
pNone
None
1
N/A
HC
Macro
VF
pNone
None
1
N/A
HC
Macro
VF
pYes
None
MS
N/A
THC
1p
None
None
None
1
N/A
Macro
C
Macro
none
None
None
None
MCMS
No
C(R)
Macro, Fracture
1P
pInput
No
MCMS
No
C(R)
Macro, Fracture
1P
pInputNo
MCMS
No
THC(R)
Macro, Fracture
VF
pYes
No
MCMS
No
THC(R)
Marco, Fracture
3P
pYes
Solid Mech
MS
N/A
THMC
Macro, Fracture
3P
pYes
No
MCMS
No
THC(R)
Macro, Fracture
VF
pYes
No
MCMS
No
THC(R)
Macro
1P
input
Input
No
MS
N/A
C(R)
Macro
1P
pNo
No
MS
N/A
C
Macro
VF
pYes
No
MS
N/A
THC(R)
Macro
3P
ffYes
Cmp
MCSS
No
THC
Macro, Fracture
2P
pYes
No
MCMS
No
THC(R)
Macro, Meso, Fracture
3P
pYes
Solid MechMCMS
No
THMC(R)
Macro, Fracture
VF
pYes
No
MCMS
No
THC(R)
Macro, Fracture
2P
pYes
No
MCMS
No
THC(R)
Macro, FractureVF
pYes
GeoMechMCMS
No
THMC(R)
Macro, FractureVF
pYes
Solid MechMCMS
No
THMC(R)
Fracture
1p
pYes
None
MCMS
No
THC(R)
Macro, Fracture
1P
pYes
No
MCMS
No
THC(R )
Macro, FractureVF
pYes
Solid MechMCMS
No
THMC(R)
Macro, Fracture
2P
pYes
GeoMechNo
No
THM
Macro, Meso, Micro, Fracture
VF
pNo
Cmp
1
No
HC
Macro, FractureVF
pNo
GeoMechNo
No
HM
Macro, Meso, Micro, Fracture
VF
pNo
Cmp
MCMS
No
C(R)
Macro, Fracture
3P
ffNo
No
No
No
H
Macro, Fracture
3P
Input
Yes
Cmp
MCMS
Yes
THC(R)
Macro, Fracture
VF
pYes
Cmp
MCMS
Yes
THC(R)
Macro, FractureVF
pYes
Cmp, GeoMech
MCMS
Yes
THC(R),THMC(R)
Macro, FractureMP
ffYes
GeoMechMCMS
Yes
THMC(R)a

Macro = Macroscale, Meso = Mesoscale, Micro = Microscale, Fracture = Fracture scale; b1P = 1 phase, VF = Variably saturated flow. 2P = 2 phases, 3P = 3 phases, MP =Multiple phase (more than 3 phases); cp- = Pressure-based approach, ff- = fractional flow-based approach, input = Must obtain flow field from other models; dYes = modeled,no = not included, input = temperature is input parameter; eCmp = Compressibility is input, geomechanics not explicitly modeled, No = not included, Solid Mech = Solid
Mechanics, GeoMech = Geo-Mechanics; fMCMS = Multiple components-multiple species, MS = Multiple species (No concept of components), MCSS = Multiple
Components-single species; gC(R) = Reaction-based Chemistry

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

Table 1 List of Partially or Fully Coupled THMC Process Models (Continued)
0

1

8

Model
MODFLOW (USGS)
GMS-FEMWATER (PSU, Army Corps)
SEWAT (USGS)FEFLOW (DHI)
AT123D (ORNL, Seview)
PHREEQC (USGS)
OS3D/GIMRT (PNNL)
CrunchFlow (LBNL)NUFT-C (LLNL)
FEHM (LANL)STORM/STOMP/ECKEChem (PNNL)PARSSIM1 (UT Texas)
RT3D (Auburn Univ)
MT3D (Univ Alabama)
HP1/HP2/HP3 (UC Riverside )UTCHEM (UT Austin)
MULTI-FLO/PFLOTRAN (LANL)TOUGH2/TOUGHREACT (LBNL)
MIN3P (Canada)
CORE (Spain)
COMSOL Multiphysics-PHREEQC (SKB)
OpenGeoSys (UFZ, Germany)
FRACHem (Switzerland)
HYTEC-CHESS-R2D2 (France)
COUPLYS (PNC, Japan)
THM (Korea)
MURF-MURT (PSU-ORNL)
COWADE123D (PSU)
HBGC123 (PSU - ORNL)MPS (PSU-UCF)
HYDROBIOGEOCHEM (PSU-UCF)
HYDROGEOCHEM 4.0 5.0 (UCF-ORNL)
HGC 4.1 to 4.6 & 5.1 to 5.6 (NCU-TPC)
HGC 6.0 to 6.2 & 7.0 to 7.2 (NCU)

Geochemical Reactions
I/E
P/D
Bio
N/A
N/A
N/A
No
No
KI
No
No
KI
No
No
KI
No
No
KI
EQ
EK
KI
No
EK
No
EQ
EK
KI
EQ
EK
No
No
EQ
No
No
KI
No
EK
EK
KI
EQ
KI
KI
EQ
KI
KI
EQ
EK
KI
No
EK
KI
No
KI
No
KI
EK
KI
EQ
EK
KI
No
KI
No
EQ
EK
No
EQ
EK
No
EQ
EK
No
EQ
EK
No
EQ
EK
No
No
No
No
No
No
No
No
No
No
EQ
EQ
KI
No
No
No
EK
EK
EK
No
No
No
EK
EK
EK
EK
EK
EK

9
h

No
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

A/C
N/A
No
No
No
No
EQ
EQ
EQ
EQ
EQ
EQ
EK
EQ
EQ
EQ
No
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
No
No
No
EQ
No
EK
No
EK
EK

A/D
N/A
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EK
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
EQ
No
No
No
EQ
No
EK
No
EK
EK

10
i

UDA
N/A
N/A
N/A
N/A
N/A
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
None
All
None
All
All

UDR
N/A
N/A
N/A
N/A
N/A
All
None
All
None
None
None
None
Bio
Bio
All
None
None
All
None
None
All
All
None
None
All
None
None
None
None
None
All
None
All
All

Numerical Methods
ACRTTransport
N/A
N/A
N/A
FEM, LE-FEM
N/A
FEM
N/A
FEM
N/A
AnalyticalOS
Mixed CellOS, DSA(GIM), SIA
FDM
OS, DSA(GIM)
FDM/TVDDSA(GIM)
IFDM
OS
FEM
OS
FDM/FVM
OS
FDM
OS
FDM
N/A
FDM/FVMOS
FEM
N/A
FDM/TVDOS
FDM/FEMOS, SIA
IFDM
SIA
FVM
SIA
FEM
OS
FEM
OS
FEM
OS
FEM
OS, SIA
FVM
OS
FEM
N/A
FEM
N/A
FEM, LE-FEM
N/A
FEM
SIA
FEM, LE-FEM
N/A
FEM, LE-FEM
OS, PC, PC-OS, SIA FEM, LE-FEM
OS, PC, PC-OS, SIA FEM, LE-FEM
OS, PC, PC-OS, SIA FEM, LE-FEM
OS, PC, PC-OS, SIA FEM, LE-FEM

GUI and
MeshGen
Yes
Yes
No
Yes
Yes
Yes
No
No
No
No
No
No
No
Yes
Yes
No
Yes
Yes
No
No
Yes
Yes
Yes
Yes
Yes
No

No
No

No
No
No

No
No
No

h

A/C = Aqueous Complexation, A/D = Adsorption-Desorption, I/E = Ion Exchange, P/D = Precipitation-Dissolution, Bio = Biomediated/Biodegradation, UDA = Allow
User’s Defined Algebraic Equations for Equilibrium Reactions, UDR = Allow User’s Defined Rate Equations for Kinetic Reactions, KI = kinetic reaction, EK =Equilibrium or kinetic reaction, No = not modeled, EQ = Equilibrium reaction, BIO = User’s defined rates are allowed only for biomediated reactions, All = User’s defined
rates are allowed for all reactions; iACRT = Approach to Coupled Reaction and Transport, OS = Operator splitting, DSA(GIM) = Direct substitution approach(Global implicit
method), PC = Predictor-corrector, SIA = Sequential iteration approach, FDM = Finite difference method, FVM = Finite volume method, IDFM = Integrated FDM, FEM =Finite element method, TVD = Total Variation Diminishing, LE-FEM = Lagrangian-Eulerian FEM.

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

Table 2 Reaction Network: 28 Reactions, 42 Species

Aqueous Phase Reactions

No

Equilibrium Constants

H 2O(l)  H + OH

R1

Log K1  14.00

Ca 2+ + CO32-  CaCO3l 

R2

Log K2  3.22

Ca 2+ + H + + CO32-  CaHCO3+

R3

Log K3  11.43+

R4

Log K4  12.85

H + + CO32-  HCO3-

R5

Log K5  10.32

 H 2CO3

R6

Log K6  16.67

NpO+2 + H 2Ol  H +  NpO2  OH 

R7

Log K7  8.85

NpO+2 + CO32-  NpO2  CO3 

R8

Log K8  5.60

R9Log K9  3.5

Gaseous Phase Reactions

No

Equilibrium Constants

O(g) + CH 4(g)  OH(g) + CH3(g)R10

Log K10  0.45

O(g) + HCO(g)  OH(g) + CO(g)R11

Log K11  21.41

Inter Phase Isotherm Reactions

No

S = Kd C

 NpO

R12

Kd, 12  0.2

NpO+2 l   NpO+2  g 

R13

Kd, 13  0.3

H 2O  H 2O NPAL

R14

Kd, 14  0.001

H 2 O  H 2 O g 

R15

Kd, 15  0.01

NAPL  NAPLl

R16

Kd, 16  0.01

NAPL  NAPLgR17

Kd, 17  0.001

Air  Airl 

R18

Kd, 18  0.0001

Air  Air NPAL

R19

Kd, 19  0.0001

+

Ca

2+

-

 H + CaOH
+

+

23

2H + CO

NpO+2 + 2CO32-  NpO2  CO3 NpO+
2 l 

3
2+
2  NAPL 

Inter-phase Kinetic Reactions

No

H 2O  H 2Ok NPALR20

H 2 O  H 2O k  g 

R21

NAPL  NAPLkl

R22

NAPL  NAPLkg

R23

Air  Airkl 

R24

Air  Airk NPAL

R25

Reaction rates

R 20 =10-5  H 2O  10-2  H 2Ok(NAPL) 

R 21 =10-5  H 2O  10-3  H 2Ok(g) 

R 22 =10-5  NAPL   10-2  NAPLk(l) 

R 23 =10-5  NAPL   10-2  NAPLk(g) 
R 24 =10-5  Air   10-1  Airk(l) 

R 25 =10-5  Air   10-1  Airk(NAPL) 

Adsorption-Desorption Reactions
=SOH  H + =SO-=SOH + H   =SOH +2

NpO-2 + H 2O + =SOH  H  +  NpO2   OH  =SOH 

No

Equilibrium Constants

R26
R27

Log K26  10.30

R28

Log K28  3.5

Log K27  5.40

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan