% _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ % % file: equic.m % % % Richard E. Zeebe & Dieter A. Wolf-Gladrow % % % Alfred Wegener Institute for % Polar and Marine Research % P.O. Box 12 01 61 % D-27515 Bremerhaven % Germany % e-mail: rzeebe@awi-bremerhaven.de wolf@awi-bremerhaven.de % www : http://www.awi-bremerhaven.de/ % % _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ % % purpose: CO2 equilibrium constants. Compare book by % % Zeebe & Wolf-Gladrow, CO2 in Seawater: % Equilibrium, Kinetics, Isotopes. % % % updates: % 15.02.01 density of seawater added % 04.01.01 K1- and K2-Mehrbach as refit by Lueker added % 04.01.01 phosphoric acid + pressure term added % 03.01.01 k -> K % 30.10.00 pressure eff. on Ksp* calcite and aragonite % 06.12.99 Kf and pressure eff. on Kf % 01.12.99 K_(calcite), K_(aragonite) % 01.10.98 include pressure effect (Millero, 1995) % 07.05.96 use constants summarized in Dickson and Goyet, 1994 % % % % MATLAB 5.3.1 % remarks: % units: mol/kg-soln % set phflag (0:Total), (1:Free) % =================================================================== tk = 273.15; % [K] (for conversion [deg C] <-> [K]) T = TC + tk; % TC [C]; T[K] Cl = S / 1.80655; % Cl = chlorinity; S = salinity (per mille) cl3 = Cl.^(1/3); ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl .* Cl; % ionic strength iom0 = 19.924*S/(1000.-1.005*S); S_T = 0.14/96.062/1.80655*S; % (mol/kg soln) total sulfate % Dickson and Goyet (1994) Ch.5 p.11 %------------------ Ks ---------------------------------------- % Dickson and Goyet (1994), Chapter 5, p.13 % (required for total2free) % Equilibrium constant for HSO4- = H+ + SO4-- % % K_S = [H+]free [SO4--] / [HSO4-] % pH-scale: free scale !!! % % the term log(1-0.001005*S) converts from % mol/kg H2O to mol/kg soln tmp1 = -4276.1 ./ T + 141.328 -23.093*log(T); tmp2 = +(-13856 ./ T + 324.57 - 47.986 * log(T)).*sqrt(iom0); tmp3 = +(35474 ./ T - 771.54 + 114.723 * log(T)).*iom0; tmp4 = -2698 ./ T .*sqrt(iom0).*iom0 + 1776 ./ T .*iom0 .*iom0; lnKs = tmp1 + tmp2 + tmp3 + tmp4 + log(1-0.001005*S); Ks = exp(lnKs); %------- total2free ----------------------------------------------- % % convert from pH_total ('total`) to pH ('free`): % pH_total = pH_free - log(1+S_T/K_S(s,tk)) total2free = 1.+S_T./Ks; % --------------------- Kf -------------------------------------------- % Kf = [H+][F-]/[HF] % % (Dickson and Riley, 1979 in Dickson and Goyet, % 1994, Chapter 5, p. 14) % pH-scale: 'total' tmp1 = 1590.2./T - 12.641 + 1.525.*sqrt(ION); tmp2 = log(1.-0.001005.*S) + log(1.+S_T./Ks); lnKf = tmp1 + tmp2; Kf = 0; if phflag == 0; Kf = exp(lnKf); end; if phflag == 1; lnKf = lnKf-log(total2free); Kf = exp(lnKf); end; %------- sws2free ----------------------------------------------- % % convert from pH_sws ('seawater scale`) to pH ('free`): % pH_sws = pH_free - log(1+S_T/K_S(S,T)+F_T/K_F(S,T)) F_T = 7.e-5.*(S./35.); sws2free = (1.+S_T./Ks+F_T./Kf); corr = sws2free./total2free; % --------------------- Kwater ----------------------------------- % % Millero (1995)(in Dickson and Goyet (1994, Chapter 5, p.18)) % $K_w$ in mol/kg-soln. % pH-scale: pH$_{total}$ ('total` scale). tmp1 = -13847.26./T + 148.96502 - 23.6521 .* log(T); tmp2 = + (118.67./T - 5.977 + 1.0495.*log(T)).*sqrt(S) - 0.01615.*S; lnKw = tmp1 + tmp2; if phflag == 0; Kw = exp(lnKw); end; if phflag == 1; lnKw = lnKw-log(total2free); Kw = exp(lnKw); end; %---------------------- Kh (K Henry) ---------------------------- % % CO2(g) <-> CO2(aq.) % Kh = [CO2]/ p CO2 % % Weiss (1974) [mol/kg/atm] % % % tmp = 9345.17 ./ T - 60.2409 + 23.3585 * log(T/100.); nKhwe74 = tmp + S.*(0.023517-0.00023656*T+0.0047036e-4*T.*T); tmpKh1= 9050.69 ./ T - 58.0931 + 22.2940 * log(T/100.); nKhwe74l = tmpKh1 + S .* (0.027766-0.00025888*T+0.0050578e-4 * T .* T); %Kh= exp(nKhwe74l); Kh= exp(nKhwe74); % --------------------- K1 --------------------------------------- % first acidity constant: % [H^+] [HCO_3^-] / [CO2] = K_1 % % (Roy et al., 1993 in Dickson and Goyet, 1994, Chapter 5, p. 14) % pH-scale: 'total'. mol/kg-soln tmp1 = 2.83655 - 2307.1266 ./ T - 1.5529413 .* log(T); tmp2 = - (0.20760841 + 4.0484 ./ T) .* sqrt(S); tmp3 = + 0.08468345 .* S - 0.00654208 .* S .* sqrt(S); tmp4 = + log(1 - 0.001005 .* S); lnK1roy = tmp1 + tmp2 + tmp3 + tmp4; K1roy = 0; if phflag == 0; K1roy = exp(lnK1roy); end; if phflag == 1; lnK1roy = lnK1roy-log(total2free); K1roy = exp(lnK1roy); end; % --------------------- K2 ---------------------------------------- % % second acidity constant: % [H^+] [CO_3^--] / [HCO_3^-] = K_2 % % (Roy et al., 1993 in Dickson and Goyet, 1994, Chapter 5, p. 15) % pH-scale: 'total'. mol/kg-soln tmp1 = -9.226508 - 3351.6106 ./ T - 0.2005743 .* log(T); tmp2 = (-0.106901773 - 23.9722 ./ T) .* sqrt(S); tmp3 = 0.1130822 .* S - 0.00846934 .* S.^1.5 + log(1 - 0.001005 * S); lnK2roy = tmp1 + tmp2 + tmp3; if phflag == 0; K2roy = exp(lnK2roy); end; if phflag == 1; lnK2roy = lnK2roy-log(total2free); K2roy = exp(lnK2roy); end; % --------------------- K1 --------------------------------------- % first acidity constant: % [H^+] [HCO_3^-] / [H_2CO_3] = K_1 % % Mehrbach et al (1973) refit by Lueker et al. (2000). % % pH-scale: 'total'. mol/kg-soln pK1mehr = 3633.86./T - 61.2172 + 9.6777.*log(T) - 0.011555.*S + 0.0001152.*S.*S; if phflag == 0; K1mehr = 10^(-pK1mehr); end; if phflag == 1; lnK1mehr = log(10^(-pK1mehr))-log(total2free); K1mehr = exp(lnK1mehr); end; % --------------------- K2 ---------------------------------------- % % second acidity constant: % [H^+] [CO_3^--] / [HCO_3^-] = K_2 % % Mehrbach et al. (1973) refit by Lueker et al. (2000). % % pH-scale: 'total'. mol/kg-soln pK2mehr = 471.78./T + 25.9290 - 3.16967.*log(T) - 0.01781.*S + 0.0001122.*S.*S; if phflag == 0; K2mehr = 10^(-pK2mehr); end; if phflag == 1; lnK2mehr = log(10^(-pK2mehr))-log(total2free); K2mehr = exp(lnK2mehr); end; %----------- Roy or Mehrbach. default: Roy K1 = K1roy; K2 = K2roy; if (exist('k1k2flag')) if k1k2flag == 0; K1 = K1roy; K2 = K2roy; end; if k1k2flag == 1; K1 = K1mehr; K2 = K2mehr; end; end; % --------------------- Kb -------------------------------------------- % Kbor = [H+][B(OH)4-]/[B(OH)3] % % (Dickson, 1990 in Dickson and Goyet, 1994, Chapter 5, p. 14) % pH-scale: 'total'. mol/kg-soln tmp1 = (-8966.90-2890.53*sqrt(S)-77.942*S+1.728*S.^(3./2.)-0.0996*S.*S); tmp2 = +148.0248+137.1942*sqrt(S)+1.62142*S; tmp3 = +(-24.4344-25.085*sqrt(S)-0.2474*S).*log(T); lnKb = tmp1 ./ T + tmp2 + tmp3 + 0.053105*sqrt(S).*T; if phflag == 0; Kb = exp(lnKb); end; if phflag == 1; lnKb = lnKb-log(total2free); Kb = exp(lnKb); end; % --------------------- Phosphoric acid --------------------- % % % (DOE, 1994) (Dickson and Goyet): pH_T, mol/(kg-soln) % Ch.5 p. 16 % % lnK1P = -4576.752 ./ T + 115.525 - 18.453*log(T) ... + (-106.736 ./ T + 0.69171) .* sqrt(S) ... + (-0.65643 ./ T - 0.01844) .* S; lnK2P = -8814.715 ./ T + 172.0883 - 27.927 * log(T) ... + (-160.34 ./ T + 1.3566) .* sqrt(S) ... + (0.37335 ./ T - 0.05778) .* S; lnK3P = -3070.75 ./ T - 18.141 ... + (17.27039 ./ T + 2.81197) .* sqrt(S) ... + (-44.99486 ./ T - 0.09984) .* S; K1P = exp(lnK1P); K2P = exp(lnK2P); K3P = exp(lnK3P); % --------------------- Silicic acid --------------------------- % % (DOE, 1994) (Dickson and Goyet): pH_T, mol/(kg-soln) % Ch.5 p. 17 % % lnKSi = -8904.2 ./ T + 117.385 - 19.334*log(T) ... + (3.5913-458.79 ./ T) .* sqrt(iom0) + (188.74 ./ T - 1.5998) .* iom0 ... + (0.07871 - 12.1652 ./ T) .*iom0.^2 + log(1-0.001005*S); KSi = exp(lnKSi); % --------------------- Kspc (calcite) ---------------------------- % % apparent solubility product of calcite % % Kspc = [Ca2+]T [CO32-]T % % where $[]_T$ refers to the equilibrium total % (free + complexed) ion concentration. % % Mucci 1983 mol/kg-soln tmp1 = -171.9065-0.077993.*T+2839.319./T+71.595.*log10(T); tmp2 = +(-0.77712+0.0028426.*T+178.34./T).*sqrt(S); tmp3 = -0.07711.*S+0.0041249.*S.^1.5; log10Kspc = tmp1 + tmp2 + tmp3; Kspc = 10.^(log10Kspc); % --------------------- Kspa (aragonite) ---------------------------- % % apparent solubility product of aragonite % % Kspa = [Ca2+]T [CO32-]T % % where $[]_T$ refers to the equilibrium total % (free + complexed) ion concentration. % % Mucci 1983 mol/kg-soln tmp1 = -171.945-0.077993.*T+2903.293./T+71.595.*log10(T); tmp2 = +(-0.068393+0.0017276.*T+88.135./T).*sqrt(S); tmp3 = -0.10018.*S+0.0059415.*S.^1.5; log10Kspa = tmp1 + tmp2 + tmp3; Kspa = 10.^(log10Kspa); %---------------------------------------------------- % % Density of seawater as function of S,T,P. % % Millero et al. 1981, Gill, 1982. % % % % %---------------------------------------------------- %------------ Density of pure water rhow = 999.842594 + 6.793952e-2*TC -9.095290e-3*TC^2 ... + 1.001685e-4*TC^3 -1.120083e-6*TC^4 + 6.536332e-9*TC^5; %------------ Density of seawater at 1 atm, P=0 A = 8.24493e-1 - 4.0899e-3*TC + 7.6438e-5*TC^2 - 8.2467e-7*TC^3 ... + 5.3875e-9*TC^4; B = -5.72466e-3 + 1.0227e-4*TC - 1.6546e-6*TC^2; C = 4.8314e-4; rho0 = rhow + A*S + B*S^(3/2) + C*S^2; %-------------- Secant bulk modulus of pure water % % The secant bulk modulus is the average change in pressure % divided by the total change in volume per unit of initial volume. Ksbmw = 19652.21 + 148.4206*TC - 2.327105*TC^2 ... + 1.360477e-2*TC^3 - 5.155288e-5*TC^4; %-------------- Secant bulk modulus of seawater at 1 atm Ksbm0 = Ksbmw ... + S*( 54.6746 - 0.603459*TC + 1.09987e-2*TC^2 ... - 6.1670e-5*TC^3) ... + S^(3/2)*( 7.944e-2 + 1.6483e-2*TC - 5.3009e-4*TC^2); %-------------- Secant bulk modulus of seawater at S,T,P Ksbm = Ksbm0 ... + P*( 3.239908 + 1.43713e-3*TC + 1.16092e-4*TC^2 ... - 5.77905e-7*TC^3) ... + P*S*( 2.2838e-3 - 1.0981e-5*TC - 1.6078e-6*TC^2) ... + P*S^(3/2)*1.91075e-4 ... + P*P*(8.50935e-5 - 6.12293e-6*TC + 5.2787e-8*TC^2) ... + P^2*S*(-9.9348e-7 + 2.0816e-8*TC + 9.1697e-10*TC^2); %------------- Density of seawater at S,T,P rho = rho0/(1.-P/Ksbm); %---------------------- Pressure effect on K's (Millero, 95) ----------% if P > 0.0 RGAS = 8.314510; % J mol-1 deg-1 (perfect Gas) R = 83.131; % mol bar deg-1 % conversion cm3 -> m3 *1.e-6 % bar -> Pa = N m-2 *1.e+5 % => *1.e-1 or *1/10 % index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8, % K1P 9, K2P 10, K3P 11 %----- note: there is an error in Table 9 of Millero, 1995. %----- The coefficients -b0 and b1 %----- have to be multiplied by 1.e-3! %----- there are some more errors! %----- the signs (+,-) of coefficients in Millero 95 do not %----- agree with Millero 79 a0 = -[25.5 15.82 29.48 25.60 18.03 9.78 48.76 46. ... 14.51 23.12 26.57]; a1 = [0.1271 -0.0219 0.1622 0.2324 0.0466 -0.0090 0.5304 0.5304 ... 0.1211 0.1758 0.2020]; a2 = [0.0 0.0 2.608 -3.6246 0.316 -0.942 0.0 0.0 ... -0.321 -2.647 -3.042]*1.e-3; b0 = -[3.08 -1.13 2.84 5.13 4.53 3.91 11.76 11.76 ... 2.67 5.15 4.08]*1.e-3; b1 = [0.0877 -0.1475 0.0 0.0794 0.09 0.054 0.3692 0.3692 ... 0.0427 0.09 0.0714]*1.e-3; b2 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0]; for ipc=1:length(a0); deltav(ipc) = a0(ipc) + a1(ipc).*TC + a2(ipc).*TC.*TC; deltak(ipc) = (b0(ipc) + b1(ipc).*TC + b2(ipc).*TC.*TC); lnkpok0(ipc) = -(deltav(ipc)./(R.*T)).*P + (0.5*deltak(ipc)./(R.*T)).*P.*P; end; K1 = K1*exp(lnkpok0(1)); K2 = K2*exp(lnkpok0(2)); Kb = Kb*exp(lnkpok0(3)); Kw = Kw*exp(lnkpok0(4)); Ks = Ks*exp(lnkpok0(5)); Kf = Kf*exp(lnkpok0(6)); Kspc = Kspc*exp(lnkpok0(7)); Kspa = Kspa*exp(lnkpok0(8)); K1P = K1P*exp(lnkpok0(9)); K2P = K2P*exp(lnkpok0(10)); K3P = K3P*exp(lnkpok0(11)); end; return;