Function HCO3_T(Sal As Double, TempC As Double, Pdbar As Double, TP As Double, TSi As Double, TA As Double, TC As Double) As Double ' This code was adapted from the work of others: ' CO2SYS.exe was originally written by Ernie Lewis and Doug Wallace. ' It was converted to Matlab by Richard E. Zeebe & Dieter A. Wolf-Gladrow. ' Parts were extracted and modified by Rachel M. Dunk. ' ' I am indebted to all of them for their initial development of this code. ' ' Function calculates HCO3 ion conc (umole/kg) from TC and TA where ... ' K1 & K2: MEHRBACH refit BY DICKSON AND MILLERO, 1987; ' SO4: Dickson; ' Total pH scale ' 'Dim TempC, Sal, Pdbar, TP, TSi, TA, TC as Double Dim TempK As Double Dim RGasConstant, RT, sqrSal, logTempK, Pbar As Double Dim TB, TF, TS As Double Dim SWStoTOT, pHfactor As Double Dim K1, K2, KW, KB, KF, KS, KP1, KP2, KP3, KSi As Double Dim K(1 To 10) As Double Dim T(1 To 5) As Double Dim TTT, lnK0, K0 As Double Dim lnKBtop, lnKB, lnKW As Double Dim lnKP1, lnKP2, lnKP3, lnKSi, lnK1, lnK2 As Double Dim lnK1fac, lnK2fac, lnKBfac, lnKWfac, lnKSifac As Double Dim lnKFfac, lnKSfac, lnKP1fac, lnKP2fac, lnKP3fac As Double Dim pHGuess, pHTol, pH, deltapH, H, OH As Double Dim ln10, Denom, CAlk, BAlk, PAlk, SiAlk, HSO4, HF As Double Dim FREEtoTOT, HFree, Residual, Slope, PhosTop, PhosBot As Double Dim HCO3, CO3, CO2Ast, fCO2 As Double Dim i As Integer TP = TP * 10 ^ -6 TSi = TSi * 10 ^ -6 TA = TA * 10 ^ -6 TC = TC * 10 ^ -6 TempK = TempC RGasConstant = 83.1451: 'bar-cm3/(mol-K): ' = 8.31451 N-m/(mol-K) TempK = TempC + 273.15 RT = RGasConstant * TempK sqrSal = Sqr(Sal) logTempK = Log(TempK) Pbar = Pdbar / 10 '=========Total Parameters================ 'CalculateTB: ' Uppstrom, L., Deep-Sea Research 21:161-162, 1974: TB = (0.000232 / 10.811) * (Sal / 1.80655): ' in mol/kg-SW ' this is .000416 * Sal / 35. = .0000119 * Sal 'CalculateTF: ' Riley, J. P., Deep-Sea Research 12:219-220, 1965: TF = (0.000067 / 18.998) * (Sal / 1.80655): ' in mol/kg-SW ' this is .000068 * Sal / 35. = .00000195 * Sal 'CalculateTS: ' Morris, A. W., and Riley, J. P., Deep-Sea Research 13:699-705, 1966: TS = (0.14 / 96.062) * (Sal / 1.80655): ' in mol/kg-SW ' this is .02824 * Sal / 35. = .0008067 * Sal '=========Eq Values================ 'CalculateK0: ' Weiss, R. F., Marine Chemistry 2:203-215, 1974. TTT = TempK / 100 lnK0 = -60.2409 + 93.4517 / TTT + 23.3585 * Log(TTT) lnK0 = lnK0 + Sal * (0.023517 - 0.023656 * TTT + 0.0047036 * TTT * TTT) K0 = Exp(lnK0): ' this is in mol/kg-SW/atm 'CalculateIonS: ' This is from the DOE handbook, Chapter 5, p. 13/22, eq. 7.2.4: IonS = 19.924 * Sal / (1000 - 1.005 * Sal) ' ' '*************************************************************************** 'CalculateKS: ' Dickson, A. G., J. Chemical Thermodynamics, 22:113-127, 1990 ' The goodness of fit is .021. ' It was given in mol/kg-H2O. I convert it to mol/kg-SW. ' TYPO!!!!!! on p. 121: the constant e9 should be e8. ' ' This is from eqs 22 and 23 on p. 123, and Table 4 on p 121: lnKS = -4276.1 / TempK + 141.328 - 23.093 * logTempK lnKS = lnKS + (-13856 / TempK + 324.57 - 47.986 * logTempK) * Sqr(IonS) lnKS = lnKS + (35474 / TempK - 771.54 + 114.723 * logTempK) * IonS lnKS = lnKS + (-2698 / TempK) * Sqr(IonS) * IonS lnKS = lnKS + (1776 / TempK) * IonS * IonS KS = Exp(lnKS): ' this is on the free pH scale in mol/kg-H2O KS = KS * (1 - 0.001005 * Sal): ' convert to mol/kg-SW 'CalculateKF: ' Dickson, A. G. and Riley, J. P., Marine Chemistry 7:89-99, 1979: lnKF = 1590.2 / TempK - 12.641 + 1.525 * Sqr(IonS) KF = Exp(lnKF): ' this is on the free pH scale in mol/kg-H2O KF = KF * (1 - 0.001005 * Sal): ' convert to mol/kg-SW SWStoTOT = (1 + TS / KS) / (1 + TS / KS + TF / KF) 'CalculateKB: ' Dickson, A. G., Deep-Sea Research 37:755-766, 1990: lnKBtop = -8966.9 - 2890.53 * sqrSal - 77.942 * Sal lnKBtop = lnKBtop + 1.728 * sqrSal * Sal - 0.0996 * Sal * Sal lnKB = lnKBtop / TempK lnKB = lnKB + 148.0248 + 137.1942 * sqrSal + 1.62142 * Sal lnKB = lnKB + (-24.4344 - 25.085 * sqrSal - 0.2474 * Sal) * logTempK lnKB = lnKB + 0.053105 * sqrSal * TempK KB = Exp(lnKB): ' this is on the total pH scale in mol/kg-SW KB = KB / SWStoTOT: ' convert to SWS pH scale 'CalculateKW: ' Millero, Geochemica et Cosmochemica Acta 59:661-677, 1995. ' his check value of 1.6 umol/kg-SW should be 6.2 lnKW = 148.9802 - 13847.26 / TempK - 23.6521 * logTempK lnKW = lnKW + (-5.977 + 118.67 / TempK + 1.0495 * logTempK) * sqrSal lnKW = lnKW - 0.01615 * Sal KW = Exp(lnKW): ' this is on the SWS pH scale in (mol/kg-SW)^2 'CalculateKP1KP2KP3KSi: ' Yao and Millero, Aquatic Geochemistry 1:53-88, 1995 ' KP1, KP2, KP3 are on the SWS pH scale in mol/kg-SW. ' KSi was given on the SWS pH scale in molal units. ' lnKP1 = -4576.752 / TempK + 115.54 - 18.453 * logTempK lnKP1 = lnKP1 + (-106.736 / TempK + 0.69171) * sqrSal lnKP1 = lnKP1 + (-0.65643 / TempK - 0.01844) * Sal KP1 = Exp(lnKP1) ' ' lnKP2 = -8814.715 / TempK + 172.1033 - 27.927 * logTempK lnKP2 = lnKP2 + (-160.34 / TempK + 1.3566) * sqrSal lnKP2 = lnKP2 + (0.37335 / TempK - 0.05778) * Sal KP2 = Exp(lnKP2) ' ' lnKP3 = -3070.75 / TempK - 18.126 lnKP3 = lnKP3 + (17.27039 / TempK + 2.81197) * sqrSal lnKP3 = lnKP3 + (-44.99486 / TempK - 0.09984) * Sal KP3 = Exp(lnKP3) ' ' lnKSi = -8904.2 / TempK + 117.4 - 19.334 * logTempK lnKSi = lnKSi + (-458.79 / TempK + 3.5913) * Sqr(IonS) lnKSi = lnKSi + (188.74 / TempK - 1.5998) * IonS lnKSi = lnKSi + (-12.1652 / TempK + 0.07871) * IonS * IonS KSi = Exp(lnKSi): ' this is on the SWS pH scale in mol/kg-H2O KSi = KSi * (1 - 0.001005 * Sal): ' convert to mol/kg-SW 'CalculateK1K2: '************************************************ ' MEHRBACH refit BY DICKSON AND MILLERO ' Dickson and Millero, Deep-Sea Research, 34(10):1733-1743, 1987 ' (see also Corrigenda, Deep-Sea Research, 36:983, 1989) ' refit data of Mehrbach et al, Limn Oc, 18(6):897-907, 1973 ' on the SWS pH scale in mol/kg-SW. ' Mehrbach et al gave results on the NBS scale. ' The 2s precision in pK1 is .011, or 2.6% in K1. ' The 2s precision in pK2 is .020, or 4.6% in K2. ' ' This is in Table 4 on p. 1739. pK1 = 3670.7 / TempK - 62.008 + 9.7944 * logTempK pK1 = pK1 - 0.0118 * Sal + 0.000116 * Sal * Sal K1 = 10 ^ (-pK1): ' this is on the SWS pH scale in mol/kg-SW ' ' This is in Table 4 on p. 1739. pK2 = 1394.7 / TempK + 4.777 - 0.0184 * Sal + 0.000118 * Sal * Sal ' If TCO2m >= 2050 Then ' pK2 = pK2 - 0.00015 * (TCO2m - 2050) ' End If K2 = 10 ^ (-pK2): ' this is on the SWS pH scale in mol/kg-SW '==============PCorrections============= 'CorrectK1K2KBForPressure: 'PressureEffectsOnK1: ' These are from Millero, 1995. ' They are the same as Millero, 1979 and Millero, 1992. ' They are from data of Culberson and Pytkowicz, 1968. deltaV = -25.5 + 0.1271 * TempC 'deltaV = deltaV - .151 * (Sal - 34.8): ' Millero, 1979 Kappa = (-3.08 + 0.0877 * TempC) / 1000 'Kappa = Kappa - .578 * (Sal - 34.8)/1000.: ' Millero, 1979 lnK1fac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' The fits given in Millero, 1983 are somewhat different. ' ' PressureEffectsOnK2: ' These are from Millero, 1995. ' They are the same as Millero, 1979 and Millero, 1992. ' They are from data of Culberson and Pytkowicz, 1968. deltaV = -15.82 - 0.0219 * TempC 'deltaV = deltaV + .321 * (Sal - 34.8): ' Millero, 1979 Kappa = (1.13 - 0.1475 * TempC) / 1000 'Kappa = Kappa - .314 * (Sal - 34.8) / 1000: ' Millero, 1979 lnK2fac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' The fit given in Millero, 1983 is different. ' Not by a lot for deltaV, but by much for Kappa!!!. ' ' ' 'PressureEffectsOnKB: ' This is from Millero, 1979. ' It is from data of Culberson and Pytkowicz, 1968. deltaV = -29.48 + 0.1622 * TempC + 0.002608 * TempC * TempC ' Millero, 1983 has: 'deltaV = -28.56 + .1211 * TempC - .000321 * TempC * TempC ' Millero, 1992 has: 'deltaV = -29.48 + .1622 * TempC + .295 * (Sal - 34.8) ' Millero, 1995 has: 'deltaV = -29.48 - .1622 * TempC - .002608 * TempC * TempC 'deltaV = deltaV + .295 * (Sal - 34.8): ' Millero, 1979 Kappa = -2.84 / 1000: ' Millero, 1979 ' Millero, 1992 and Millero, 1995 also have this. 'Kappa = Kappa + .354 * (Sal - 34.8) / 1000: ' Millero,1979 ' Millero, 1983 has: 'Kappa = (-3 + .0427 * TempC) / 1000 lnKBfac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT 'CorrectKWForPressure: ' GEOSECS doesn't include OH term, so this won't matter. ' Peng et al didn't include pressure, but here I assume that the KW correction ' is the same as for the other seawater cases. 'PressureEffectsOnKW: ' This is from Millero, 1983 and his programs CO2ROY(T).BAS. deltaV = -20.02 + 0.1119 * TempC - 0.001409 * TempC * TempC ' Millero, 1992 and Millero, 1995 have: 'deltaV = -25.6 + .2324*TempC - .0036246*TempC*TempC ' This is the freshwater value listed in Millero, 1983. ' The difference is about 4 to 5 over the range 0 < TempC < 20, ' which corresponds to a change in KW(P) of 3% at 200 bar, ' 8% at 500 bar, and 18% at 1000 bar. ' This is probably correct since in Millero, 1983 values of ' -deltaVs are less in seawater than pure water in all cases. Kappa = (-5.13 + 0.0794 * TempC) / 1000: ' Millero, 1983 ' Millero, 1995 has this too, but Millero, 1992 is different. lnKWfac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' Millero, 1979 does not list values for these. 'PressureEffectsOnKF: ' This is from Millero, 1995, which is the same as Millero, 1983. ' It is assumed that KF is on the free pH scale. deltaV = -9.78 - 0.009 * TempC - 0.000942 * TempC * TempC Kappa = (-3.91 + 0.054 * TempC) / 1000 lnKFfac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' ' 'PressureEffectsOnKS: ' This is from Millero, 1995, which is the same as Millero, 1983. ' It is assumed that KS is on the free pH scale. deltaV = -18.03 + 0.0466 * TempC + 0.000316 * TempC * TempC Kappa = (-4.53 + 0.09 * TempC) / 1000 lnKSfac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT CorrectKP1KP2KP3KSiForPressure: ' These corrections don't matter for the GEOSECS choice (WhichKs% = 6) and ' the freshwater choice (WhichKs% = 8). For the Peng choice I assume ' that they are the same as for the other choices (WhichKs% = 1 to 5). ' The corrections for KP1, KP2, and KP3 are from Millero, 1995, which are the ' same as Millero, 1983. ' ' 'PressureEffectsOnKP1: deltaV = -14.51 + 0.1211 * TempC - 0.000321 * TempC * TempC Kappa = (-2.67 + 0.0427 * TempC) / 1000 lnKP1fac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' ' 'PressureEffectsOnKP2: deltaV = -23.12 + 0.1758 * TempC - 0.002647 * TempC * TempC Kappa = (-5.15 + 0.09 * TempC) / 1000 lnKP2fac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' ' 'PressureEffectsOnKP3: deltaV = -26.57 + 0.202 * TempC - 0.003042 * TempC * TempC Kappa = (-4.08 + 0.0714 * TempC) / 1000 lnKP3fac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' ' 'PressureEffectsOnKSi: ' !!!! The only mention of this is Millero, 1995 where it is stated that the ' values have been estimated from the values of boric acid. HOWEVER, ' there is no listing of the values in the table. ' I used the values for boric acid from above. deltaV = -29.48 + 0.1622 * TempC - 0.002608 * TempC * TempC Kappa = -2.84 / 1000 lnKSifac = (-deltaV + 0.5 * Kappa * Pbar) * Pbar / RT ' ' '************************************************ CorrectKsForPressureHere: K1fac = Exp(lnK1fac): K1 = K1 * K1fac K2fac = Exp(lnK2fac): K2 = K2 * K2fac KWfac = Exp(lnKWfac): KW = KW * KWfac KBfac = Exp(lnKBfac): KB = KB * KBfac KFfac = Exp(lnKFfac): KF = KF * KFfac KSfac = Exp(lnKSfac): KS = KS * KSfac KP1fac = Exp(lnKP1fac): KP1 = KP1 * KP1fac KP2fac = Exp(lnKP2fac): KP2 = KP2 * KP2fac KP3fac = Exp(lnKP3fac): KP3 = KP3 * KP3fac KSifac = Exp(lnKSifac): KSi = KSi * KSifac '=============================================== K(1) = K1: K(2) = K2: K(3) = KW: K(4) = KB: K(5) = KF K(6) = KS: K(7) = KP1: K(8) = KP2: K(9) = KP3: K(10) = KSi 'ConvertFromSWSpHScaleToChosenScale: SWStoTOT = (1 + TS / KS) / (1 + TS / KS + TF / KF) pHfactor = SWStoTOT For i = 1 To 4 K(i) = K(i) * pHfactor Next i ' KS and KF remain on the free pH scale For i = 7 To 10 K(i) = K(i) * pHfactor Next i ' SUB CalculatepHfromTATC, version 04.01, 10-13-96, written by Ernie Lewis. ' Inputs: TA, TC, K(), T() ' Output: pH ' This calculates pH from TA and TC using K1 and K2 by Newton's method. ' It tries to solve for the pH at which Residual = 0. ' The starting guess is pH = 8. ' Though it is coded for H on the total pH scale, for the pH values occuring ' in seawater (pH > 6) it will be equally valid on any pH scale (H terms ' negligible) as long as the K constants are on that scale. ' ' K1 = K(1): K2 = K(2): KW = K(3): KB = K(4): KF = K(5) KS = K(6): KP1 = K(7): KP2 = K(8): KP3 = K(9): KSi = K(10) ' TB = T(1): TF = T(2): TS = T(3): TP = T(4): TSi = T(5) ' ' pHGuess = 8!: ' this is the first guess pHTol = 0.0001: ' this is .0001 pH units ln10 = Log(10) pH = pHGuess Do H = 10 ^ (-pH) Denom = (H * H + K1 * H + K1 * K2) CAlk = TC * K1 * (H + 2 * K2) / Denom BAlk = TB * KB / (KB + H) OH = KW / H PhosTop = KP1 * KP2 * H + 2 * KP1 * KP2 * KP3 - H * H * H PhosBot = H * H * H + KP1 * H * H + KP1 * KP2 * H + KP1 * KP2 * KP3 PAlk = TP * PhosTop / PhosBot SiAlk = TSi * KSi / (KSi + H) FREEtoTOT = (1 + TS / KS): ' pH scale conversion factor HFree = H / FREEtoTOT: ' for H on the total scale HSO4 = TS / (1 + KS / HFree): ' since KS is on the free scale HF = TF / (1 + KF / HFree): ' since KF is on the free scale Residual = TA - CAlk - BAlk - OH - PAlk - SiAlk + HFree + HSO4 + HF ' ' find Slope dTA/dpH: ' (this is not exact, but keeps all important terms): Slope = ln10 * (TC * K1 * H * (H * H + K1 * K2 + 4 * H * K2) / Denom / Denom + BAlk * H / (KB + H) + OH + H) deltapH = Residual / Slope: ' this is Newton's method ' to keep the jump from being too big: Do While Abs(deltapH) > 1: deltapH = deltapH / 2!: Loop pH = pH + deltapH Loop While Abs(deltapH) > pHTol HCO3 = TC * K1 * H / (H ^ 2 + K1 * H + K1 * K2) CO3 = TC * K1 * K2 / (H ^ 2 + K1 * H + K1 * K2) CO2Ast = H * HCO3 / K1 fCO2 = CO2Ast / K0 HCO3_T = HCO3 * 1000000 End Function