import numpy as np # Calculate k0ext for each calibration point - from k0ext_from_Vext_pHcal.py # https://github.com/SUPScientist/SCS_Analysis/blob/main/k0ext_from_Vext_pHcal.py # modified to include slainity as an input def k0ext_from_Vext_pHcal(Vext, pHcal, T_C, calsal): ''' Inputs Vext (V): External ref voltage pHcal (total scale): Calibration pH T_C (ºC): Measured T during calibration calsal (PSU): Measured S during calibration Outputs k0ext (V): defined as calibration constant at 0 C k0ext_insitu (V): temperature-corrected calibration value ''' # Constants # Universal gas constant, (R) , http://physics.nist.gov/cgi-bin/cuu/Value?r R = 8.31451 # J/(mol K) F = 96487 # Faraday constant Coulomb / mol k2_int = -0.00125 k2_ext = -0.001048 # Lumped, converted quantities T_K = T_C+273.15 # used often enough to warrant new variable S_Nernst = (R*T_K)/F*np.log(10) # ditto ^^^ # Physical chemistry calcs based on temp and salinity Z = 19.924*calsal/(1000-1.005*calsal) # Ionic strength, Dickson et al. 2007 SO4_tot = (0.14/96.062)*(calsal/1.80655) # Total conservative sulfate cCl = 0.99889/35.453*calsal/1.80655 # Conservative chloride mCl = cCl*1000/(1000-calsal*35.165/35) # mol/kg-H2O K_HSO4 = np.exp(-4276.1/T_K+141.328-23.093*np.log(T_K) +(-13856/T_K+324.57-47.986*np.log(T_K))*Z**0.5 +(35474/T_K-771.54+114.723*np.log(T_K))*Z-2698/T_K*Z**1.5 +1776/T_K*Z**2+np.log(1-0.001005*calsal)) # Bisulfate equilibrium const., Dickson et al. 2007 pHcal_free = pHcal+np.log10(1+SO4_tot/K_HSO4) cHfree = 10**(-pHcal_free) # mol/kg-sw, not used pHcal_free = pHcal_free+np.log10((1000-calsal*35.165/35)/1000) # mol/kg-H2O mHfree = 10**(-pHcal_free) # mol/kg-H2O DHconst = 0.00000343*T_C**2+0.00067524*T_C+0.49172143 # Debye-Huckel, Khoo et al. 1977 log10gamma_HCl = 2*(-DHconst*np.sqrt(Z)/(1+1.394*np.sqrt(Z))+(0.08885-0.000111*T_C)*Z) aHfree_aCl = mHfree*mCl*10**(log10gamma_HCl) # Nernst calibration coefficients, standard-ish potentials k0ext_insitu = Vext+S_Nernst*np.log10(aHfree_aCl) k0ext = k0ext_insitu-k2_ext*T_C # k0 return k0ext_insitu, k0ext #--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- #https://github.com/SUPScientist/SCS_Analysis/blob/main/pHext_from_Vext_k0ext.py def pHext_from_Vext_k0ext(k0ext, Vext, T_C, sal): # See Martz et al. 2010, Bresnahan et al. 2014 for greater detail ''' Inputs External ref voltage (V) Calibration coefficient, k0ext (0 ºC value) Measured T (C) Salinity (PSU) Outputs pHext (total scale) ''' # Constants # Universal gas constant, (R) , http://physics.nist.gov/cgi-bin/cuu/Value?r R = 8.31451 # J/(mol K) F = 96487 # Faraday constant Coulomb / mol k2_int = -0.00125 k2_ext = -0.001048 # Lumped, converted quantities T_K = T_C+273.15 # used often enough to warrant new variable S_Nernst = (R*T_K)/F*np.log(10) # ditto ^^^ # Physical chemistry calcs based on temp and salinity Z = 19.924*sal/(1000-1.005*sal) # Ionic strength, Dickson et al. 2007 SO4_tot = (0.14/96.062)*(sal/1.80655) # Total conservative sulfate cCl = 0.99889/35.453*sal/1.80655 # Conservative chloride mCl = cCl*1000/(1000-sal*35.165/35) # mol/kg-H2O K_HSO4 = np.exp(-4276.1/T_K+141.328-23.093*np.log(T_K) +(-13856/T_K+324.57-47.986*np.log(T_K))*Z**0.5 +(35474/T_K-771.54+114.723*np.log(T_K))*Z-2698/T_K*Z**1.5 +1776/T_K*Z**2+np.log(1-0.001005*sal)) # Bisulfate equilibrium const., Dickson et al. 2007 DHconst = 0.00000343*T_C**2+0.00067524*T_C+0.49172143 # Debye-Huckel, Khoo et al. 1977 log10gamma_HCl = 2*(-DHconst*np.sqrt(Z)/(1+1.394*np.sqrt(Z))+(0.08885-0.000111*T_C)*Z) # Nernstian pH calculations pHext_free = -(((k0ext+k2_ext*(T_C-0))-Vext)-S_Nernst*(np.log10(mCl)+log10gamma_HCl))/S_Nernst # mol/kg-H2O pHext_free = pHext_free-np.log10((1000-sal*35.165/35)/1000) # mol/kg-sw pHext_tot = pHext_free-np.log10(1+SO4_tot/K_HSO4) return pHext_tot #--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Calculate tris pH def pHtris_from_T(T_C, S): T_K = T_C+273.15 # Eq. 18 DelValls and Dickson 1998 pH = (11911.08 - 18.2499*S - 0.039336*S**2)/T_K - 366.27059 + 0.53993607*S + 0.00016329*S**2 + (64.52243 - 0.084041*S)*np.log(T_K) - 0.11149858*T_K return pH #--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- def KT_HCl(Vext, T_C, P, HCl=0.01): '''Calculate KT for K2 calibrations done in 0.01 N HCl and at atmospheric pressure (e.g. not under pressure). Inputs Vext (V): External ref voltage T_C: Temperature of Solution (Celsius) P: Pressure (dBar) HCl concentraion: Set to 0.01 N by default Outputs KT: the temperature dependent calibration coefficient ''' R = 8.3143 # J/(mol K) Gas constant F = 96487 #Faraday constant Coulomb/mol T_K = 273.15 + T_C #Temperature Kelvin ln10 = np.log(10) #natural log of 10 S_Nernst = (R*T_K)/F*np.log(10) #Calculate temperature dependent parameters from matlab script dV0 = 16.38 +0.0896*T_C**2 #partial molar vol HCl DI aHCl = -0.000113*T_C +0.9067 #reference vct = -(R*T_K/F)*np.log((HCl*aHCl)**2) vcp = -dV0 /23061*0.0242*(P-10)/10 *1/1.01 KT = Vext - vcp - vct return KT #--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- def phcalc_jp(Vrs, Press, Temp, Salt, k0, k2, Pcoefs): """ Calculates the pH (free and total) based on the input parameters. Parameters: Vrs : Voltage between reference electrode and ISFET source Press : Pressures in decibars Temp : Temperature in degrees Celsius Salt : Salinity (usually CTD salinity on the PSS) k0 : Sensor reference potential (intercept at Temp = 0C) k2 : Linear temperature coefficient (slope) Pcoefs : Sensor dependent pressure coefficients Returns: phfree : pH on the free scale (mol/kg H2O) phtot : pH on the total scale (mol/kg seawater) Translated from matlab to python """ # Constants R = 8.31446 # J/(mol K), Universal gas constant F = 96485 # C/mol, Faraday constant Tk = 273.15 + Temp # Temperature in Kelvin ln10 = np.log(10) # Natural log of 10 # Ionic strength of seawater (mol/kg H2O) IonS = 19.924 * Salt / (1000 - 1.005 * Salt) # Mean seawater sulfate concentration (mol/kg solution) Stotal = (0.14 / 96.062) * (Salt / 1.80655) # Mean seawater chloride concentration (mol/kg H2O) Cltotal = 0.99889 / 35.453 * Salt / 1.80655 # mol/kg solution Cltotal = Cltotal / (1 - 0.001005 * Salt) # mol/kg H2O # Bisulfate dissociation constant at T, S, and ionic strength (mol/kg solution) Khso4 = np.exp(-4276.1 / Tk + 141.328 - 23.093 * np.log(Tk) + (-13856 / Tk + 324.57 - 47.986 * np.log(Tk)) * np.sqrt(IonS) + (35474 / Tk - 771.54 + 114.723 * np.log(Tk)) * IonS - 2698 / Tk * IonS ** 1.5 + 1776 / Tk * IonS ** 2 + np.log(1 - 0.001005 * Salt)) # Pressure correction factors deltaVHSO4 = -18.03 + 0.0466 * Temp + 0.000316 * Temp ** 2 KappaHSO4 = (-4.53 + 0.09 * Temp) / 1000 lnKhso4fac = (-deltaVHSO4 + 0.5 * KappaHSO4 * (Press / 10)) * \ (Press / 10) / (R * 10 * Tk) Khso4TPS = Khso4 * np.exp(lnKhso4fac) # Gamma +/- HCl at T and S ADH = (3.4286e-6 * Temp ** 2 + 6.7524e-4 * Temp + 0.49172143) log10gammaHCl = -ADH * np.sqrt(IonS) / (1 + 1.394 * np.sqrt(IonS)) + \ (0.08885 - 0.000111 * Temp) * IonS deltaVHcl = 17.85 + 0.1044 * Temp - 0.001316 * Temp ** 2 ThermoPress = -deltaVHcl * 0.0242 / (23061 * 1.01) * Press / 10 log10gammaHCLtP = log10gammaHCl # Sensor constant E0T = k0 + k2 * Temp # Calculate pressure correction (polynomial function of pressure) pc = np.append(Pcoefs, 0) # Add 0 for the n+1 term (descending powers) pcorr = np.polyval(pc, Press) E0TP = E0T + pcorr # Calculate pH on the free scale (mol/kg H2O) phfree = (Vrs - E0TP - ThermoPress) / (R * Tk / F * ln10) + \ np.log(Cltotal) / ln10 + 2 * log10gammaHCLtP # Convert pH free to mol/kg seawater scale phfree = phfree - np.log10(1 - 0.001005 * Salt) # Calculate pH on the total scale (mol/kg seawater) phtot = phfree - np.log10(1 + Stotal / Khso4TPS) return phfree, phtot # Calculate k0ext for each calibration point - from k0ext_from_Vext_pHcal.py # https://github.com/SUPScientist/SCS_Analysis/blob/main/k0ext_from_Vext_pHcal.py # modified to include slainity as an input def k0ext_from_Vext_and_k2_pHcal(Vext, pHcal, k2_ext, T_C, calsal): ''' Inputs Vext (V): External ref voltage pHcal (total scale): Calibration pH T_C (ºC): Measured T during calibration calsal (PSU): Measured S during calibration Outputs k0ext (V): defined as calibration constant at 0 C k0ext_insitu (V): temperature-corrected calibration value ''' # Constants # Universal gas constant, (R) , http://physics.nist.gov/cgi-bin/cuu/Value?r R = 8.31451 # J/(mol K) F = 96487 # Faraday constant Coulomb / mol #k2_int = -0.00125 #k2_ext = -0.001048 # Lumped, converted quantities T_K = T_C+273.15 # used often enough to warrant new variable S_Nernst = (R*T_K)/F*np.log(10) # ditto ^^^ # Physical chemistry calcs based on temp and salinity Z = 19.924*calsal/(1000-1.005*calsal) # Ionic strength, Dickson et al. 2007 SO4_tot = (0.14/96.062)*(calsal/1.80655) # Total conservative sulfate cCl = 0.99889/35.453*calsal/1.80655 # Conservative chloride mCl = cCl*1000/(1000-calsal*35.165/35) # mol/kg-H2O K_HSO4 = np.exp(-4276.1/T_K+141.328-23.093*np.log(T_K) +(-13856/T_K+324.57-47.986*np.log(T_K))*Z**0.5 +(35474/T_K-771.54+114.723*np.log(T_K))*Z-2698/T_K*Z**1.5 +1776/T_K*Z**2+np.log(1-0.001005*calsal)) # Bisulfate equilibrium const., Dickson et al. 2007 pHcal_free = pHcal+np.log10(1+SO4_tot/K_HSO4) cHfree = 10**(-pHcal_free) # mol/kg-sw, not used pHcal_free = pHcal_free+np.log10((1000-calsal*35.165/35)/1000) # mol/kg-H2O mHfree = 10**(-pHcal_free) # mol/kg-H2O DHconst = 0.00000343*T_C**2+0.00067524*T_C+0.49172143 # Debye-Huckel, Khoo et al. 1977 log10gamma_HCl = 2*(-DHconst*np.sqrt(Z)/(1+1.394*np.sqrt(Z))+(0.08885-0.000111*T_C)*Z) aHfree_aCl = mHfree*mCl*10**(log10gamma_HCl) # Nernst calibration coefficients, standard-ish potentials k0ext_insitu = Vext+S_Nernst*np.log10(aHfree_aCl) k0ext = k0ext_insitu-k2_ext*T_C # k0 return k0ext_insitu, k0ext