Contents

clear
clc

% CO2 SOLUBILITY IN PURE WATER AND AQUEOUS SOLUTIONS

% Model written based on "An improved model calculating CO2 solubility in
% pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to
% 2000 bar" by Zhenhao Duan and Rui Sun (2003), and "An improved model for
% the calculation of CO2 solubility in aqueous solutions containing Na+,
% K+, Ca2+, Cl-, and SO42-" by Zhenhao Duan and Rui Sun (2005). Valid range
% T[K]: 273 < T < 533 , P[bar]: 0 < P < 2000, [m/kg]: 0 < m < 4.3

INPUT

check = 1;
% (1/0) 1 = Only calculate value. 0 = Produce surface plot of whole valid interval.


Pmin = 0;     % Minimum pressure [Bar]
Pmax = 2000;  % Maximum pressure [Bar]
Tmin = 273;   % Minimum Temperature [Kelvin]
Tmax = 533;   % Maximum Temperature [Kelvin]

% SSW composition
mNa    = 0.45;      % Molality [mol/kg]
mK     = 0.01;      % Molality [mol/kg]
mCa    = 0.013;     % Molality [mol/kg]
mMg    = 0.045;     % Molality [mol/kg]
mCl    = 0.525;     % Molality [mol/kg]
mSO4   = 0.024;     % Molality [mol/kg]

if check == 1;
    P = 50;     %[bar]
    T = 298.15; %[K]
    solubility = CO2sol(P,T,mNa,mCa,mK,mMg,mCl,mSO4)
    return
end

% Creating solubility matrix of the whole valid interval

Prange = linspace(Pmin,Pmax,(Pmax-Pmin)+1);
Trange = linspace(Tmin,Tmax,(Tmax-Tmin)+1);
for j = 1:(Tmax-Tmin)+1
    T = Trange(j);
    for i = 1:(Pmax-Pmin)+1
        P = Prange(i);
        solubility(i,j) = CO2sol(P,T,mNa,mCa,mK,mMg,mCl,mSO4);
    end
end

solubility(imag(solubility)~=0) = 0;
mesh(Trange,Prange,solubility)

Published with MATLAB® R2014b

Contents

function[CO2insolution] = CO2sol(P,T,mNa,mCa,mK,mMg,mCl,mSO4)

% CHECK IF WITHIN VALID RANGE

if T < 273 || T > 533
       disp('Temperature out of range!')
       return;
end
if P < 0 || P > 2000
       disp('Pressure out of range!')
       return;
end
if mNa < 0 || mNa > 4.3
       disp('Consentration of Na out of range!')
       return;
end
if mCa < 0 || mCa > 4.3
       disp('Consentration of Ca out of range!')
       return;
end
if mK < 0 || mK > 4.3
       disp('Consentration of K out of range!')
       return;
end
if mMg < 0 || mMg > 4.3
       disp('Consentration of Mg out of range!')
       return;
end
if mCl < 0 || mCl > 4.3
       disp('Consentration of Cl out of range!')
       return;
end


% Calulation of fugacity coefficient (Non-iterative). Fitted to the iterative EoS from 2003 paper, but with updated parameters for increased accuracy.

% Calculating test parameter
if T < 305
    P1 = (1.1617*T^2-545.9*T+65929)*10^(-2);
    % CO2 saturation pressure fitted to data from http://www.linde-gas.ro. R^2 = 1

    elseif T >= 305 && T < 405;
        P1 = 75+(T-305)*1.25;
    elseif T >= 405
        P1 = 200;
end

% Fitting region 1
if T >= 273 && T < 573 && P<P1
    c1 = 1; c2 = 4.7586835*10^(-3); c3 = -3.3569963*10^(-6);
    c4 = 0; c5 = -1.3179396; c6 = -3.8389101*10^(-6); c7 = 0;
    c8 = 2.2815104*10^(-3); c9 = 0; c10 = 0; c11 = 0; c12 = 0;
    c13 = 0; c14 = 0; c15 = 0;
end

% Fitting region 2
if T >= 273 && T < 340 && P >= P1 && P < 1000
    c1 = -7.1734882*10^(-1); c2 = 1.5985379*10^(-4);
    c3 = -4.9286471*10^(-7); c4 = 0;  c5 = 0;
    c6 = -2.7855285*10^(-7);  c7 = 1.1877015*10^(-9);
    c8 = 0;  c9 = 0; c10 = 0; c11 = 0; c12 = -96.539512;
    c13 = 4.4774938*10^(-1); c14 = 101.81078;
    c15 = 5.3783879*10^(-6);
end

% Fitting region 3
if T >= 273 && T < 340 && P >= 1000
    c1 = -6.5129019*10^(-2); c2 = -2.1429977*10^(-4);
    c3 = -1.1444930*10^(-6); c4 = 0; c5 = 0; c6 = -1.1558081*10^(-7);
    c7 = 1.1952370*10^(-9); c8 = 0; c9 = 0; c10 = 0; c11 = 0;
    c12 = -221.34306; c13 = 0; c14 = 71.820393; c15 = 6.6089246*10^(-6);
end

