(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:

1. 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).

2. 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.

3. 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.

4. 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

5. Finalize the code by including the enthalpy change in the heat flow code.

6. Run a range of scenarios which you create to test certain conditions.

3 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

4 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.

5 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.