CO2 Dissolution in Water
ENGF0003 Project
24-25
Guidelines:
• Type your project in Word or LaTeX. Follow UCL Accessibility Guidelines to format your document. Include a table of contents, page numbers, and use built-in styles (Heading 1, Heading 2) to structure your document.
• All figures and tables must be numbered and contain informative captions. All the main equations throughout your work must be numbered and typed appropriately.
• Submit a single PDF document. Do not write down your name, student number, or any information that might help identify you in any part of the project. Do not copy and paste the coursework brief into your submission – Rewrite information where necessary for the sake of your argument.
This project counts towards 30% of your final ENGF0003 grade.
Introduction
In your ENGF0003 coursework you have taken a data-driven approach to studying ocean acidification via summarising, describing, visualising and generalising data into mathematical models.
In this project, you will work with two theoretical models of the dissolution of CO2 in water, implement them, and discover how theory compliments real-world data in engineering mathematics.
1.1 Phase Equilibrium
During your ENGF0003 journey,you have learned about stationary points, which are those where the derivative of a function is zero and the function does not change with time. The mathematical model of CO2 equilibrium in water is similar, where we assume that variables such as temperature and pressure do not change with time, or change so slow that we can say that their time derivative is sufficiently close to zero.
To create a mathematical model of the solubility of carbon dioxide in seawater, we start from a law of physics known as Henry’s Law. Henry’s Law states that the amount of gas dissolved in a liquid at constant temperature increases as the pressure of the gas above the surface of the liquid is raised.
In this project, two main simplifications will be made to model surface ocean water, the real-world system we wish to model:
i. We will assume that surface seawater behaves like pure water. This assumption is made because there are well-documented empirical relations for CO2 dissolution in pure water.
ii. Although the atmosphere is composed of H2O vapour, N2, O2 , Ar, CO2 , Ne, He, CH4 , Kr, H2 , NO, Xe, O3 , I2 , CO and NH3 , we will focus on modelling a system that is formed only of H2O and CO2 .
1.1.1 Henry’s Law
Figure 1 shows a closed container with a gas and a liquid phase, where the gas phase is represented in orange and the liquid phase is represented in blue. Suppose that this system represents a mixture of water and CO2 both at temperature T [K] and total pressure P [Pa].
• In the gas phase, there is a mixture of water vapour and CO2 gas.
• In the liquid phase, there is a mixture of liquid water and CO2 particles dissolved into the water.
Figure 1. Closed container where a mixture is in phase equilibrium.
Henry’s law states that the following relations are valid in equilibrium (Carey, 1988; Carroll, Slupsky and Mather, 1991):
x1 pv, 1 = y1 φ1 p, (1.1)
x2 H21 = y2 φ2 p. (1.2)
In equations 1.1 and 1.2 xi is the molar fraction of the component with index i in the liquid phase, and yi is the molar fraction of the same component in the vapour phase. Water is the solvent of the system, being represented by an index i = 1 and carbon dioxide is the solute represented by an index i = 2.
• We will model the gas phase in terms of the mass fraction of CO2 in gas formy2 and the mass fraction of water vapour y1 .
• Likewise, we model the mass fraction of water in liquid phase as x1 and the mass fraction of CO2 dissolved in water as x2 .
This results in a system of equations such as
Finally, Pv,1 is the vapour pressure of water described empirically by Equation 3 and H21 is the temperature-dependent Henry’s constant described empirically by Equation 4.
The coefficients φ1 and φ2 are the fugacity ratio in vapour and liquid phases, denoted as φ1 = φl1/φv1 and φ2 = φ . The fugacity ratios allow us to approximate non-ideal gases more accurately while using idealised equations. Appendix A summarises the approach of Peng and Robinson (1976) in obtaining these quantities.
1.1.2 Vapour Pressure of Water
The vapour pressure of water is represented as Pv,1 . This quantity can be computed with the empirical approximation given by the International Association for the Properties of Water and Steam IAPWS (Wagner and Pruß, 2002) as
where x = 1 − Tr , and Tr = T/Tc is a non-dimensional temperature variable and values for αi in Eq. 3 can be found in Table 1. In Eq. 3, Tc = 647.096 [K] is the critical temperature of water, and Pc = 22.064 [MPa] is the critical pressure of water.
Table 1. Coefficients αi for Equation 3
1.1.3 Henry’s Constant
H21 is the Henry constant for the system H2O + CO2 . An empirical temperature- dependent expression for H21 (T) is given by (Carroll, Slupsky and Mather, 1991) as
where the coefficients ℎ i can be found in Table 2.
Table 2. Coefficients ℎ i for Equation 4.
1.2 Reaction Kinetics
When considering engineering mathematical models, we need to pay close attention to time scales. This means that some parts of our problem might be better represented as stationary phenomena, like in 1.1, but other parts of happen faster, and it is important to understand how they develop in time.
We will now explore a time-dependent model of the solubility of CO2 in water. In specific, we will focus on mathematical models of the cascade of chemical reactions triggered when CO2 is dissolved in water.
Carbon dioxide does not stay inert when it dissolves in water. It undergoes a chain of reversible chemical reactions producing carbonic acid (H2 CO3 ), bicarbonate ion (HCO3(−)), carbonate ion (CO3(2) − ) and hydrogen ion (H+ ).
This chain is expressed through Reactions 1 – 4 as:
The kinetics of these reactions can be modelled by a system of nonlinear ordinary differential equations given by:
In these equations, t is time in seconds and a square bracket around the name of a chemical species indicates its concentration. ki represents rate constants for the ith reaction. Since these reactions are reversible, ki denotes a forward rate and k −i a reverse rate.
The rate constants ki and k −i are given in inverse second [s-1]. The set of equations 5.1 – 5.6 can be interpreted as follows:
• Eq. 5.1 models R1 and accounts for the time it takes for carbon dioxide gas, represented as CO2 (g) to dissolve into water and become aqueous CO2 (aq).
• Eq. 5.2 models both R1 and R2. It accounts for:
o Production of CO2 (aq) by the dissolution of CO2 (g) at rate k1 ,
o Loss of CO2 (aq) by reverse reactions back to CO2 (g) with rate k −1 ,
o Loss/production of CO2 (aq) by forward/reverse reactions to/from H2 CO3 .
A similar logic can be used for Eqs. 5.3 – 5.6. The concentrations in Eqs. 5.1 – 5.6 are given in molarity [M] which is equivalent to one mol of the chemical species per litre of water. The best case-study to understand the dynamics of CO2 dissolution in water is to imagine a system where all species are in equilibrium at t = 0. Table 3 shows the equilibrium molarities for all species in this reaction.
Table 3. Forward and reverse rate constants.
If at t = 0 we inject gaseous carbon dioxide at concentration [CO2 (g)]0 = 0.065 [M], this will trigger a cascade of reactions that will cause all other species to react and then reach equilibrium within 100 seconds.
The equilibrium conditions in [M] for all species which should be used as initial conditions in solving this system of differential equations at 25˚C and 1 atm pressure are:
[CO2 (aq)]0 = 5.41 × 10 −4
[H2 CO3 ]0 = 1.64 × 10 −6
[HCO3(−)]0 = 3.28 × 10 −4
[CO3(2) − ]0 = 1.97 × 10 −8
[H+]0 = 1 × 10 −6
Implementation tips:
• Use the self-paced course Solving Ordinary Differential Equations with MATLAB to learn how to write solutions to systems of non-linear ODEs in MATLAB.
• Since the scales of the rate parameters in Table 3 vary widely from 10 −2 to 1010 , the best MATLAB solvers for this system are stiff solvers such as ode15s or ode23s.
• You are also recommended to set an Absolute Tolerance ‘AbsTol’ of 10 −12 and Relative Tolerance ‘RelTol’ of 10 −6 for this problem because the initial conditions are as low as 10 −8 and might get smaller.
[40 marks] Task 1
Your task is to study, describe, and operationalise the mathematical model in section 1.1 of this document. [3 PAGE MAXIMUM]
A. Mathematical Task: Express Eqs. 1.1, 1.2, 2.1, and 2.2 as a linear system in terms of matrix-vector multiplication and matrix inversion. Find the conditions where the model represented by Eqs. 1.1, 1.2, 2.1, and 2.2 cannot be solved. Discuss the implications of this condition in a real- life system.
B. Communication Task: Create a schematic showing how all equations in Appendix A and pages 4 through to 8 in this brief are connected. Your summary should communicate clearly how they can be solved in a logical order. Use this schematic to describe how you will structure MATLAB code to solve this problem.
C. Modelling Task: Propose an extended version of the model in Equations 1.1 through to 2.2 that also includes gases other than CO2 . Describe the specific quantities that you would need before implementing this model computationally.
[30 marks] Task 2
Your task is to implement and validate the mathematical model in 1.1 in MATLAB. [2 PAGE MAXIMUM]
A. Coding Task: Use MATLAB to calculate and plot H21 , Pv,1 , φ1 and φ2 when T ∈ [10, 80] ˚C and P = 101.325 kPa. Use the subplot function to produce a 2x2 grid of figures.
B. Validation Task: Use MATLAB to calculate the equilibrium concentration of CO2 in water x2 (T, P) for T ∈ [10, 80] ˚C and four distinct pressures such that P ∈ {50, 101.325, 200} kPa. Contrast and compare your results with those of Table 3 in (Carroll, Slupsky and Mather, 1991) .
C. Modelling Task: Attached to this brief is a comprehensive dataset of the Great Barrier Reef. Use data on pressure, temperature and air concentration of CO2 to model the equilibrium concentration of CO2 in water over time.
[30 marks] Task 3
Your task is to study, describe, and operationalise the mathematical model in section 1.2 of this document. [2 PAGE MAXIMUM]
A. Coding Task: Use MATLAB to solve the system of ordinary differential equations shown in Eqs. 5.1 – 5.6 using the reaction rates given in Table 3 and the initial conditions outlined in page 11. Plot your time-dependent solutions for all chemical species in a single figure.
B. Analysis Task: Contrast and compare the speed at which each of the reactions in page 9 takes place. Use the constants in Table 3 and your numerical solution of the system 5.1 – 5.6 to support your argument. Estimate the time it takes for the full system to reach steady-state equilibrium.
C. Summary Task: Contrast and compare the Phase Equilibrium and Reaction Kinetics models in this project. Explain what phenomena they represent and how these are connected. Discuss their fundamental assumptions and how these assumptions shape and limit their possibilities.