% Fitting region 4
if T >= 340 && T < 435 && P >= P1 && P < 1000
    c1 = 5.0383896; c2 = -4.4257744*10^(-3); c3 = 0; c4 = 1.95727333;
    c5 = 0; c6 = 2.4223436*10^(-6); c7 = 0; c8 = -9.3796135*10^(-4);
    c9 = -1.5026030; c10 = 3.0272240*10^(-3); c11 = -31.3777342;
    c12 = -12.847063; c13 = 0; c14 = 0; c15 = -1.5056648*10^(-5);
end

% Fitting region 5
if T >= 340 && T < 435 && P >= 1000
    c1 = -16.063152; c2 = -2.7057990*10^(-3); c3 = 0;
    c4 = 1.4119239*10^(-1); c5 = 0; c6 = 8.1132965*10^(-7);
    c7 = 0; c8 = -1.1453082*10^(-4); c9 = 2.3895671;
    c10 = 5.0527457*10^(-4); c11 = -17.763460; c12 = 985.92232;
    c13 = 0; c14 = 0; c15 = -5.4965256*10^(-7);
end

% Fitting region 6
if T >= 435 && P >= P1
    c1 = -1.5693490*10^(-1); c2 = 4.4621407*10^(-4);
    c3 = -9.1080591*10^(-7); c4 = 0; c5 = 0; c6 = 1.0647399*10^(-7);
    c7 = 2.4273357*10^(-10); c8 = 0; c9 = 3.5874255*10^(-1);
    c10 = 6.3319710*10^(-5); c11 = -249.89661; c12 = 0; c13 = 0;
    c14 = 888.76800; c15 = -6.6348003*10^(-7);
end

% Table for calculating parameters
a1 = (c2+c3*T+c4/T+c5/(T-150))*P;
a2 = (c6+(c7*T)+(c8/T))*P^2;
a3 = (c9+(c10*T)+(c11/T))*log(P);
a4 = (c12+(c13*T))/P;
a5 = c14/T;
a6 = c15*(T^2);

% Fugacity coefficient
phi = c1 + a1 + a2 + a3 + a4 + a5 + a6;


% Emprical calculation of water vapour pressure and mole frac CO2 in gas.

% Assumes water vapour pressure in mixtures are the same as for pure water.

Tcw = 647.29;       % Critical temperature of water [K]
Pcw = 220.85;       % Critical pressure of water [Bar]
t = (T-Tcw)/Tcw;

b1 = -38.640844;
b2 = 5.8948420;
b3 = 59.876516;
b4 = 26.654627;
b5 = 10.637097;

PH2O = ((Pcw*T)/Tcw)*(1+b1*(-t)^1.9 + b2*t + b3*t^2 + b4*t^3 + b5*t^4);
yco2 = (P-PH2O)/P; % mole frac CO2 in vapour.

% Calculating interaction parameters
% Par(T,P) = f1 + f2T + f3/T + f4T^2 + f5/(630-T) + f6P + f7PlnT+ f8P/T + f9P(630-T) + f10P^2/(630-T)^2 + c11TlnP

% Myu-CO2
m1 = 28.9447706; m2 = -0.0354581768; m3 = -4770.67077; m4 = 1.02782768*10^(-5);
m5 = 33.8126098; m6 = 9.04037140*10^(-3); m7 = -1.14934031*10^(-3);
m8 = -0.307405726; m9 = -0.0907301486; m10 = 9.32713393*10^(-4);

% lambdaco2Na
l1 = -0.411370585; l2 = 6.07632013*10^(-4); l3 = 97.5347708; l4 = 0;
l5 = 0; l6 = 0; l7 = 0; l8 = -0.0237622469; l9 = 0.0170656236; l10 = 0;
l11 = 1.41335834*10^(-5);

% zetaco2NaCl
z1 = 3.36389723*10^(-4); z2 = -1.98298980*10^(-5); z3 = 0; z4 = 0;
z5 = 0; z6 = 0; z7 = 0; z8 = 2.12220830*10^-3; z9 = -5.24873303*10^-(3);
z10 = 0; z11 = 0;

%

myuco2RT = m1 + m2*T + m3/T + m4*T^2 + m5/(630-T) + m6*P + m7*P*log(T) + (m8*P)/T + m9*P/(630-T) + (m10*P^2)/((630-T)^2);
lambdaco2Na = l1 + l2*T + l3/T + l4*T^2 + l5/(630-T) + l6*P + l7*P*log(T) + (l8*P)/T + l9*P/(630-T) + (l10*P^2)/((630-T)^2) + l11*T*log(P);
zetaco2NaCl = z1 + z2*T + z3/T + z4*T^2 + z5/(630-T) + z6*P + z7*P*log(T) + (z8*P)/T + z9*P/(630-T) + (z10*P^2)/((630-T)^2) + z11*T*log(P);

Calulating CO2 solved in water [moles/kg]

d1 = -myuco2RT;
d2 = -2*lambdaco2Na*(mNa+mK+2*mCa+2*mMg);
d3 = -zetaco2NaCl*mCl*(mNa+mK+mMg+mCa);
d4 = 0.07*mSO4;
lnmco2 = log(yco2*phi*P)+d1+d2+d3+d4;
CO2insolution = exp(lnmco2);

end

Published with MATLAB® R2014b