Personal Assignment: Reactive Transport, CIE4365 Coupled Processes
(Deadline Tuesday June 22, 2021)
6th June 2021
1 Coupling Richards Equation, Heat Flow and Reactive Solute Trans-
port
Write a MATLAB or Python code where you expand the code for unsaturated water and heat flow problem with a
code for reactive transport. We will solve the simplified problem of microbial mineralization of glucose in soil that
was introduced during the lecture (see lecture slides Lecture_biogeochemistry_2020.pdf). You need to implement
the problem, write a short report including a summary of the important equations, your assumptions and the
scenarios you implemented.
Please note, this is one of the individual assignments for the final grade.
2 Approach
In order to solve this problem we need to consider the states which will change during the simulation. When
modeling reactive transport all chemical species involved are states, so for this problem, in addition to the pressure
head (hw) and temperature (T), we also need to consider the concentrations of oxygen (O2), glucose (C6H12O6),
bacterial biomass (X or C5H9O2.5), CO2 and water (H2O). In order to simplify the problem we will initially assume
that only oxygen is transported with the water. Both glucose and biomass are immobile and present, dissolved in
the water phase, everywhere in the column. As a consequence, the total mass of glucose and biomass present in the
nodes can change if the amount of water decreases, so the first version of the model should be implemented with
initial and boundary conditions that ensure that glucose concentrations remain constant. The initial concentration
of oxygen in the column is zero and it enters the column dissolved in the water entering the soil column at the
top. We initially assume that the growth and death rates of biomass are equal so the dissolved concentration of of
living biomass remains constant (steady-state). Glucose is consumed and therefore depletes in time. As we are not
interested in the evolution of the CO2 concentrations in time we disregard this state. Therefore, we only need to
simulate the rate of change in the oxygen and glucose concentrations.
In order to program this problem it is important to increase the complexity of the program step by step, while
ensuring correct functioning of the code at every step. I suggest you follow the following steps:
- Use the implementation of the heat flow as a template for the transport single solute (one extra state)
without reaction. Please note that the diffusion constant is temperature dependent (see section on equations).
Implement the code with some very simple boundary conditions so you can check the functioning of the code
(constant flow, and temperature, fully saturated column). - Expand the solute transport code with extra solutes (i.e. extra states). Please make this expansion as general
as possible so adding extra solutes is as easy as possible. - Include reaction terms in order to add the reaction between the different compounds (in this case you start
by calculating the rate of change in oxygen, and from this the rate of change in glucose. If you want you can
also include the rate of change in bacterial biomass and water content as well. - Use the complex derivative version of the Jacobian matrices. Check the number of non-zero diagonals and see
if they indeed fit your problem. Please note that some times the Jacobian contains some very small numbers
due to round off errors, sometimes you have to remove this before analyzing the Jacobian. You can plot the
position of the non-zero values in a matrix with the matplotlib or matlab function spy.
1 - Finalize the code by including the enthalpy change in the heat flow code.
- Run a range of scenarios which you create to test certain conditions.
- Equations
As is given in the slides, the water balance, heat balance and solute balance equations are:
?θw
?t
= C ′(hw)
?hw
?t
= ?? · qw +Rw
ζb
?T
?t
= ?? · qH + θwRO2RH ? T
?ζb(θw)
?θw
?θw
?t
θw
?Cα
?t
= ?? · qSα + θwRSα ? Cα
?θw
?t
where
qw = K
w
satkrw(?hw +?z)
qH = ?λ(θw)?T + qwζwT
qSα = ?Dα(θw, T)?Cα + qwCα
where the parameters and variables for the water and heat balance equation are the same as in the previous
assignments, Cα is the concentration in the water phase of a solute α, Dα(θw, T) is the diffusivity of solute α in the
water phase which is a function of water content and temperature.
In this problem we add local sink-source terms (Rw, RH , and RSα which includes RO2 ,) to the balance equations
which are dependent on the chemical stoichiometry and Monod-functions:
RO2 = ?5.4
μmax(T)
Yxs
CX
CO2
CO2 +KS
Seff
RC6H12O6 =
1
5.4
RO2
RH2O =
?5.55
5.4
RO2
RH = ?2800 [kJ/molO2]
The diffusion constant for the solutes depends on water content and temperature dependent and can be described
with an Arrhenius equation:
Dα(θw) =
D0θw
τ
exp
(?Ea,α
RT
)
where D0 is the diffusion constant, Ea,α is the activation energy of diffusion, and τ is the tortuosity factor accounting
for the increased path length in a tortuous porous system and R is the gas constant [8.3145 J/(mol K)]. The maximum
growth rate also depends on temperature according to an Arrhenius equation:
μmax(T) = μmax25 exp
(
64
R
(
1
25
? 1
T ? 273.15
))
where μmax25 is the maximum growth at 25oC.
2 - Parameter values
You can use the following parameters for the model:
parameter unit value
CX mol/L 0.001
CO2 inflow mol/L 0.630× 10?3
CC6H12O6 initial mol/L 0.05
molar mass H2O g/mol 18
μmax25 mol/(L day) 1
Ks mol/L 0.1
D0,O2 m
2/day 0.3034
Ea,O2 J/mol 18.23× 103
τ [-] 0.5
Yxs molX/molgluc 0.1
Of course you can also vary these parameter in the final step where you test different scenarios. - Boundary conditions
Ensure that your model works for all boundary conditions you applied in the previous assignments. Clearly explain
boundary conditions you used to test the assignment.
The final problem is a scenario where you start with a fully saturated column with the above described initial
conditions (no oxygen, biomass and glucose present in the liquid phase). Then the column is allowed to drain while
water is added to the top of the column at a rate of 5 mm/day. This water has an oxygen concentration of 10 mg/L
(or 0.630 × 10?3mol/L). Assume your soil is a 1 meter column of fine sandy soil with a capillary rise of about 0.5
m. The final steady state is in equilibrium with a water table at 0.5 m below the bottom of the column.
WX：codehelp