% _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ % % file: csys.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: given two components of CO2 system, calculate all other. % % Compare book by % % Zeebe & Wolf-Gladrow, CO2 in Seawater: % Equilibrium, Kinetics, Isotopes. % updates: % % % 03.01.01 k -> K % 01.10.98 include pressure % 08.05.95 set phflag Total or free % % MATLAB 5.3.1 % % remarks: specify S, T, and P, flag-number for the two variables % pH flag, choose constants (Roy or Mehrbach), % and values for two variables. % % the variable s is equivalent to [CO2] % % use: save file in directory. open matlab. type 'csys' % % --------------------------------------------------------------- %============= INPUT SECTION ===============================% % flag = 1 h and CO2 given % flag = 101 h and pco2 given % flag = 2 CO2 and HCO3 given % flag = 3 CO2 and CO3 given % flag = 4 CO2 and ALK given % flag = 5 CO2 and DIC given % flag = 6 h and HCO3 given % flag = 7 h and CO3 given % flag = 8 h and ALK given % flag = 9 h and DIC given % flag = 10 HCO3 and CO3 given % flag = 11 HCO3 and ALK given % flag = 12 HCO3 and DIC given % flag = 13 CO3 and ALK given % flag = 14 CO3 and DIC given % flag = 15 ALK and DIC given %============= Set flag ======================================% flag = 101; %-------------- pH Total or pH free ------------% % % Specify on which pH scale we are working on. % DIC, ALK, pCO2, CO2, CO3, and HCO3 are % physical quantities which do not depend on the pH % scale! However, the pH does depend on the scale % used. A certain set of equilibrium % constants belongs to a certain pH scale. % The scale conversion of the constants is % done in equic.m. % % Note that from [H+]_total > [H+]_free % follows pH_total < pH_free % % The total scale is recommended. % phflag = 0; % 0: Total scale % 1: Free scale %----- choose K1 and K2: Roy or Mehrbach. k1k2flag = 1; % 0: Roy et al. (1993) % 1: Mehrbach et al (1973) as % refit by Lueker et al. (2000) on total scale. %----- set temperature, salinity and pressure TC = 25.; % deg C S = 35.; P = 0.; % bar %---- load constants equic; bor = 1.*(416.*(S/35.))* 1.e-6; % (mol/kg), DOE94 %---- set values for carbonate system %---- specify two variables according to flag (see above) ph1 = 7.9; % pH h1 = 10^(-ph1); % [H+] %s1 = 18.e-6; % [CO2] (mol/kg) %hco31 = 1750.e-6; % [HCO3-] (mol/kg) %co31 = 225.e-6; % [CO3--] (mol/kg) %alk1 = 3800.e-6; % ALK (mol/kg) %dic1 = 3900.e-6; % [DIC] (mol/kg) pco21 = 1300.e-6; % pCO2 (atm) %----------- These values are test values, set test= 0/1 -------% test = 0; if test == 1; TC = 25.; S = 35.; P = 0.; equic; bor = (416.*(S/35.))* 1.e-6; % (mol/kg), DOE94 ph1 = 8.0902; % pH (8.0902 Total, 8.1979 Free) h1 = 10^(-ph1); % [H+] s1 = 10.1304e-6; % [CO2] hco31 = 1735.9e-06; % [HCO3-] co31 = 253.9924e-6; % [CO3--] dic1 = 2000.e-6; % [DIC] alk1 = 2350.e-6; % Alk pco21 = 356.8058e-6; % pCO2 end; %============= END INPUT SECTION ===============================% % % rest of file contains the numerical routines. %==================================================================% % ----------------- case 1.) h and s given if flag == 1 disp('flag = 1, pH and CO2 given'); h = h1; s = s1; dic = s*(1.+K1/h1+K1*K2/h1/h1); hco3 = dic/(1+h1/K1+K2/h1); co3 = dic/(1+h1/K2+h1*h1/K1/K2); alk = s*(K1/h1+2.*K1*K2/h1/h1)+Kb*bor/(Kb+h1)+Kw/h1-h1; pco2 = s/Kh; % ----------- change units: mumol/kg CO2 = 1.e6*s1 HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 pco2 = pco2*1.e6 %alktest = (2*co31 + hco31 + bor/(1+h1/Kb) + Kw/h1 - h1)*1.e6 end; % ----------------- case 101.) h and pco2 given if flag == 101 disp('flag = 101, pH and pCO2 given'); h = h1; pco2 = pco21; s = Kh*pco2; dic = s*(1.+K1/h1+K1*K2/h1/h1); hco3 = dic/(1+h1/K1+K2/h1); co3 = dic/(1+h1/K2+h1*h1/K1/K2); alk = s*(K1/h1+2.*K1*K2/h1/h1)+Kb*bor/(Kb+h1)+Kw/h1-h1; % ----------- change units: mumol/kg CO2 = 1.e6*s HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 pco2 = pco2*1.e6 %alktest = (2*co31 + hco31 + bor/(1+h1/Kb) + Kw/h1 - h1)*1.e6 end; % ------------ s and HCO3 given ------------------ if flag == 2 disp('flag = 2, CO2 and HCO3 given'); s = s1; hco3 = hco31; p3 = -hco3/K1; p2 = s - hco3; p1 = s*K1 - hco3*K2; p0 = s*K1*K2; p = [p3 p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; dic = s*(1.+K1/h+K1*K2/h/h); % hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ s and CO3 given ------------------ if flag == 3 disp('flag = 3, CO2 and CO3 given'); s = s1; co3 = co31; p4 = -co3/K1/K2; p3 = -co3/K2; p2 = s-co3; p1 = s*K1; p0 = s*K1*K2; p = [p4 p3 p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); % co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ s and ALK given ------------------ if flag == 4 disp('flag = 4, CO2 and ALK given'); s = s1; alk = alk1; p4 = 1.; p3 = Kb+alk; p2 = alk*Kb-s*K1-Kb*bor-Kw; p1 = -s*Kb*K1-s*2.*K1*K2-Kw*Kb; p0 = -2.*s*Kb*K1*K2; p = [p4 p3 p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); % alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ s and DIC given ------------------ if flag == 5 disp('flag = 5, CO2 and DIC given'); s = s1; dic = dic1; p2 = dic - s; p1 = -s*K1; p0 = -s*K1*K2; p = [p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; % dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ h and HCO3 given ------------------ if flag == 6 disp('flag = 6, pH and HCO3 given'); h = h1; hco3 = hco31; dic = hco3 * (1+h/K1+K2/h); s = dic / (1.+K1/h+K1*K2/h/h); h*1.e12; % dic = s*(1.+K1/h+K1*K2/h/h); % hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ h and CO3 given ------------------ if flag == 7 disp('flag = 7, pH and CO3 given'); h = h1; co3 = co31; dic = co3 * (1+h/K2+h*h/K1/K2); s = dic / (1.+K1/h+K1*K2/h/h); h*1.e12; % dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); % co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ h and ALK given ------------------ if flag == 8 disp('flag = 8, pH and ALK given'); h = h1; alk = alk1; s = (alk-Kw/h+h-Kb*bor/(Kb+h)) / (K1/h+2.*K1*K2/h/h); h*1.e12; dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); % alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ h and DIC given ------------------ if flag == 9 disp('flag = 9, pH and DIC given'); h = h1; dic = dic1; s = dic / (1.+K1/h+K1*K2/h/h); h*1.e12; % dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ HCO3 and CO3 given ------------------ if flag == 10 disp('flag = 10, HCO3 and CO3 given'); hco3 = hco31; co3 = co31; p3 = -co3/K1/K2; p2 = -co3/K2 + hco3/K1; p1 = -co3 + hco3; p0 = hco3*K2; p = [p3 p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; dic = hco3 * (1+h/K1+K2/h); s = dic / (1.+K1/h+K1*K2/h/h); % hco3 = dic/(1+h/K1+K2/h); % co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ HCO3 and ALK given ------------------ % % don't split the lines, % matlab does not understand it. % if flag == 11 disp('flag = 11, HCO3 and ALK given'); hco3 = hco31; alk = alk1; p5 = 1.; p4 = alk - hco3 + K1 + Kb; p3 = alk*(Kb+K1)-hco3*(K1+Kb+2.*K2)-Kw+K1*Kb+K1*K2-Kb*bor; tmp = alk*(Kb*K1+K1*K2)-hco3*((Kb+2.*K2)*K1+2.*Kb*K2+K1*K2); p2 = tmp +(-K1*Kb*bor-Kw*Kb-K1*Kw+K1*K2*Kb); tmp = alk*Kb*K1*K2-hco3*(2.*Kb*K1*K2+K2*K1*(Kb+2.*K2)); p1 = tmp +(-K1*K2*Kb*bor-K1*Kw*Kb-K1*K2*Kw); p0 = -hco3*2.*K2*Kb*K1*K2-K1*K2*Kw*Kb; p = [p5 p4 p3 p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; dic = hco3 * (1+h/K1+K2/h); s = dic / (1.+K1/h+K1*K2/h/h); % dic = s*(1.+K1/h+K1*K2/h/h); % hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ HCO3 and DIC given ------------------ if flag == 12 disp('flag = 12, HCO3 and DIC given'); hco3 = hco31; dic = dic1; p2 = hco3/K1; p1 = hco3-dic; p0 = hco3*K2; p = [p2 p1 p0]; r = roots(p); h = min(real(r)); % min instead of max !!!!! h*1.e12; s = dic / (1.+K1/h+K1*K2/h/h); dic = s*(1.+K1/h+K1*K2/h/h); % hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ CO3 and ALK given ------------------ if flag == 13 disp('flag = 13, CO3 and ALK given'); co3 = co31; alk = alk1; p5 = -co3/K2+1.; p4 = alk - co3*(K1/K2+(Kb+2.*K2)/K2) + Kb + K1; tmp = alk*(Kb+K1)-co3*(K1+K1*(Kb+2.*K2)/K2+2.*Kb); p3 = tmp+(-Kb*bor-Kw+K1*Kb+K1*K2); tmp = alk*(Kb*K1+K1*K2)-co3*(K1*(Kb+2.*K2)+2.*Kb*K1); p2 = tmp+(-Kw*Kb-K1*Kb*bor-K1*Kw+K1*K2*Kb); tmp = alk*Kb*K1*K2-co3*2.*Kb*K1*K2-K1*Kw*Kb; p1 = tmp+(-K1*K2*Kb*bor-K1*K2*Kw); p0 = -K1*K2*Kw*Kb; p = [p5 p4 p3 p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; % % NOTE : calculate dic from dic = co3*(1+h/K2+h^2/K1/K2); % not from dic = hco3*(1+h/K1+K2/h); % dic = co3 * (1+h/K2+h^2/K1/K2); s = dic / (1.+K1/h+K1*K2/h/h); % dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); % co3 = dic/(1+h/K2+h*h/K1/K2); % alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ CO3 and DIC given ------------------ if flag == 14 disp('flag = 14, CO3 and DIC given'); co3 = co31; dic = dic1; p2 = co3/K1/K2; p1 = co3/K2; p0 = co3-dic; p = [p2 p1 p0]; r = roots(p); h = max(real(r)); h*1.e12; s = dic / (1.+K1/h+K1*K2/h/h); % dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); % co3 = dic/(1+h/K2+h*h/K1/K2); alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); % 6.3096e-09 CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end % ------------ ALK and DIC given ------------------ if flag == 15 disp('flag = 15, ALK and DIC given'); alk = alk1; dic = dic1; p5 = -1.; p4 = -alk-Kb-K1; p3 = dic*K1-alk*(Kb+K1)+Kb*bor+Kw-Kb*K1-K1*K2; tmp = dic*(Kb*K1+2.*K1*K2)-alk*(Kb*K1+K1*K2)+Kb*bor*K1; p2 = tmp+(Kw*Kb+Kw*K1-Kb*K1*K2); tmp = 2.*dic*Kb*K1*K2-alk*Kb*K1*K2+Kb*bor*K1*K2; p1 = tmp+(+Kw*Kb*K1+Kw*K1*K2); p0 = Kw*Kb*K1*K2; p = [p5 p4 p3 p2 p1 p0]; r = roots(p); h = max(real(r)); % test = p5*h^5+p4*h^4+p3*h^3+p2*h^2+p1*h+p0 h*1.e12; s = dic / (1.+K1/h+K1*K2/h/h); % dic = s*(1.+K1/h+K1*K2/h/h); hco3 = dic/(1+h/K1+K2/h); co3 = dic/(1+h/K2+h*h/K1/K2); % alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h; % ----------- change units: mumol/kg %h1 = 10^(-ph1); % CO2 = s*1.e6 pCO2 = s*1.e6/Kh HCO3 = hco3*1.e6 CO3 = co3*1.e6 DIC = dic*1.e6 ALK = alk*1.e6 end ; % ====================================================== if phflag == 0; disp('Total pH scale used.'); phtotal = -log10(h) oh = Kw./h; phfree = -log10(h/total2free) end; if phflag == 1; disp('Free pH scale used.'); phfree = -log10(h) oh = Kw./h; phtotal = -log10(h*total2free) end; return;