function[pco2,pH,co2,hco3,co3] = co3eq(temp,s,z,alk,dic) % Function to calculate pCO2 (microatm), bicarbonate ion (micromol/kg), and % carbonate ion (micromol/kg) concentrations from Alkalinity (micromol/kg) and % DIC (micromol/kg). This function requires temperature (degrees C), % salinity (ppt), and depth (m), and uses the equations in Zeebe and % Wolf-Gladrow (2000). t = temp + 273.15; Pr = z/10; alk = alk*.000001; dic = dic*.000001; R = 83.131; % Calculate total borate from chlorinity tbor = .000416 * s / 35.0; % Calculate Henry's Law coeffieicent (Weiss, 1974) U1 = -60.2409 + 93.4517 * (100/t) + 23.3585*log(t/100); % note that *log* in matlab is the natural log, ln U2 = s * (.023517 - .023656 * (t/100) + .0047036 * (t/100)^2); KH = exp(U1+U2); % Calculate KB from temp and salinity (Dickson, 1990) KB = exp((-8966.9 - 2890.53 * s^0.5 - 77.942 * s + 1.728 * s^1.5 - 0.0996 * s^2)/t ... + 148.0248 + 137.1942 * s^0.5 + 1.62142 * s - (24.4344 + 25.085 * s^0.5 + ... 0.2474 * s) * log(t) + 0.053105 * s^0.5 * t); % Calculate K1 and K2 (Luecker et al., 2000) K1 = 10^(-(3633.86/t - 61.2172 + 9.67770 * log(t) - 0.011555 * s + 0.0001152 * s^2)); K2 = 10^(-(471.78/t + 25.92990 - 3.16967 * log(t) - 0.01781*s + 0.0001122 * s^2)); % Pressure variation of K1, K2, and KB (Millero, 1995) dvB = -29.48 + 0.1622 * temp - .002608 * (temp)^2; dv1 = -25.50 + 0.1271 * temp; dv2 = -15.82 - 0.0219 * temp; dkB = -.00284; dk1 = -.00308 + 0.0000877 * temp; dk2 = .00113 - .0001475 * temp; KB = (exp(-(dvB / (R * t)) * Pr + (0.5 * dkB / (R * t)) * Pr^2)) * KB; K1 = (exp(-(dv1 / (R * t)) * Pr + (0.5 * dk1 / (R * t)) * Pr^2)) * K1; K2 = (exp(-(dv2 / (R * t)) * Pr + (0.5 * dk2 / (R * t)) * Pr^2)) * K2; % temperature dependence of KW (DoE, 1994) KW1 = 148.96502 - 13847.26 / t - 23.65218 * log(t); KW2 = (118.67 / t - 5.977 + 1.0495 * log(t)) * s^.5 - 0.01615 * s; KW = exp(KW1 + KW2); % solve for H ion (Zeebe and Wolf-Gladrow, 2000) a1 = 1; a2 = (alk + KB + K1); a3 = (alk * KB - KB * tbor - KW + alk * K1 + K1 * KB + K1 * K2 - dic * K1); a4 = (-KW * KB + alk * KB * K1 - KB * tbor * K1 - KW * K1 + alk * K1 * K2 ... + KB * K1 * K2 - dic * KB * K1 - 2 * dic * K1 * K2); a5 = (-KW * KB * K1 + alk * KB * K1 * K2 - KW * K1 * K2 - KB * tbor * K1 * K2 ... - 2 * dic * KB * K1 * K2); a6 = -KB * KW * K1 * K2; p = [a1 a2 a3 a4 a5 a6]; r = roots(p); h = max(real(r)); % calculate bicarbonate, carbonate, and aqueous CO2 using DIC, Alk, and H+ format short g; hco3 = dic/ (1 + h/K1 + K2/h) * 1e6; co3 = dic / (1 + h/K2 + h * h / (K1 * K2)) * 1e6; co2 = dic / (1 + K1/h + K1 * K2 / (h * h)) * 1e6; pco2 = co2 / KH ; pH = -log10(h); % calculate B(OH)4 and OH BOH4 = KB * tbor / (h + KB); OH = KW / h; % recalculate DIC and Alk to check calculations Ct = (hco3 + co3 + co2) * 1e6; At = (hco3 + 2*co3 + BOH4 + OH - h) * 1e6;