Function pH_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 pH from TC and TA where ...
' K1 & K2: MEHRBACH refit BY DICKSON AND MILLERO, 1987;
' SO4: Dickson;
' Total pH scale
'Dim Sal, TempC, 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
pH_T = pH
End Function