Nothing
R/chem_const.R
In SolveSAPHE: Solver Suite for Alkalinity-PH Equations
#
# Copyright 2013, 2014, 2020, 2021 Guy Munhoven
#
# This file is part of SolveSAPHE v. 2
# SolveSAPHE is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# SolveSAPHE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with SolveSAPHE. If not, see <http://www.gnu.org/licenses/>.
#
# --------------------------------------------------------
# List of subroutines for the chemical constants (PRIVATE)
# --------------------------------------------------------
# AK_CARB_0_WEIS74
# AK_CARB_1_MILL95, AK_CARB_2_MILL95
# AK_CARB_1_LUEK00, AK_CARB_2_LUEK00
# AK_CARB_1_ROYE93, AK_CARB_2_ROYE93
# AK_BORA_DICK90
# AK_W_MILL95
# AK_PHOS_1_MILL95, AK_PHOS_2_MILL95, AK_PHOS_3_MILL95
# AK_SILI_1_MILL95
# AK_H2S_1_MILL95
# AK_AMMO_1_YAMI95
# ABETA_HF_DIRI79
# AK_HF_PEFR87
# AK_HSO4_DICK90
# A_IONSTRENGTH_SALIN
# ASP_CALC_MUCC83, ASP_ARAG_MUCC83
# A_BTOT_SALIN, A_CATOT_SALIN, A_FTOT_SALIN, A_SO4TOT_SALIN
# ACVT_HSWS_O_HTOT, ACVT_HTOT_O_HFREE, ACVT_HSWS_O_HFREE
# A_RHOSW1_MUNH97, A_RHOSW2_MUNH97
# --------------------------------------
# Parameters for usage within the module
# --------------------------------------
# Gas constant
# ------------
#gasconst_bar_cm3_o_mol_k = 83.14510 # libthdyct
gasconst_bar_cm3_o_mol_k = 83.14472 # Handbook (2007)
# 0 degrees centigrade in Kelvin
# ------------------------------
t_k_zerodegc = 273.15 # Handbook (2007)
#=======================================================================
compute_dissoc_constants <- function (p_temp, p_sal, p_pres, p_pHscale="T")
#=======================================================================
#
# This function computes and returns dissociation constants
# of main chemical elements dissolved in sea water
#
# Input parameters :
# p_temp : temperature in degree Celsius
# p_sal : salinity in psu (Practical Salinity Unit)
# p_pres : pressure in bar
# p_pHscale : pH scale: "T" for Total scale (default), "SWS" for SeaWater, "F" for Free
#
# Output : this function returns a named list of all dissociation constants. Member names are :
#
# K1_DIC : First dissociation constant of carbonic acid (mol/kg) on chosen scale
# K2_DIC : Second dissociation constant of carbonic acid (mol/kg) on chosen scale
# K_BT : Dissociation constant of boric acid (mol/kg) on chosen scale
# K1_PO4 : First dissociation constant of phosphoric acid (mol/kg) on chosen scale
# K2_PO4 : Second dissociation constant of phosphoric acid (mol/kg) on chosen scale
# K3_PO4 : third dissociation constant of phosphoric acid (mol/kg) on chosen scale
# K_Sil : Dissociation constant of sillicic acid (mol/kg) on chosen scale
# K_NH4 : Dissociation constant of ammonium (mol/kg) on chosen scale
# K_H2S : Dissociation constant of hydrogen sulfide (mol/kg) on chosen scale
# K_HSO4 : Dissociation constant of hydrogen sulfate (mol/kg) on free scale
# K_HF : Dissociation constant of hydrogen fluoride (mol/kg) on free scale
# K_H2O : Dissociation constant of water (mol/kg) on chosen scale
#
# Note : K2_Sil, second dissociation constant of sillicic acid is not computed.
{
t_k = p_temp + t_k_zerodegc
dissoc <- list()
if (p_pHscale == "SWS")
{
zcvt_hsws_o_htot = ACVT_HSWS_O_HTOT(t_k, p_sal, p_pres)
dissoc$K1_DIC <- AK_CARB_1_MILL95(t_k, p_sal, p_pres)
dissoc$K2_DIC <- AK_CARB_2_MILL95(t_k, p_sal, p_pres)
dissoc$K_BT <- AK_BORA_DICK90(t_k, p_sal, p_pres) * zcvt_hsws_o_htot
dissoc$K1_PO4 <- AK_PHOS_1_MILL95(t_k, p_sal, p_pres)
dissoc$K2_PO4 <- AK_PHOS_2_MILL95(t_k, p_sal, p_pres)
dissoc$K3_PO4 <- AK_PHOS_3_MILL95(t_k, p_sal, p_pres)
dissoc$K_Sil <- AK_SILI_1_MILL95(t_k, p_sal)
dissoc$K_NH4 <- AK_AMMO_1_YAMI95(t_k, p_sal, p_pres)
dissoc$K_H2S <- AK_H2S_1_MILL95(t_k, p_sal, p_pres)
dissoc$K_HSO4 <- AK_HSO4_DICK90(t_k, p_sal, p_pres)
dissoc$K_HF <- 1. / ABETA_HF_DIRI79(t_k, p_sal, p_pres)
dissoc$K_H2O <- AK_W_MILL95(t_k, p_sal, p_pres)
}
else if (p_pHscale == "F")
{
zcvt_hfree_o_htot = 1. / ACVT_HTOT_O_HFREE(t_k, p_sal, p_pres)
zcvt_hfree_o_hsws = 1. / ACVT_HSWS_O_HFREE(t_k, p_sal, p_pres)
dissoc$K1_DIC <- AK_CARB_1_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K2_DIC <- AK_CARB_2_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_BT <- AK_BORA_DICK90(t_k, p_sal, p_pres) * zcvt_hfree_o_htot
dissoc$K1_PO4 <- AK_PHOS_1_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K2_PO4 <- AK_PHOS_2_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K3_PO4 <- AK_PHOS_3_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_Sil <- AK_SILI_1_MILL95(t_k, p_sal) * zcvt_hfree_o_hsws
dissoc$K_NH4 <- AK_AMMO_1_YAMI95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_H2S <- AK_H2S_1_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_HSO4 <- AK_HSO4_DICK90(t_k, p_sal, p_pres)
dissoc$K_HF <- 1. / ABETA_HF_DIRI79(t_k, p_sal, p_pres)
dissoc$K_H2O <- AK_W_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
}
else
{
zcvt_htot_o_hsws = 1. / ACVT_HSWS_O_HTOT(t_k, p_sal, p_pres)
zcvt_htot_o_hfree = ACVT_HTOT_O_HFREE(t_k, p_sal, p_pres)
dissoc$K1_DIC <- AK_CARB_1_LUEK00(t_k, p_sal, p_pres)
dissoc$K2_DIC <- AK_CARB_2_LUEK00(t_k, p_sal, p_pres)
dissoc$K_BT <- AK_BORA_DICK90(t_k, p_sal, p_pres)
dissoc$K1_PO4 <- AK_PHOS_1_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K2_PO4 <- AK_PHOS_2_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K3_PO4 <- AK_PHOS_3_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K_Sil <- AK_SILI_1_MILL95(t_k, p_sal) * zcvt_htot_o_hsws
dissoc$K_NH4 <- AK_AMMO_1_YAMI95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K_H2S <- AK_H2S_1_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K_HSO4 <- AK_HSO4_DICK90(t_k, p_sal, p_pres)
dissoc$K_HF <- AK_HF_PEFR87(t_k, p_sal, p_pres) / zcvt_htot_o_hfree
dissoc$K_H2O <- AK_W_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
}
return (dissoc)
}
#=======================================================================
AK_CARB_0_WEIS74 <- function (t_k, s)
#=======================================================================
{
# Function calculates K0 in (mol/kg-SW)/atmosphere
# References: Weiss (1979) [(mol/kg-SW)/atm]
# pH scale : N/A
# Note : currently no pressure correction
# ------------------
# Argument variables
# ------------------
# s : salinity
# t_k : temperature in K
zt_k_o_100 = t_k/100.
AK_CARB_0_WEIS74 = exp( -60.2409 + 93.4517/zt_k_o_100
+ 23.3585*log(zt_k_o_100)
+ ( 0.023517 - 0.023656*zt_k_o_100
+ 0.0047036*zt_k_o_100*zt_k_o_100)*s )
return (AK_CARB_0_WEIS74)
}
#=======================================================================
AK_CARB_1_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates first dissociation constant of carbonic acid
# in mol/kg-SW on the SWS pH-scale.
# References: Millero (1995, eq 50 -- ln K1(COM))
# Millero (1982) pressure correction
# pH scale: SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kc1_p0 : ln(K_C1) at p_bar = 0
# zln_kc1_pp : pressure correction for p_bar /= 0
# ln(K_C1) value at p_bar = 0
zsqrts = sqrt(s)
zln_kc1_p0 = ( 2.18867 - 2275.0360/t_k - 1.468591*log(t_k)
+ (-0.138681 - 9.33291/t_k)*zsqrts
+ 0.0726483*s
- 0.00574938*s*zsqrts)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -25.50 - 0.151*zds + 0.1271*zt_degc
zdki = ( -3.08 - 0.578*zds + 0.0877*zt_degc)*1.0E-03
zln_kc1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C1 value
AK_CARB_1_MILL95 = exp( zln_kc1_p0 + zln_kc1_pp )
return (AK_CARB_1_MILL95)
}
#=======================================================================
AK_CARB_2_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates second dissociation constant K1
# in mol/kg-SW on the SWS pH-scale.
# References: Millero (1995, eq 51 -- ln K2(COM))
# Millero (1979) pressure correction
# pH scale: SWS
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kc2_p0 : ln(K_C2) at p_bar = 0
# zln_kc2_pp : pressure correction for p_bar /= 0
# ln(K_C2) value at p_bar = 0
zsqrts = sqrt(s)
zln_kc2_p0 = ( -0.84226 - 3741.1288/t_k - 1.437139*log(t_k)
+ (-0.128417 - 24.41239/t_k)*zsqrts
+ 0.1195308*s
- 0.00912840*s*zsqrts)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -15.82 + 0.321*zds - 0.0219*zt_degc
zdki = ( 1.13 - 0.314*zds - 0.1475*zt_degc)*1.0E-03
zln_kc2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C2 value
AK_CARB_2_MILL95 = exp( zln_kc2_p0 + zln_kc2_pp )
return (AK_CARB_2_MILL95)
}
#=======================================================================
AK_CARB_1_LUEK00 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates first dissociation constant of carbonic acid
# in mol/kg-SW on the Total pH-scale.
# References: Luecker et al. (2000) -- also Handbook (2007)
# Millero (1979) pressure correction
# pH scale: Total
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zds : salinity-34.8
# zlog10_kc1_p0 : log_10(k_C1) at p_bar = 0
# zln_kc1_pp : pressure correction for p_bar /= 0
# log_10(K_C1) value at p_bar = 0
zlog10_kc1_p0 = ( 61.2172 - 3633.86/t_k - 9.67770*log(t_k)
+ s*(0.011555 - s*0.0001152))
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -25.50 - 0.151*zds + 0.1271*zt_degc
zdki = ( -3.08 - 0.578*zds + 0.0877*zt_degc)*1.0E-03
zln_kc1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C1 value
AK_CARB_1_LUEK00 = 10.**zlog10_kc1_p0 * exp(zln_kc1_pp)
return (AK_CARB_1_LUEK00)
}
#=======================================================================
AK_CARB_2_LUEK00 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates second dissociation constant K1
# in mol/kg-SW on the Total pH-scale.
# References: Luecker et al. (2000) -- also Handbook (2007)
# Millero (1979) pressure correction
# pH scale: Total
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zlog10_kc2_p0 : log_10(K_C2) at p_bar = 0
# zln_kc2_pp : pressure correction for p_bar /= 0
# log_10(K_C2) value at p_bar = 0
zlog10_kc2_p0 = (-25.9290 - 471.78/t_k + 3.16967*log(t_k)
+ s*(0.01781 - s*0.0001122))
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -15.82 + 0.321*zds - 0.0219*zt_degc
zdki = ( 1.13 - 0.314*zds - 0.1475*zt_degc)*1.0E-03
zln_kc2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C2 value
AK_CARB_2_LUEK00 = 10.**zlog10_kc2_p0 *exp(zln_kc2_pp)
return (AK_CARB_2_LUEK00)
}
#=======================================================================
AK_CARB_1_ROYE93 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates first dissociation constant of carbonic acid
# in mol/kg-SW on the Total pH-scale.
# References: Roy et al. (1993) -- also Handbook (1994)
# Millero (1979) pressure correction
# pH scale : Total
# Note : converted here from mol/kg-H2O to mol/kg-SW
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zds : salinity-34.8
# zln_kc1_p0 : ln(k_C1) at p_bar = 0
# zln_kc1_pp : pressure correction for p_bar /= 0
# ln(K_C1) value at p_bar = 0
zsqrts = sqrt(s)
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zln_kc1_p0 = (-2307.1266/t_k + 2.83655 - 1.5529413*log(t_k)
+ (-4.0484/t_k - 0.20760841)*zsqrts
+ 0.08468345*s
- 0.00654208*zsqrts*s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -25.50 - 0.151*zds + 0.1271*zt_degc
zdki = ( -3.08 - 0.578*zds + 0.0877*zt_degc)*1.0E-03
zln_kc1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C1 value
AK_CARB_1_ROYE93 = exp(zln_kc1_p0 + zln_kc1_pp) * zcvt_to_kgsw
return (AK_CARB_1_ROYE93)
}
#=======================================================================
AK_CARB_2_ROYE93 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates second dissociation constant K1
# in mol/kg-SW on the Total pH-scale.
# References: Roy et al. (1993) -- also Handbook (1994)
# Millero (1979) pressure correction
# pH scale : Total
# Note : converted here from mol/kg-H2O to mol/kg-SW
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kc2_p0 : ln(K_C2) at p_bar = 0
# zln_kc2_pp : pressure correction for p_bar /= 0
# ln(K_C2) value at p_bar = 0
zsqrts = sqrt(s)
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zln_kc2_p0 = (-3351.6106/t_k - 9.226508 - 0.2005743*log(t_k)
+ ( -23.9722/t_k - 0.106901773)*zsqrts
+ 0.1130822*s
- 0.00846934*zsqrts*s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -15.82 + 0.321*zds - 0.0219*zt_degc
zdki = ( 1.13 - 0.314*zds - 0.1475*zt_degc)*1.0E-03
zln_kc2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C2 value
AK_CARB_2_ROYE93 = exp(zln_kc2_p0 + zln_kc2_pp) * zcvt_to_kgsw
return (AK_CARB_2_ROYE93)
}
#=======================================================================
AK_BORA_DICK90 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates boric acid dissociation constant KB
# in mol/kg-SW on the total pH-scale.
# References: Dickson (1990, eq. 23) -- also Handbook (2007, eq. 37)
# Millero (1979) pressure correction
# pH scale : total
# ------------------
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kb_p0 : K_b at p_bar = 0
# zln_kb_pp : pressure correction for p_bar /= 0
# ln(K_B) value at p_bar = 0
zsqrts = sqrt(s)
zln_kb_p0 = (( -8966.90
+ zsqrts*( -2890.53
+ zsqrts*( -77.942
+ zsqrts*( 1.728 - 0.0996*zsqrts)))) / t_k
+ 148.0248 + zsqrts*(137.1942 + zsqrts*1.62142)
+ (-24.4344 + zsqrts*(-25.085 - zsqrts*0.2474)) * log(t_k)
+ 0.053105*zsqrts*t_k)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -29.48 + 0.295*zds + 0.1622*zt_degc - 0.002608*zt_degc*zt_degc
zdki = (-2.84 + 0.354*zds)*1.0E-03
zln_kb_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_B value
AK_BORA_DICK90 = exp( zln_kb_p0 + zln_kb_pp )
return (AK_BORA_DICK90)
}
#=======================================================================
AK_W_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates water dissociation constant Kw in (mol/kg-SW)^2
# References: Millero (1995) for value at p_bar = 0
# Millero (pers. comm. 1996) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kw_p0 : ln(K_w) at p_bar = 0
# zln_kw_pp : pressure correction for p_bar /= 0
# ln(K_w) value at p_bar = 0
zln_kw_p0 = ( 148.9802
- 13847.26/t_k
- 23.6521*log(t_k)
+ ( -5.977 + 118.67/t_k + 1.0495*log(t_k))*sqrt(s)
- 0.01615*s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -20.02 + 0.1119*zt_degc - 0.1409E-02*zt_degc*zt_degc
zdki = ( -5.13 + 0.0794*zt_degc)*1.0E-03
zln_kw_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_w value
AK_W_MILL95 = exp( zln_kw_p0 + zln_kw_pp )
return (AK_W_MILL95)
}
#=======================================================================
AK_PHOS_1_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the first dissociation constant
# of phosphoric acid (H3PO4) in seawater
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp1_p0 : ln(K_p1) at p_bar = 0
# zln_kp1_pp : pressure correction for p_bar /= 0
# ln(K_P1) for p_bar = 0
zln_kp1_p0 = ( 115.54 - 4576.752/t_k - 18.453*log(t_k)
+ ( 0.69171 - 106.736/t_k)* sqrt(s)
+ (-0.01844 - 0.65643/t_k)*s )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -14.51 + 0.1211*zt_degc - 0.321E-03*zt_degc*zt_degc
zdki = ( -2.67 + 0.0427*zt_degc)*1.0E-03
zln_kp1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of K_P1
AK_PHOS_1_MILL95 = exp(zln_kp1_p0 + zln_kp1_pp)
return (AK_PHOS_1_MILL95)
}
#=======================================================================
AK_PHOS_2_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the second dissociation constant
# of phosphoric acid (H3PO4) in seawater
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp2_p0 : ln(K_P2) at p_bar = 0
# zln_kp2_pp : pressure correction for p_bar /= 0
# ln(K_P2) for p_bar = 0
zln_kp2_p0 = ( 172.1033
- 8814.715/t_k
- 27.927*log(t_k)
+ ( 1.3566 - 160.340/t_k)*sqrt(s)
+ (-0.05778 + 0.37335/t_k)*s )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -23.12 + 0.1758*zt_degc -2.647E-03*zt_degc*zt_degc
zdki = ( -5.15 + 0.09*zt_degc)*1.0E-03
zln_kp2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_P2 value
AK_PHOS_2_MILL95 = exp( zln_kp2_p0 + zln_kp2_pp )
return (AK_PHOS_2_MILL95)
}
#=======================================================================
AK_PHOS_3_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the third dissociation constant
# of phosphoric acid (H3PO4) in seawater
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp3_p0 : ln(K_P3) at p_bar = 0
# zln_kp3_pp : pressure correction for p_bar /= 0
# ln(K_P3) for p_bar = 0
zln_kp3_p0 = ( -18.126 - 3070.75/t_k
+ ( 2.81197 + 17.27039/t_k)*sqrt(s)
+ (-0.09984 - 44.99486/t_k)*s )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -26.57 + 0.2020*zt_degc -3.042E-03*zt_degc*zt_degc
zdki = ( -4.08 + 0.0714*zt_degc)*1.0E-03
zln_kp3_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_P3 value
AK_PHOS_3_MILL95 = exp( zln_kp3_p0 + zln_kp3_pp )
return (AK_PHOS_3_MILL95)
}
#=======================================================================
AK_SILI_1_MILL95 <- function (t_k, s)
#=======================================================================
{
# Function returns the first dissociation constant
# of silicic acid (H4SiO4) in seawater
# References: Yao and Millero (1995) cited by Millero (1995)
# pH scale : SWS (according to Dickson et al, 2007)
# Note : No pressure correction available
# Note : converted here from mol/kg-H2O to mol/kg-sw
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# ---------------
# Local variables
# ---------------
# zcvt_to_kgsw: fraction of pure water in 1 kg seawater at salinity s
# zionst : ionic strength [mol/kg-H2O]
# zln_ksi1_p0 : ln(K_Si1) at p_bar = 0
# zln_ksi1_pp : pressure correciotn for p_bar /= 0
# K_Si1 value at p_bar = 0
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zionst = A_IONSTRENGTH_SALIN(s)/zcvt_to_kgsw # mol/kg-H2O ##
zln_ksi1_p0 = ( 117.40 - 8904.2/t_k - 19.334 * log(t_k)
+ ( 3.5913 - 458.79/t_k) * sqrt(zionst)
+ (-1.5998 + 188.74/t_k) * zionst
+ (0.07871 - 12.1652/t_k) * zionst*zionst )
# Pressure correction : currently none
zln_ksi1_pp = 0.
# Final value
AK_SILI_1_MILL95 = exp( zln_ksi1_p0 + zln_ksi1_pp ) * zcvt_to_kgsw
return (AK_SILI_1_MILL95)
}
#=======================================================================
AK_H2S_1_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the dissociation constant of hydrogen sulfide in sea-water
# References: Millero et al. (1988) (cited by Millero (1995)
# Millero (1995) for pressure correction
# pH scale : - SWS (according to Yao and Millero, 1995, p. 82: "refitted if necessary")
# - Total (according to Lewis and Wallace, 1998)
# Note : we stick to SWS here for the time being
# Note : the fits from Millero (1995) and Yao and Millero (1995)
# derive from Millero et al. (1988), with all the coefficients
# multiplied by -ln(10)
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zt_degc : temperature in degrees Celsius
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kh2s_p0 : ln(K_H2S) at p_bar = 0
# zln_kh2s_pp : pressure correction for p_bar /= 0
# K_H2S value at p_bar = 0
# ------------------------
zln_kh2s_p0 = ( 225.838
- 13275.3/t_k
- 34.6435 * log(t_k)
+ 0.3449*sqrt(s)
- 0.0274*s )
# Pressure correction
# -------------------
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -14.80 + zt_degc*(0.0020 - zt_degc*0.400E-03)
zdki = ( 2.89 + zt_degc*0.054)*1.0E-03
zln_kh2s_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_H2S value
# -----------------
AK_H2S_1_MILL95 = exp( zln_kh2s_p0 + zln_kh2s_pp )
return (AK_H2S_1_MILL95)
}
#=======================================================================
AK_AMMO_1_YAMI95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the dissociation constant
# of ammonium in sea-water [mol/kg-SW]
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zt_degc : temperature in degrees Celsius
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_knh4_p0 : ln(K_NH4) at p_bar = 0
# zln_knh4_pp : pressure correction for p_bar /= 0
# K_NH4 value at p_bar = 0
# ------------------------
zln_knh4_p0 = ( -0.25444 - 6285.33/t_k + 0.0001635*t_k
+ ( 0.46532 - 123.7184/t_k) * sqrt(s)
+ (-0.01992 + 3.17556/t_k) * s )
# Pressure correction
# -------------------
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -26.43 + zt_degc*(0.0889 - zt_degc*0.905E-03)
zdki = ( -5.03 + zt_degc*0.0814)*1.0E-03
zln_knh4_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_NH4 value
# -----------------
AK_AMMO_1_YAMI95 = exp( zln_knh4_p0 + zln_knh4_pp )
return (AK_AMMO_1_YAMI95)
}
#=======================================================================
ACVT_KGH2O_O_KGSW <- function (s)
#=======================================================================
{
# Function returns the mass of pure water in one kg of seawater
# of salinity s
# References: "libthdyct" -- derived by Munhoven (1997) from data by Millero (1982)
# "Handbook (2007)" -- Handbook (2007)
# pH scale: N/A
#ACVT_KGH2O_O_KGSW = 1. - 0.0010049*s # libthdyct
ACVT_KGH2O_O_KGSW = 1. - 0.001005*s # Handbook (2007)
return (ACVT_KGH2O_O_KGSW)
}
#=======================================================================
A_IONSTRENGTH_SALIN <- function (s)
#=======================================================================
{
# Function calculates ionic strength in mol/kg-SW, for given salinity.
# References: "libthdyct" -- derived by Munhoven (1997) from data by Millero (1982)
# "Handbook (2007)" -- Handbook (2007)
# pH scale: N/A
#A_IONSTRENGTH_SALIN = (0.019920*s) # libthdyct
A_IONSTRENGTH_SALIN = (0.019924*s) # Handbook (2007)
return (A_IONSTRENGTH_SALIN)
}
#=======================================================================
ABETA_HF_DIRI79 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates association constant \beta_{HF} [(mol/kg-SW)^{-1}]
# in (mol/kg-SW)^{-1}, where
# \beta_{HF} = \frac{ [HF] }{ [H^{+}] [F^{-}] }
# References: Dickson and Riley (1979)
# Millero (1995) for pressure correction
# pH scale : free
# Note : converted here from mol/kg-H2O to mol/kg-SW
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zionst : ionic strength [mol/kg-H2O]
# zcvt_to_kgsw : mass of pure water in 1kg of seawater as a fct. of salinity
# zln_bhf_p0 : \beta_HF at p_bar = 0
# zln_khf_pp : pressure correction for k_HF = 1/\beta_HF at p_bar /= 0
# \beta_HF at p_bar = 0
# ---------------------
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zionst = A_IONSTRENGTH_SALIN(s)/zcvt_to_kgsw
zln_bhf_p0 = -1590.2/t_k + 12.641 - 1.525*sqrt(zionst)
# Pressure correction
# -------------------
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -9.78 + zt_degc*(-0.0090 - zt_degc*0.942E-03)
zdki = ( -3.91 + zt_degc*0.054)*1.0E-03
zln_khf_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final \beta_HF value
# --------------------
# notice that ln(k_HF(P)) = ln(k_HF(0)) + zln_khf_pp
# <=> -ln(\beta_HF(P)) = -ln(\beta_HF(0)) + zln_khf_pp
# <=> ln(\beta_HF(P)) = ln(\beta_HF(0)) - zln_khf_pp
ABETA_HF_DIRI79 = exp(zln_bhf_p0 - zln_khf_pp ) / zcvt_to_kgsw
return (ABETA_HF_DIRI79)
}
#=======================================================================
AK_HF_PEFR87 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates dissociation constant for hydrogen fluoride
# in mol/kg-SW
# References: Perez and Fraga (1987)
# Millero (1995) for pressure correction
# pH scale : Total (according to Handbook, 2007)
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_khf_p0 : ln(K_HF) at p_bar = 0
# zln_khf_pp : pressure correction for p_bar /= 0
# ln(K_HF) at p_bar = 0
zln_khf_p0 = 874./t_k - 9.68 + 0.111*sqrt(s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -9.78 + zt_degc*(-0.0090 - zt_degc*0.942E-03)
zdki = ( -3.91 + zt_degc*0.054)*1.0E-03
zln_khf_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of K_HF
AK_HF_PEFR87 = exp( zln_khf_p0 + zln_khf_pp )
return (AK_HF_PEFR87)
}
#=======================================================================
AK_HSO4_DICK90 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the dissociation constant of hydrogen sulfate (bisulfate)
# References: Dickson (1990) -- also Handbook (2007)
# Millero (1995) for pressure correction
# pH scale : free
# Note : converted here from mol/kg-H2O to mol/kg-SW
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zionst : ionic strength in mol/-kg-H2O
# zsqrti : square root og ion strength
# zcvt_to_kgsw : mass of pure water in 1kg of seawater as a fct. of salinity
# zln_khso4_p0 : K_HSO4 at p_bar = 0
# zln_khso4_pp : pressure correction for p_bar /= 0
# ln(K_HSO4) at p_bar = 0
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zionst = A_IONSTRENGTH_SALIN(s)/zcvt_to_kgsw
zsqrti = sqrt(zionst)
zln_khso4_p0 = ( -4276.1/t_k + 141.328 - 23.093*log(t_k)
+ (-13856./t_k + 324.57 - 47.986*log(t_k)) * zsqrti
+ ( 35474./t_k - 771.54 + 114.723*log(t_k)) * zionst
- ( 2698./t_k)*zsqrti * zionst
+ ( 1776./t_k)*zionst*zionst )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -18.03 + zt_degc*(0.0466 + zt_degc*0.316E-03)
zdki = ( -4.53 + zt_degc*0.0900)*1.0E-03
zln_khso4_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# ln(K_HSO4) at p_bar = 0
AK_HSO4_DICK90 = zcvt_to_kgsw * exp( zln_khso4_p0 + zln_khso4_pp )
return (AK_HSO4_DICK90)
}
#=======================================================================
ASP_CALC_MUCC83 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns stoechiometric solubility product
# of calcite in seawater
# References: Mucci (1983)
# Millero (1995) for pressure correction
# pH scale : N/A
# Units : (mol/kg-SW)^2
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zsqrts : square root of salinity
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp1_p0 : ln(K_p1) at p_bar = 0
# zln_kp1_pp : pressure correction for p_bar /= 0
zsqrts = sqrt(s)
# log10(Ksp_Calc) for p_bar = 0
zlog10_kspcalc_p0 = ( -171.9065 - 0.077993*t_k
+ 2839.319/t_k + 71.595*log10(t_k)
+ (-0.77712 + 0.0028426*t_k + 178.34/t_k)*zsqrts
- 0.07711*s
+ 0.0041249*s*zsqrts )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -48.76 + 0.5304*zt_degc
zdki = (-11.76 + 0.3692*zt_degc)*1.0E-03
zln_kspcalc_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of Ksp_Calc
ASP_CALC_MUCC83 = 10.**(zlog10_kspcalc_p0) * exp(zln_kspcalc_pp)
return (ASP_CALC_MUCC83)
}
#=======================================================================
ASP_ARAG_MUCC83 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns stoechiometric solubility product
# of aragonite in seawater
# References: Mucci (1983)
# Millero (1979) for pressure correction
# pH scale : N/A
# Units : (mol/kg-SW)^2
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zsqrts : square root of salinity
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp1_p0 : ln(K_p1) at p_bar = 0
# zln_kp1_pp : pressure correction for p_bar /= 0
zsqrts = sqrt(s)
# log10(Ksp_Arag) for p_bar = 0
zlog10_ksparag_p0 = ( -171.945 - 0.077993*t_k
+ 2903.293/t_k + 71.595*log10(t_k)
+ (-0.068393 + 0.0017276*t_k + 88.135/t_k)*zsqrts
- 0.10018*s
+ 0.0059415*s*zsqrts )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -48.76 + 0.5304*zt_degc + 2.8
zdki = (-11.76 + 0.3692*zt_degc)*1.0E-03
zln_ksparag_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of Ksp_Arag
ASP_ARAG_MUCC83 = 10.**(zlog10_ksparag_p0) * exp(zln_ksparag_pp)
return (ASP_ARAG_MUCC83)
}
#=======================================================================
A_BTOT_SALIN <- function (s)
#=======================================================================
{
# Function returns total borate concentration in mol/kg-SW
# given the salinity of a sample
# References: Uppström (1974), cited by Dickson et al. (2007, chapter 5, p 10)
# Millero (1982) cited in Millero (1995)
# pH scale : N/A
A_BTOT_SALIN = 0.000416*(s/35.)
return (A_BTOT_SALIN)
}
#=======================================================================
A_CATOT_SALIN <- function (s)
#=======================================================================
{
# Function returns total calcium concentration in mol/kg-SW
# given the salinity of a sample
# References: Culkin and Cox (1966)
# pH scale : N/A
# (g Ca/kg)/Cl_permil values
# Culkin (1967): 0.0213
# Culkin and Cox (DSR 13 1966): 0.02126 +/- 0.00004 (stdev)
# Riley and Tongudai (Chem Geol 2 1967): 0.02128 +/- 0.00006 (stdev)
# Handbook (2007): 0.02127
# with reference to Riley and Tongudai (1967) (???)
# Here:
# (g Ca/kg)/Cl_permil = 0.02127
# (g Ca)/(mol Ca) = 40.078
# Cl_permil = S/1.80655
# mol Ca/kg = (0.02127/40.078) * (35/1.80655)
#A_CATOT_SALIN = 0.010282*(s/35.)
A_CATOT_SALIN = (0.02127/40.078) * (s/1.80655)
return (A_CATOT_SALIN)
}
#=======================================================================
A_FTOT_SALIN <- function (s)
#=======================================================================
{
# Function returns total calcium concentration in mol/kg-SW
# given the salinity of a sample
# References: Culkin (1965) (???)
# pH scale : N/A
A_FTOT_SALIN = 0.000068*(s/35.)
return (A_FTOT_SALIN)
}
#=======================================================================
A_SO4TOT_SALIN <- function(s)
#=======================================================================
{
# Function returns total sulfate concentration in mol/kg-SW
# given the salinity of a sample
# References: Morris, A.W. and Riley, J.P. (1966) quoted in Handbook (2007)
# pH scale : N/A
#A_SO4TOT_SALIN = 0.028234*(s/35.) # in libthdyct and Thesis
#A_SO4TOT_SALIN = 0.02824*(s/35.) # Handbook (2007, chap 6, p 10, tab 2, col 3)
A_SO4TOT_SALIN = (0.1400/96.062)*(s/1.80655) # Handbook (2007, chap 6, p 10)
return (A_SO4TOT_SALIN)
}
#=======================================================================
ACVT_HSWS_O_HTOT <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the ratio H_SWS/H_Tot as a function of salinity s
# Reference: Munhoven
# pH scale: all
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zso4_tot: total sulfate concentration in mol/kg-SW
# zf_tot : total fluoride concentration in mol/kg-SW
#-----------------------------------------------------------------------
zso4_tot = A_SO4TOT_SALIN(s)
zf_tot = A_FTOT_SALIN(s)
ACVT_HSWS_O_HTOT = (1. + (zf_tot*ABETA_HF_DIRI79(t_k, s, p_bar))
/(1. + zso4_tot/AK_HSO4_DICK90(t_k,s, p_bar)))
return (ACVT_HSWS_O_HTOT)
}
#=======================================================================
ACVT_HTOT_O_HFREE <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the ratio H_Tot/H_free as a function of salinity s
# Reference: Munhoven
# pH scale: N/A
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zso4_tot: total sulfate concentration in mol/kg-SW
#-----------------------------------------------------------------------
zso4_tot = A_SO4TOT_SALIN(s)
ACVT_HTOT_O_HFREE = 1. + zso4_tot/AK_HSO4_DICK90(t_k,s, p_bar)
return (ACVT_HTOT_O_HFREE)
}
#=======================================================================
ACVT_HSWS_O_HFREE <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the ratio H_SWS/H_free as a function
# of salinity s
# Reference: Munhoven
# pH scale: N/A
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zso4_tot: total sulfate concentration in mol/kg-SW
# zf_tot : total fluoride concentration in mol/kg-SW
#-----------------------------------------------------------------------
zso4_tot = A_SO4TOT_SALIN(s)
zf_tot = A_FTOT_SALIN(s)
ACVT_HSWS_O_HFREE = ( 1. + zf_tot*ABETA_HF_DIRI79(t_k, s, p_bar)
+ zso4_tot/AK_HSO4_DICK90(t_k,s, p_bar) )
return (ACVT_HSWS_O_HFREE)
}
#=======================================================================
A_RHOSW1_MUNH97 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns first order approximation of \rho in (kg-SW)/(m^3-SW)
# References: Munhoven (1997)
# after EOS80 (UNESCO, 1981, 1983)
# ------------------
# Argument variables
# ------------------
# s : salinity
# tk : temperature in K
# p_bar : depth in m
# ---------------
# Local variables
# ---------------
s0 = 35.5
t_k0 = 285.16
p_bar0 = 300.0
A_RHOSW1_MUNH97 = ( 1039.9044 + 0.77629393*(s-s0)
- 0.19692738*(t_k-t_k0) )
return (A_RHOSW1_MUNH97)
}
#=======================================================================
A_RHOSW2_MUNH97 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns second order approximation of \rho in (kg-SW)/(m^3-SW)
# References: Munhoven (1997)
# after EOS80 (UNESCO, 1981, 1983)
# ------------------
# Argument variables
# ------------------
# s : salinity
# tk : temperature in K
# p_bar : depth in m
# ---------------
# Local variables
# ---------------
s0 = 35.5
t_k0 = 285.16
p_bar0 = 300.0
A_RHOSW2_MUNH97 =( 1040.0145
+ 0.77629393*(s-s0)
- 0.25013591*(t_k-t_k0)
+ 4.2026266E-02*(p_bar-p_bar0)
- 4.7473116E-03*(t_k-t_k0)*(t_k-t_k0)
- 4.7974224E-06*(p_bar-p_bar0)*(p_bar-p_bar0)
- 2.1404592E-04*(t_k-t_k0)*(p_bar-p_bar0) )
return (A_RHOSW2_MUNH97)
}
#=======================================================================
CHECKCONSTANTS <- function()
#=======================================================================
{
# ---------------
# Local variables
# ---------------
# s : salinity
# tk : temperature in K
# p_bar : applied pressure in bar
print ('Checking constant values generated from MOD_CHEMCONST')
print ('')
print (' % indicates checking against the Handbook (1994);')
print (' $ indicates checking against the Lewis and Wallace (1998);')
print (' * indicates checking against the Handbook (2007);')
print (' target values are quoted in brackets')
print ('')
print ('')
print (' For S = 35, P = 0 and T/K = 298.15:')
print ('')
s = 35.
p_bar = 0.
t_k = 298.15
zkc0 = AK_CARB_0_WEIS74(t_k, s)
print ('')
print ('K_0 -- Weiss (1974)')
print ('===================')
print ('')
print (c(' K_0 :', zkc0))
print (c(' ln(K_0) :', log(zkc0)))
print (c(' pK_0 :', -log10(zkc0)))
print ('')
zkhso4 = AK_HSO4_DICK90(t_k, s, p_bar)
print ('')
print ('K_HSO4 -- Dickson (1990) -- pH_free')
print ('===================================')
print ('')
print (c(' K_HSO4 :', zkhso4))
print (c(' ln(K_HSO4) :', log(zkhso4)))
print (c(' pK_HSO4 :', -log10(zkhso4)))
print ('')
zkb = AK_BORA_DICK90(t_k, s, p_bar)
print ('')
print ('K_b -- Dickson (1990) -- pH_tot')
print ('===============================')
print ('')
print (c(' K_b :', zkb))
print (c(' ln(K_b) :', log(zkb)))
print (c(' pK_b :', -log10(zkb)))
print ('')
zkc1 = AK_CARB_1_LUEK00(t_k, s, p_bar)
print ('')
print ('K_1 -- Luecker et al (2000) -- pH_tot')
print ('=====================================')
print ('')
print (c(' K_1 :', zkc1))
print (c(' ln(K_1) :', log(zkc1)))
print (c(' pK_1 :', -log10(zkc1)))
print ('')
zkc2 = AK_CARB_2_LUEK00(t_k, s, p_bar)
print ('')
print ('K_2 -- Luecker et al (2000) -- pH_tot')
print ('=====================================')
print ('')
print (c(' K_2 :', zkc2))
print (c(' ln(K_2) :', log(zkc2)))
print (c(' pK_2 :', -log10(zkc2)))
print ('')
zkc1 = AK_CARB_1_ROYE93(t_k, s, p_bar)
print ('')
print ('K_1 -- Roy et al (1993) -- pH_tot')
print ('=================================')
print ('')
print (c(' K_1 :', zkc1))
print (c(' ln(K_1) :', log(zkc1)))
print (c(' pK_1 :', -log10(zkc1)))
print ('')
zkc2 = AK_CARB_2_ROYE93(t_k, s, p_bar)
print ('')
print ('K_2 -- Roy et al (1993) -- pH_tot')
print ('=================================')
print ('')
print (c(' K_2 :', zkc2))
print (c(' ln(K_2) :', log(zkc2)))
print (c(' pK_2 :', -log10(zkc2)))
print ('')
zkhf = AK_HF_PEFR87(t_k, s, p_bar)
print ('')
print ('K_HF -- Perez and Fraga (1987) -- pH_tot')
print ('========================================')
print ('')
print (c(' K_HF :', zkhf))
print (c(' ln(K_HF) :', log(zkhf)))
print (c(' pK_HF :', -log10(zkhf)))
print ('')
zkp1 = AK_PHOS_1_MILL95(t_k, s, p_bar)
print ('')
print ('K_P1 -- Millero (1995) -- pH_SWS')
print ('================================')
print ('')
print (c(' K_P1 :', zkp1))
print (c(' ln(K_P1) :', log(zkp1)))
print (c(' pK_1 :', -log10(zkp1)))
print ('')
zkp2 = AK_PHOS_2_MILL95(t_k, s, p_bar)
print ('')
print ('K_P2 -- Millero (1995) -- pH_SWS')
print ('================================')
print ('')
print (c(' K_2 :', zkp2))
print (c(' ln(K_P2) :', log(zkp2)))
print (c(' pK_2 :', -log10(zkp2)))
print ('')
zkp3 = AK_PHOS_3_MILL95(t_k, s, p_bar)
print ('')
print ('K_P3 -- Millero (1995) -- pH_SWS')
print ('================================')
print ('')
print (c(' K_P3 :', zkp3))
print (c(' ln(K_P3) :', log(zkp3)))
print (c(' pK_P3 :', -log10(zkp3)))
print ('')
zksi1 = AK_SILI_1_MILL95(t_k, s)
print ('')
print ('K_Si1 -- Millero (1995) -- pH_SWS')
print ('=================================')
print ('')
print (c(' K_Si1 :', zksi1))
print (c(' ln(K_Si1) :', log(zksi1)))
print (c(' pK_Si1 :', -log10(zksi1)))
print ('')
zkw = AK_W_MILL95(t_k, s, p_bar)
print ('')
print ('K_w -- Millero (1995) -- pH_SWS')
print ('===============================')
print ('')
print (c(' K_w :', zkw))
print (c(' ln(K_w) :', log(zkw)))
print (c(' pK_w :', -log10(zkw)))
print ('')
zkh2s = AK_H2S_1_MILL95(t_k, s, p_bar)
print ('')
print ('K_H2S -- Millero (1995) -- pH_SWS')
print ('=================================')
print ('')
print (c(' K_H2S :', zkh2s))
print (c(' ln(K_H2S) :', log(zkh2s)))
print (c(' pK_H2S :', -log10(zkh2s)))
print ('')
zknh4 = AK_AMMO_1_YAMI95(t_k, s, p_bar)
print ('')
print ('K_NH4 -- Yao and Millero (1995) -- pH_SWS')
print ('=========================================')
print ('')
print (c(' K_NH4 :', zknh4))
print (c(' ln(K_NH4) :', log(zknh4)))
print (c(' pK_NH4 :', -log10(zknh4)))
print ('')
return()
}
Try the SolveSAPHE package in your browser
Any scripts or data that you put into this service are public.
SolveSAPHE documentation built on May 2, 2021, 1:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#
# Copyright 2013, 2014, 2020, 2021 Guy Munhoven
#
# This file is part of SolveSAPHE v. 2
# SolveSAPHE is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# SolveSAPHE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with SolveSAPHE. If not, see <http://www.gnu.org/licenses/>.
#
# --------------------------------------------------------
# List of subroutines for the chemical constants (PRIVATE)
# --------------------------------------------------------
# AK_CARB_0_WEIS74
# AK_CARB_1_MILL95, AK_CARB_2_MILL95
# AK_CARB_1_LUEK00, AK_CARB_2_LUEK00
# AK_CARB_1_ROYE93, AK_CARB_2_ROYE93
# AK_BORA_DICK90
# AK_W_MILL95
# AK_PHOS_1_MILL95, AK_PHOS_2_MILL95, AK_PHOS_3_MILL95
# AK_SILI_1_MILL95
# AK_H2S_1_MILL95
# AK_AMMO_1_YAMI95
# ABETA_HF_DIRI79
# AK_HF_PEFR87
# AK_HSO4_DICK90
# A_IONSTRENGTH_SALIN
# ASP_CALC_MUCC83, ASP_ARAG_MUCC83
# A_BTOT_SALIN, A_CATOT_SALIN, A_FTOT_SALIN, A_SO4TOT_SALIN
# ACVT_HSWS_O_HTOT, ACVT_HTOT_O_HFREE, ACVT_HSWS_O_HFREE
# A_RHOSW1_MUNH97, A_RHOSW2_MUNH97
# --------------------------------------
# Parameters for usage within the module
# --------------------------------------
# Gas constant
# ------------
#gasconst_bar_cm3_o_mol_k = 83.14510 # libthdyct
gasconst_bar_cm3_o_mol_k = 83.14472 # Handbook (2007)
# 0 degrees centigrade in Kelvin
# ------------------------------
t_k_zerodegc = 273.15 # Handbook (2007)
#=======================================================================
compute_dissoc_constants <- function (p_temp, p_sal, p_pres, p_pHscale="T")
#=======================================================================
#
# This function computes and returns dissociation constants
# of main chemical elements dissolved in sea water
#
# Input parameters :
# p_temp : temperature in degree Celsius
# p_sal : salinity in psu (Practical Salinity Unit)
# p_pres : pressure in bar
# p_pHscale : pH scale: "T" for Total scale (default), "SWS" for SeaWater, "F" for Free
#
# Output : this function returns a named list of all dissociation constants. Member names are :
#
# K1_DIC : First dissociation constant of carbonic acid (mol/kg) on chosen scale
# K2_DIC : Second dissociation constant of carbonic acid (mol/kg) on chosen scale
# K_BT : Dissociation constant of boric acid (mol/kg) on chosen scale
# K1_PO4 : First dissociation constant of phosphoric acid (mol/kg) on chosen scale
# K2_PO4 : Second dissociation constant of phosphoric acid (mol/kg) on chosen scale
# K3_PO4 : third dissociation constant of phosphoric acid (mol/kg) on chosen scale
# K_Sil : Dissociation constant of sillicic acid (mol/kg) on chosen scale
# K_NH4 : Dissociation constant of ammonium (mol/kg) on chosen scale
# K_H2S : Dissociation constant of hydrogen sulfide (mol/kg) on chosen scale
# K_HSO4 : Dissociation constant of hydrogen sulfate (mol/kg) on free scale
# K_HF : Dissociation constant of hydrogen fluoride (mol/kg) on free scale
# K_H2O : Dissociation constant of water (mol/kg) on chosen scale
#
# Note : K2_Sil, second dissociation constant of sillicic acid is not computed.
{
t_k = p_temp + t_k_zerodegc
dissoc <- list()
if (p_pHscale == "SWS")
{
zcvt_hsws_o_htot = ACVT_HSWS_O_HTOT(t_k, p_sal, p_pres)
dissoc$K1_DIC <- AK_CARB_1_MILL95(t_k, p_sal, p_pres)
dissoc$K2_DIC <- AK_CARB_2_MILL95(t_k, p_sal, p_pres)
dissoc$K_BT <- AK_BORA_DICK90(t_k, p_sal, p_pres) * zcvt_hsws_o_htot
dissoc$K1_PO4 <- AK_PHOS_1_MILL95(t_k, p_sal, p_pres)
dissoc$K2_PO4 <- AK_PHOS_2_MILL95(t_k, p_sal, p_pres)
dissoc$K3_PO4 <- AK_PHOS_3_MILL95(t_k, p_sal, p_pres)
dissoc$K_Sil <- AK_SILI_1_MILL95(t_k, p_sal)
dissoc$K_NH4 <- AK_AMMO_1_YAMI95(t_k, p_sal, p_pres)
dissoc$K_H2S <- AK_H2S_1_MILL95(t_k, p_sal, p_pres)
dissoc$K_HSO4 <- AK_HSO4_DICK90(t_k, p_sal, p_pres)
dissoc$K_HF <- 1. / ABETA_HF_DIRI79(t_k, p_sal, p_pres)
dissoc$K_H2O <- AK_W_MILL95(t_k, p_sal, p_pres)
}
else if (p_pHscale == "F")
{
zcvt_hfree_o_htot = 1. / ACVT_HTOT_O_HFREE(t_k, p_sal, p_pres)
zcvt_hfree_o_hsws = 1. / ACVT_HSWS_O_HFREE(t_k, p_sal, p_pres)
dissoc$K1_DIC <- AK_CARB_1_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K2_DIC <- AK_CARB_2_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_BT <- AK_BORA_DICK90(t_k, p_sal, p_pres) * zcvt_hfree_o_htot
dissoc$K1_PO4 <- AK_PHOS_1_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K2_PO4 <- AK_PHOS_2_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K3_PO4 <- AK_PHOS_3_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_Sil <- AK_SILI_1_MILL95(t_k, p_sal) * zcvt_hfree_o_hsws
dissoc$K_NH4 <- AK_AMMO_1_YAMI95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_H2S <- AK_H2S_1_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
dissoc$K_HSO4 <- AK_HSO4_DICK90(t_k, p_sal, p_pres)
dissoc$K_HF <- 1. / ABETA_HF_DIRI79(t_k, p_sal, p_pres)
dissoc$K_H2O <- AK_W_MILL95(t_k, p_sal, p_pres) * zcvt_hfree_o_hsws
}
else
{
zcvt_htot_o_hsws = 1. / ACVT_HSWS_O_HTOT(t_k, p_sal, p_pres)
zcvt_htot_o_hfree = ACVT_HTOT_O_HFREE(t_k, p_sal, p_pres)
dissoc$K1_DIC <- AK_CARB_1_LUEK00(t_k, p_sal, p_pres)
dissoc$K2_DIC <- AK_CARB_2_LUEK00(t_k, p_sal, p_pres)
dissoc$K_BT <- AK_BORA_DICK90(t_k, p_sal, p_pres)
dissoc$K1_PO4 <- AK_PHOS_1_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K2_PO4 <- AK_PHOS_2_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K3_PO4 <- AK_PHOS_3_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K_Sil <- AK_SILI_1_MILL95(t_k, p_sal) * zcvt_htot_o_hsws
dissoc$K_NH4 <- AK_AMMO_1_YAMI95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K_H2S <- AK_H2S_1_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
dissoc$K_HSO4 <- AK_HSO4_DICK90(t_k, p_sal, p_pres)
dissoc$K_HF <- AK_HF_PEFR87(t_k, p_sal, p_pres) / zcvt_htot_o_hfree
dissoc$K_H2O <- AK_W_MILL95(t_k, p_sal, p_pres) * zcvt_htot_o_hsws
}
return (dissoc)
}
#=======================================================================
AK_CARB_0_WEIS74 <- function (t_k, s)
#=======================================================================
{
# Function calculates K0 in (mol/kg-SW)/atmosphere
# References: Weiss (1979) [(mol/kg-SW)/atm]
# pH scale : N/A
# Note : currently no pressure correction
# ------------------
# Argument variables
# ------------------
# s : salinity
# t_k : temperature in K
zt_k_o_100 = t_k/100.
AK_CARB_0_WEIS74 = exp( -60.2409 + 93.4517/zt_k_o_100
+ 23.3585*log(zt_k_o_100)
+ ( 0.023517 - 0.023656*zt_k_o_100
+ 0.0047036*zt_k_o_100*zt_k_o_100)*s )
return (AK_CARB_0_WEIS74)
}
#=======================================================================
AK_CARB_1_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates first dissociation constant of carbonic acid
# in mol/kg-SW on the SWS pH-scale.
# References: Millero (1995, eq 50 -- ln K1(COM))
# Millero (1982) pressure correction
# pH scale: SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kc1_p0 : ln(K_C1) at p_bar = 0
# zln_kc1_pp : pressure correction for p_bar /= 0
# ln(K_C1) value at p_bar = 0
zsqrts = sqrt(s)
zln_kc1_p0 = ( 2.18867 - 2275.0360/t_k - 1.468591*log(t_k)
+ (-0.138681 - 9.33291/t_k)*zsqrts
+ 0.0726483*s
- 0.00574938*s*zsqrts)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -25.50 - 0.151*zds + 0.1271*zt_degc
zdki = ( -3.08 - 0.578*zds + 0.0877*zt_degc)*1.0E-03
zln_kc1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C1 value
AK_CARB_1_MILL95 = exp( zln_kc1_p0 + zln_kc1_pp )
return (AK_CARB_1_MILL95)
}
#=======================================================================
AK_CARB_2_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates second dissociation constant K1
# in mol/kg-SW on the SWS pH-scale.
# References: Millero (1995, eq 51 -- ln K2(COM))
# Millero (1979) pressure correction
# pH scale: SWS
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kc2_p0 : ln(K_C2) at p_bar = 0
# zln_kc2_pp : pressure correction for p_bar /= 0
# ln(K_C2) value at p_bar = 0
zsqrts = sqrt(s)
zln_kc2_p0 = ( -0.84226 - 3741.1288/t_k - 1.437139*log(t_k)
+ (-0.128417 - 24.41239/t_k)*zsqrts
+ 0.1195308*s
- 0.00912840*s*zsqrts)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -15.82 + 0.321*zds - 0.0219*zt_degc
zdki = ( 1.13 - 0.314*zds - 0.1475*zt_degc)*1.0E-03
zln_kc2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C2 value
AK_CARB_2_MILL95 = exp( zln_kc2_p0 + zln_kc2_pp )
return (AK_CARB_2_MILL95)
}
#=======================================================================
AK_CARB_1_LUEK00 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates first dissociation constant of carbonic acid
# in mol/kg-SW on the Total pH-scale.
# References: Luecker et al. (2000) -- also Handbook (2007)
# Millero (1979) pressure correction
# pH scale: Total
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zds : salinity-34.8
# zlog10_kc1_p0 : log_10(k_C1) at p_bar = 0
# zln_kc1_pp : pressure correction for p_bar /= 0
# log_10(K_C1) value at p_bar = 0
zlog10_kc1_p0 = ( 61.2172 - 3633.86/t_k - 9.67770*log(t_k)
+ s*(0.011555 - s*0.0001152))
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -25.50 - 0.151*zds + 0.1271*zt_degc
zdki = ( -3.08 - 0.578*zds + 0.0877*zt_degc)*1.0E-03
zln_kc1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C1 value
AK_CARB_1_LUEK00 = 10.**zlog10_kc1_p0 * exp(zln_kc1_pp)
return (AK_CARB_1_LUEK00)
}
#=======================================================================
AK_CARB_2_LUEK00 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates second dissociation constant K1
# in mol/kg-SW on the Total pH-scale.
# References: Luecker et al. (2000) -- also Handbook (2007)
# Millero (1979) pressure correction
# pH scale: Total
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zlog10_kc2_p0 : log_10(K_C2) at p_bar = 0
# zln_kc2_pp : pressure correction for p_bar /= 0
# log_10(K_C2) value at p_bar = 0
zlog10_kc2_p0 = (-25.9290 - 471.78/t_k + 3.16967*log(t_k)
+ s*(0.01781 - s*0.0001122))
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -15.82 + 0.321*zds - 0.0219*zt_degc
zdki = ( 1.13 - 0.314*zds - 0.1475*zt_degc)*1.0E-03
zln_kc2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C2 value
AK_CARB_2_LUEK00 = 10.**zlog10_kc2_p0 *exp(zln_kc2_pp)
return (AK_CARB_2_LUEK00)
}
#=======================================================================
AK_CARB_1_ROYE93 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates first dissociation constant of carbonic acid
# in mol/kg-SW on the Total pH-scale.
# References: Roy et al. (1993) -- also Handbook (1994)
# Millero (1979) pressure correction
# pH scale : Total
# Note : converted here from mol/kg-H2O to mol/kg-SW
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zds : salinity-34.8
# zln_kc1_p0 : ln(k_C1) at p_bar = 0
# zln_kc1_pp : pressure correction for p_bar /= 0
# ln(K_C1) value at p_bar = 0
zsqrts = sqrt(s)
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zln_kc1_p0 = (-2307.1266/t_k + 2.83655 - 1.5529413*log(t_k)
+ (-4.0484/t_k - 0.20760841)*zsqrts
+ 0.08468345*s
- 0.00654208*zsqrts*s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -25.50 - 0.151*zds + 0.1271*zt_degc
zdki = ( -3.08 - 0.578*zds + 0.0877*zt_degc)*1.0E-03
zln_kc1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C1 value
AK_CARB_1_ROYE93 = exp(zln_kc1_p0 + zln_kc1_pp) * zcvt_to_kgsw
return (AK_CARB_1_ROYE93)
}
#=======================================================================
AK_CARB_2_ROYE93 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates second dissociation constant K1
# in mol/kg-SW on the Total pH-scale.
# References: Roy et al. (1993) -- also Handbook (1994)
# Millero (1979) pressure correction
# pH scale : Total
# Note : converted here from mol/kg-H2O to mol/kg-SW
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kc2_p0 : ln(K_C2) at p_bar = 0
# zln_kc2_pp : pressure correction for p_bar /= 0
# ln(K_C2) value at p_bar = 0
zsqrts = sqrt(s)
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zln_kc2_p0 = (-3351.6106/t_k - 9.226508 - 0.2005743*log(t_k)
+ ( -23.9722/t_k - 0.106901773)*zsqrts
+ 0.1130822*s
- 0.00846934*zsqrts*s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -15.82 + 0.321*zds - 0.0219*zt_degc
zdki = ( 1.13 - 0.314*zds - 0.1475*zt_degc)*1.0E-03
zln_kc2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_C2 value
AK_CARB_2_ROYE93 = exp(zln_kc2_p0 + zln_kc2_pp) * zcvt_to_kgsw
return (AK_CARB_2_ROYE93)
}
#=======================================================================
AK_BORA_DICK90 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates boric acid dissociation constant KB
# in mol/kg-SW on the total pH-scale.
# References: Dickson (1990, eq. 23) -- also Handbook (2007, eq. 37)
# Millero (1979) pressure correction
# pH scale : total
# ------------------
# Argument variables
# ------------------
# t_k : temperature in Kelvin
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zsqrts : square root of salinity
# zds : salinity-34.8
# zln_kb_p0 : K_b at p_bar = 0
# zln_kb_pp : pressure correction for p_bar /= 0
# ln(K_B) value at p_bar = 0
zsqrts = sqrt(s)
zln_kb_p0 = (( -8966.90
+ zsqrts*( -2890.53
+ zsqrts*( -77.942
+ zsqrts*( 1.728 - 0.0996*zsqrts)))) / t_k
+ 148.0248 + zsqrts*(137.1942 + zsqrts*1.62142)
+ (-24.4344 + zsqrts*(-25.085 - zsqrts*0.2474)) * log(t_k)
+ 0.053105*zsqrts*t_k)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zds = s - 34.8
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -29.48 + 0.295*zds + 0.1622*zt_degc - 0.002608*zt_degc*zt_degc
zdki = (-2.84 + 0.354*zds)*1.0E-03
zln_kb_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_B value
AK_BORA_DICK90 = exp( zln_kb_p0 + zln_kb_pp )
return (AK_BORA_DICK90)
}
#=======================================================================
AK_W_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates water dissociation constant Kw in (mol/kg-SW)^2
# References: Millero (1995) for value at p_bar = 0
# Millero (pers. comm. 1996) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kw_p0 : ln(K_w) at p_bar = 0
# zln_kw_pp : pressure correction for p_bar /= 0
# ln(K_w) value at p_bar = 0
zln_kw_p0 = ( 148.9802
- 13847.26/t_k
- 23.6521*log(t_k)
+ ( -5.977 + 118.67/t_k + 1.0495*log(t_k))*sqrt(s)
- 0.01615*s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -20.02 + 0.1119*zt_degc - 0.1409E-02*zt_degc*zt_degc
zdki = ( -5.13 + 0.0794*zt_degc)*1.0E-03
zln_kw_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_w value
AK_W_MILL95 = exp( zln_kw_p0 + zln_kw_pp )
return (AK_W_MILL95)
}
#=======================================================================
AK_PHOS_1_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the first dissociation constant
# of phosphoric acid (H3PO4) in seawater
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp1_p0 : ln(K_p1) at p_bar = 0
# zln_kp1_pp : pressure correction for p_bar /= 0
# ln(K_P1) for p_bar = 0
zln_kp1_p0 = ( 115.54 - 4576.752/t_k - 18.453*log(t_k)
+ ( 0.69171 - 106.736/t_k)* sqrt(s)
+ (-0.01844 - 0.65643/t_k)*s )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -14.51 + 0.1211*zt_degc - 0.321E-03*zt_degc*zt_degc
zdki = ( -2.67 + 0.0427*zt_degc)*1.0E-03
zln_kp1_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of K_P1
AK_PHOS_1_MILL95 = exp(zln_kp1_p0 + zln_kp1_pp)
return (AK_PHOS_1_MILL95)
}
#=======================================================================
AK_PHOS_2_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the second dissociation constant
# of phosphoric acid (H3PO4) in seawater
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp2_p0 : ln(K_P2) at p_bar = 0
# zln_kp2_pp : pressure correction for p_bar /= 0
# ln(K_P2) for p_bar = 0
zln_kp2_p0 = ( 172.1033
- 8814.715/t_k
- 27.927*log(t_k)
+ ( 1.3566 - 160.340/t_k)*sqrt(s)
+ (-0.05778 + 0.37335/t_k)*s )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -23.12 + 0.1758*zt_degc -2.647E-03*zt_degc*zt_degc
zdki = ( -5.15 + 0.09*zt_degc)*1.0E-03
zln_kp2_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_P2 value
AK_PHOS_2_MILL95 = exp( zln_kp2_p0 + zln_kp2_pp )
return (AK_PHOS_2_MILL95)
}
#=======================================================================
AK_PHOS_3_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the third dissociation constant
# of phosphoric acid (H3PO4) in seawater
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp3_p0 : ln(K_P3) at p_bar = 0
# zln_kp3_pp : pressure correction for p_bar /= 0
# ln(K_P3) for p_bar = 0
zln_kp3_p0 = ( -18.126 - 3070.75/t_k
+ ( 2.81197 + 17.27039/t_k)*sqrt(s)
+ (-0.09984 - 44.99486/t_k)*s )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -26.57 + 0.2020*zt_degc -3.042E-03*zt_degc*zt_degc
zdki = ( -4.08 + 0.0714*zt_degc)*1.0E-03
zln_kp3_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_P3 value
AK_PHOS_3_MILL95 = exp( zln_kp3_p0 + zln_kp3_pp )
return (AK_PHOS_3_MILL95)
}
#=======================================================================
AK_SILI_1_MILL95 <- function (t_k, s)
#=======================================================================
{
# Function returns the first dissociation constant
# of silicic acid (H4SiO4) in seawater
# References: Yao and Millero (1995) cited by Millero (1995)
# pH scale : SWS (according to Dickson et al, 2007)
# Note : No pressure correction available
# Note : converted here from mol/kg-H2O to mol/kg-sw
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# ---------------
# Local variables
# ---------------
# zcvt_to_kgsw: fraction of pure water in 1 kg seawater at salinity s
# zionst : ionic strength [mol/kg-H2O]
# zln_ksi1_p0 : ln(K_Si1) at p_bar = 0
# zln_ksi1_pp : pressure correciotn for p_bar /= 0
# K_Si1 value at p_bar = 0
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zionst = A_IONSTRENGTH_SALIN(s)/zcvt_to_kgsw # mol/kg-H2O ##
zln_ksi1_p0 = ( 117.40 - 8904.2/t_k - 19.334 * log(t_k)
+ ( 3.5913 - 458.79/t_k) * sqrt(zionst)
+ (-1.5998 + 188.74/t_k) * zionst
+ (0.07871 - 12.1652/t_k) * zionst*zionst )
# Pressure correction : currently none
zln_ksi1_pp = 0.
# Final value
AK_SILI_1_MILL95 = exp( zln_ksi1_p0 + zln_ksi1_pp ) * zcvt_to_kgsw
return (AK_SILI_1_MILL95)
}
#=======================================================================
AK_H2S_1_MILL95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the dissociation constant of hydrogen sulfide in sea-water
# References: Millero et al. (1988) (cited by Millero (1995)
# Millero (1995) for pressure correction
# pH scale : - SWS (according to Yao and Millero, 1995, p. 82: "refitted if necessary")
# - Total (according to Lewis and Wallace, 1998)
# Note : we stick to SWS here for the time being
# Note : the fits from Millero (1995) and Yao and Millero (1995)
# derive from Millero et al. (1988), with all the coefficients
# multiplied by -ln(10)
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zt_degc : temperature in degrees Celsius
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kh2s_p0 : ln(K_H2S) at p_bar = 0
# zln_kh2s_pp : pressure correction for p_bar /= 0
# K_H2S value at p_bar = 0
# ------------------------
zln_kh2s_p0 = ( 225.838
- 13275.3/t_k
- 34.6435 * log(t_k)
+ 0.3449*sqrt(s)
- 0.0274*s )
# Pressure correction
# -------------------
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -14.80 + zt_degc*(0.0020 - zt_degc*0.400E-03)
zdki = ( 2.89 + zt_degc*0.054)*1.0E-03
zln_kh2s_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_H2S value
# -----------------
AK_H2S_1_MILL95 = exp( zln_kh2s_p0 + zln_kh2s_pp )
return (AK_H2S_1_MILL95)
}
#=======================================================================
AK_AMMO_1_YAMI95 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the dissociation constant
# of ammonium in sea-water [mol/kg-SW]
# References: Yao and Millero (1995)
# Millero (1995) for pressure correction
# pH scale : SWS
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zt_degc : temperature in degrees Celsius
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_knh4_p0 : ln(K_NH4) at p_bar = 0
# zln_knh4_pp : pressure correction for p_bar /= 0
# K_NH4 value at p_bar = 0
# ------------------------
zln_knh4_p0 = ( -0.25444 - 6285.33/t_k + 0.0001635*t_k
+ ( 0.46532 - 123.7184/t_k) * sqrt(s)
+ (-0.01992 + 3.17556/t_k) * s )
# Pressure correction
# -------------------
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -26.43 + zt_degc*(0.0889 - zt_degc*0.905E-03)
zdki = ( -5.03 + zt_degc*0.0814)*1.0E-03
zln_knh4_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final K_NH4 value
# -----------------
AK_AMMO_1_YAMI95 = exp( zln_knh4_p0 + zln_knh4_pp )
return (AK_AMMO_1_YAMI95)
}
#=======================================================================
ACVT_KGH2O_O_KGSW <- function (s)
#=======================================================================
{
# Function returns the mass of pure water in one kg of seawater
# of salinity s
# References: "libthdyct" -- derived by Munhoven (1997) from data by Millero (1982)
# "Handbook (2007)" -- Handbook (2007)
# pH scale: N/A
#ACVT_KGH2O_O_KGSW = 1. - 0.0010049*s # libthdyct
ACVT_KGH2O_O_KGSW = 1. - 0.001005*s # Handbook (2007)
return (ACVT_KGH2O_O_KGSW)
}
#=======================================================================
A_IONSTRENGTH_SALIN <- function (s)
#=======================================================================
{
# Function calculates ionic strength in mol/kg-SW, for given salinity.
# References: "libthdyct" -- derived by Munhoven (1997) from data by Millero (1982)
# "Handbook (2007)" -- Handbook (2007)
# pH scale: N/A
#A_IONSTRENGTH_SALIN = (0.019920*s) # libthdyct
A_IONSTRENGTH_SALIN = (0.019924*s) # Handbook (2007)
return (A_IONSTRENGTH_SALIN)
}
#=======================================================================
ABETA_HF_DIRI79 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates association constant \beta_{HF} [(mol/kg-SW)^{-1}]
# in (mol/kg-SW)^{-1}, where
# \beta_{HF} = \frac{ [HF] }{ [H^{+}] [F^{-}] }
# References: Dickson and Riley (1979)
# Millero (1995) for pressure correction
# pH scale : free
# Note : converted here from mol/kg-H2O to mol/kg-SW
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zionst : ionic strength [mol/kg-H2O]
# zcvt_to_kgsw : mass of pure water in 1kg of seawater as a fct. of salinity
# zln_bhf_p0 : \beta_HF at p_bar = 0
# zln_khf_pp : pressure correction for k_HF = 1/\beta_HF at p_bar /= 0
# \beta_HF at p_bar = 0
# ---------------------
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zionst = A_IONSTRENGTH_SALIN(s)/zcvt_to_kgsw
zln_bhf_p0 = -1590.2/t_k + 12.641 - 1.525*sqrt(zionst)
# Pressure correction
# -------------------
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -9.78 + zt_degc*(-0.0090 - zt_degc*0.942E-03)
zdki = ( -3.91 + zt_degc*0.054)*1.0E-03
zln_khf_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final \beta_HF value
# --------------------
# notice that ln(k_HF(P)) = ln(k_HF(0)) + zln_khf_pp
# <=> -ln(\beta_HF(P)) = -ln(\beta_HF(0)) + zln_khf_pp
# <=> ln(\beta_HF(P)) = ln(\beta_HF(0)) - zln_khf_pp
ABETA_HF_DIRI79 = exp(zln_bhf_p0 - zln_khf_pp ) / zcvt_to_kgsw
return (ABETA_HF_DIRI79)
}
#=======================================================================
AK_HF_PEFR87 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function calculates dissociation constant for hydrogen fluoride
# in mol/kg-SW
# References: Perez and Fraga (1987)
# Millero (1995) for pressure correction
# pH scale : Total (according to Handbook, 2007)
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_khf_p0 : ln(K_HF) at p_bar = 0
# zln_khf_pp : pressure correction for p_bar /= 0
# ln(K_HF) at p_bar = 0
zln_khf_p0 = 874./t_k - 9.68 + 0.111*sqrt(s)
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -9.78 + zt_degc*(-0.0090 - zt_degc*0.942E-03)
zdki = ( -3.91 + zt_degc*0.054)*1.0E-03
zln_khf_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of K_HF
AK_HF_PEFR87 = exp( zln_khf_p0 + zln_khf_pp )
return (AK_HF_PEFR87)
}
#=======================================================================
AK_HSO4_DICK90 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the dissociation constant of hydrogen sulfate (bisulfate)
# References: Dickson (1990) -- also Handbook (2007)
# Millero (1995) for pressure correction
# pH scale : free
# Note : converted here from mol/kg-H2O to mol/kg-SW
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zionst : ionic strength in mol/-kg-H2O
# zsqrti : square root og ion strength
# zcvt_to_kgsw : mass of pure water in 1kg of seawater as a fct. of salinity
# zln_khso4_p0 : K_HSO4 at p_bar = 0
# zln_khso4_pp : pressure correction for p_bar /= 0
# ln(K_HSO4) at p_bar = 0
zcvt_to_kgsw = ACVT_KGH2O_O_KGSW(s)
zionst = A_IONSTRENGTH_SALIN(s)/zcvt_to_kgsw
zsqrti = sqrt(zionst)
zln_khso4_p0 = ( -4276.1/t_k + 141.328 - 23.093*log(t_k)
+ (-13856./t_k + 324.57 - 47.986*log(t_k)) * zsqrti
+ ( 35474./t_k - 771.54 + 114.723*log(t_k)) * zionst
- ( 2698./t_k)*zsqrti * zionst
+ ( 1776./t_k)*zionst*zionst )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -18.03 + zt_degc*(0.0466 + zt_degc*0.316E-03)
zdki = ( -4.53 + zt_degc*0.0900)*1.0E-03
zln_khso4_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# ln(K_HSO4) at p_bar = 0
AK_HSO4_DICK90 = zcvt_to_kgsw * exp( zln_khso4_p0 + zln_khso4_pp )
return (AK_HSO4_DICK90)
}
#=======================================================================
ASP_CALC_MUCC83 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns stoechiometric solubility product
# of calcite in seawater
# References: Mucci (1983)
# Millero (1995) for pressure correction
# pH scale : N/A
# Units : (mol/kg-SW)^2
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zsqrts : square root of salinity
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp1_p0 : ln(K_p1) at p_bar = 0
# zln_kp1_pp : pressure correction for p_bar /= 0
zsqrts = sqrt(s)
# log10(Ksp_Calc) for p_bar = 0
zlog10_kspcalc_p0 = ( -171.9065 - 0.077993*t_k
+ 2839.319/t_k + 71.595*log10(t_k)
+ (-0.77712 + 0.0028426*t_k + 178.34/t_k)*zsqrts
- 0.07711*s
+ 0.0041249*s*zsqrts )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -48.76 + 0.5304*zt_degc
zdki = (-11.76 + 0.3692*zt_degc)*1.0E-03
zln_kspcalc_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of Ksp_Calc
ASP_CALC_MUCC83 = 10.**(zlog10_kspcalc_p0) * exp(zln_kspcalc_pp)
return (ASP_CALC_MUCC83)
}
#=======================================================================
ASP_ARAG_MUCC83 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns stoechiometric solubility product
# of aragonite in seawater
# References: Mucci (1983)
# Millero (1979) for pressure correction
# pH scale : N/A
# Units : (mol/kg-SW)^2
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zrt : R*t_k, R in bar*cm3/(mol*K)
# zsqrts : square root of salinity
# zt_degc : temperature in degrees Celsius
# zdvi : volume change for ionization
# zdki : compressibility change for ionization
# zln_kp1_p0 : ln(K_p1) at p_bar = 0
# zln_kp1_pp : pressure correction for p_bar /= 0
zsqrts = sqrt(s)
# log10(Ksp_Arag) for p_bar = 0
zlog10_ksparag_p0 = ( -171.945 - 0.077993*t_k
+ 2903.293/t_k + 71.595*log10(t_k)
+ (-0.068393 + 0.0017276*t_k + 88.135/t_k)*zsqrts
- 0.10018*s
+ 0.0059415*s*zsqrts )
# Pressure correction
zt_degc = t_k - t_k_zerodegc
zrt = gasconst_bar_cm3_o_mol_k * t_k
zdvi = -48.76 + 0.5304*zt_degc + 2.8
zdki = (-11.76 + 0.3692*zt_degc)*1.0E-03
zln_ksparag_pp = (-zdvi + zdki*p_bar/2.)*p_bar/zrt
# Final value of Ksp_Arag
ASP_ARAG_MUCC83 = 10.**(zlog10_ksparag_p0) * exp(zln_ksparag_pp)
return (ASP_ARAG_MUCC83)
}
#=======================================================================
A_BTOT_SALIN <- function (s)
#=======================================================================
{
# Function returns total borate concentration in mol/kg-SW
# given the salinity of a sample
# References: Uppström (1974), cited by Dickson et al. (2007, chapter 5, p 10)
# Millero (1982) cited in Millero (1995)
# pH scale : N/A
A_BTOT_SALIN = 0.000416*(s/35.)
return (A_BTOT_SALIN)
}
#=======================================================================
A_CATOT_SALIN <- function (s)
#=======================================================================
{
# Function returns total calcium concentration in mol/kg-SW
# given the salinity of a sample
# References: Culkin and Cox (1966)
# pH scale : N/A
# (g Ca/kg)/Cl_permil values
# Culkin (1967): 0.0213
# Culkin and Cox (DSR 13 1966): 0.02126 +/- 0.00004 (stdev)
# Riley and Tongudai (Chem Geol 2 1967): 0.02128 +/- 0.00006 (stdev)
# Handbook (2007): 0.02127
# with reference to Riley and Tongudai (1967) (???)
# Here:
# (g Ca/kg)/Cl_permil = 0.02127
# (g Ca)/(mol Ca) = 40.078
# Cl_permil = S/1.80655
# mol Ca/kg = (0.02127/40.078) * (35/1.80655)
#A_CATOT_SALIN = 0.010282*(s/35.)
A_CATOT_SALIN = (0.02127/40.078) * (s/1.80655)
return (A_CATOT_SALIN)
}
#=======================================================================
A_FTOT_SALIN <- function (s)
#=======================================================================
{
# Function returns total calcium concentration in mol/kg-SW
# given the salinity of a sample
# References: Culkin (1965) (???)
# pH scale : N/A
A_FTOT_SALIN = 0.000068*(s/35.)
return (A_FTOT_SALIN)
}
#=======================================================================
A_SO4TOT_SALIN <- function(s)
#=======================================================================
{
# Function returns total sulfate concentration in mol/kg-SW
# given the salinity of a sample
# References: Morris, A.W. and Riley, J.P. (1966) quoted in Handbook (2007)
# pH scale : N/A
#A_SO4TOT_SALIN = 0.028234*(s/35.) # in libthdyct and Thesis
#A_SO4TOT_SALIN = 0.02824*(s/35.) # Handbook (2007, chap 6, p 10, tab 2, col 3)
A_SO4TOT_SALIN = (0.1400/96.062)*(s/1.80655) # Handbook (2007, chap 6, p 10)
return (A_SO4TOT_SALIN)
}
#=======================================================================
ACVT_HSWS_O_HTOT <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the ratio H_SWS/H_Tot as a function of salinity s
# Reference: Munhoven
# pH scale: all
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zso4_tot: total sulfate concentration in mol/kg-SW
# zf_tot : total fluoride concentration in mol/kg-SW
#-----------------------------------------------------------------------
zso4_tot = A_SO4TOT_SALIN(s)
zf_tot = A_FTOT_SALIN(s)
ACVT_HSWS_O_HTOT = (1. + (zf_tot*ABETA_HF_DIRI79(t_k, s, p_bar))
/(1. + zso4_tot/AK_HSO4_DICK90(t_k,s, p_bar)))
return (ACVT_HSWS_O_HTOT)
}
#=======================================================================
ACVT_HTOT_O_HFREE <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the ratio H_Tot/H_free as a function of salinity s
# Reference: Munhoven
# pH scale: N/A
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zso4_tot: total sulfate concentration in mol/kg-SW
#-----------------------------------------------------------------------
zso4_tot = A_SO4TOT_SALIN(s)
ACVT_HTOT_O_HFREE = 1. + zso4_tot/AK_HSO4_DICK90(t_k,s, p_bar)
return (ACVT_HTOT_O_HFREE)
}
#=======================================================================
ACVT_HSWS_O_HFREE <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns the ratio H_SWS/H_free as a function
# of salinity s
# Reference: Munhoven
# pH scale: N/A
# ------------------
# Argument variables
# ------------------
# t_k : temperature in K
# s : salinity
# p_bar : applied pressure in bar
# ---------------
# Local variables
# ---------------
# zso4_tot: total sulfate concentration in mol/kg-SW
# zf_tot : total fluoride concentration in mol/kg-SW
#-----------------------------------------------------------------------
zso4_tot = A_SO4TOT_SALIN(s)
zf_tot = A_FTOT_SALIN(s)
ACVT_HSWS_O_HFREE = ( 1. + zf_tot*ABETA_HF_DIRI79(t_k, s, p_bar)
+ zso4_tot/AK_HSO4_DICK90(t_k,s, p_bar) )
return (ACVT_HSWS_O_HFREE)
}
#=======================================================================
A_RHOSW1_MUNH97 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns first order approximation of \rho in (kg-SW)/(m^3-SW)
# References: Munhoven (1997)
# after EOS80 (UNESCO, 1981, 1983)
# ------------------
# Argument variables
# ------------------
# s : salinity
# tk : temperature in K
# p_bar : depth in m
# ---------------
# Local variables
# ---------------
s0 = 35.5
t_k0 = 285.16
p_bar0 = 300.0
A_RHOSW1_MUNH97 = ( 1039.9044 + 0.77629393*(s-s0)
- 0.19692738*(t_k-t_k0) )
return (A_RHOSW1_MUNH97)
}
#=======================================================================
A_RHOSW2_MUNH97 <- function (t_k, s, p_bar)
#=======================================================================
{
# Function returns second order approximation of \rho in (kg-SW)/(m^3-SW)
# References: Munhoven (1997)
# after EOS80 (UNESCO, 1981, 1983)
# ------------------
# Argument variables
# ------------------
# s : salinity
# tk : temperature in K
# p_bar : depth in m
# ---------------
# Local variables
# ---------------
s0 = 35.5
t_k0 = 285.16
p_bar0 = 300.0
A_RHOSW2_MUNH97 =( 1040.0145
+ 0.77629393*(s-s0)
- 0.25013591*(t_k-t_k0)
+ 4.2026266E-02*(p_bar-p_bar0)
- 4.7473116E-03*(t_k-t_k0)*(t_k-t_k0)
- 4.7974224E-06*(p_bar-p_bar0)*(p_bar-p_bar0)
- 2.1404592E-04*(t_k-t_k0)*(p_bar-p_bar0) )
return (A_RHOSW2_MUNH97)
}
#=======================================================================
CHECKCONSTANTS <- function()
#=======================================================================
{
# ---------------
# Local variables
# ---------------
# s : salinity
# tk : temperature in K
# p_bar : applied pressure in bar
print ('Checking constant values generated from MOD_CHEMCONST')
print ('')
print (' % indicates checking against the Handbook (1994);')
print (' $ indicates checking against the Lewis and Wallace (1998);')
print (' * indicates checking against the Handbook (2007);')
print (' target values are quoted in brackets')
print ('')
print ('')
print (' For S = 35, P = 0 and T/K = 298.15:')
print ('')
s = 35.
p_bar = 0.
t_k = 298.15
zkc0 = AK_CARB_0_WEIS74(t_k, s)
print ('')
print ('K_0 -- Weiss (1974)')
print ('===================')
print ('')
print (c(' K_0 :', zkc0))
print (c(' ln(K_0) :', log(zkc0)))
print (c(' pK_0 :', -log10(zkc0)))
print ('')
zkhso4 = AK_HSO4_DICK90(t_k, s, p_bar)
print ('')
print ('K_HSO4 -- Dickson (1990) -- pH_free')
print ('===================================')
print ('')
print (c(' K_HSO4 :', zkhso4))
print (c(' ln(K_HSO4) :', log(zkhso4)))
print (c(' pK_HSO4 :', -log10(zkhso4)))
print ('')
zkb = AK_BORA_DICK90(t_k, s, p_bar)
print ('')
print ('K_b -- Dickson (1990) -- pH_tot')
print ('===============================')
print ('')
print (c(' K_b :', zkb))
print (c(' ln(K_b) :', log(zkb)))
print (c(' pK_b :', -log10(zkb)))
print ('')
zkc1 = AK_CARB_1_LUEK00(t_k, s, p_bar)
print ('')
print ('K_1 -- Luecker et al (2000) -- pH_tot')
print ('=====================================')
print ('')
print (c(' K_1 :', zkc1))
print (c(' ln(K_1) :', log(zkc1)))
print (c(' pK_1 :', -log10(zkc1)))
print ('')
zkc2 = AK_CARB_2_LUEK00(t_k, s, p_bar)
print ('')
print ('K_2 -- Luecker et al (2000) -- pH_tot')
print ('=====================================')
print ('')
print (c(' K_2 :', zkc2))
print (c(' ln(K_2) :', log(zkc2)))
print (c(' pK_2 :', -log10(zkc2)))
print ('')
zkc1 = AK_CARB_1_ROYE93(t_k, s, p_bar)
print ('')
print ('K_1 -- Roy et al (1993) -- pH_tot')
print ('=================================')
print ('')
print (c(' K_1 :', zkc1))
print (c(' ln(K_1) :', log(zkc1)))
print (c(' pK_1 :', -log10(zkc1)))
print ('')
zkc2 = AK_CARB_2_ROYE93(t_k, s, p_bar)
print ('')
print ('K_2 -- Roy et al (1993) -- pH_tot')
print ('=================================')
print ('')
print (c(' K_2 :', zkc2))
print (c(' ln(K_2) :', log(zkc2)))
print (c(' pK_2 :', -log10(zkc2)))
print ('')
zkhf = AK_HF_PEFR87(t_k, s, p_bar)
print ('')
print ('K_HF -- Perez and Fraga (1987) -- pH_tot')
print ('========================================')
print ('')
print (c(' K_HF :', zkhf))
print (c(' ln(K_HF) :', log(zkhf)))
print (c(' pK_HF :', -log10(zkhf)))
print ('')
zkp1 = AK_PHOS_1_MILL95(t_k, s, p_bar)
print ('')
print ('K_P1 -- Millero (1995) -- pH_SWS')
print ('================================')
print ('')
print (c(' K_P1 :', zkp1))
print (c(' ln(K_P1) :', log(zkp1)))
print (c(' pK_1 :', -log10(zkp1)))
print ('')
zkp2 = AK_PHOS_2_MILL95(t_k, s, p_bar)
print ('')
print ('K_P2 -- Millero (1995) -- pH_SWS')
print ('================================')
print ('')
print (c(' K_2 :', zkp2))
print (c(' ln(K_P2) :', log(zkp2)))
print (c(' pK_2 :', -log10(zkp2)))
print ('')
zkp3 = AK_PHOS_3_MILL95(t_k, s, p_bar)
print ('')
print ('K_P3 -- Millero (1995) -- pH_SWS')
print ('================================')
print ('')
print (c(' K_P3 :', zkp3))
print (c(' ln(K_P3) :', log(zkp3)))
print (c(' pK_P3 :', -log10(zkp3)))
print ('')
zksi1 = AK_SILI_1_MILL95(t_k, s)
print ('')
print ('K_Si1 -- Millero (1995) -- pH_SWS')
print ('=================================')
print ('')
print (c(' K_Si1 :', zksi1))
print (c(' ln(K_Si1) :', log(zksi1)))
print (c(' pK_Si1 :', -log10(zksi1)))
print ('')
zkw = AK_W_MILL95(t_k, s, p_bar)
print ('')
print ('K_w -- Millero (1995) -- pH_SWS')
print ('===============================')
print ('')
print (c(' K_w :', zkw))
print (c(' ln(K_w) :', log(zkw)))
print (c(' pK_w :', -log10(zkw)))
print ('')
zkh2s = AK_H2S_1_MILL95(t_k, s, p_bar)
print ('')
print ('K_H2S -- Millero (1995) -- pH_SWS')
print ('=================================')
print ('')
print (c(' K_H2S :', zkh2s))
print (c(' ln(K_H2S) :', log(zkh2s)))
print (c(' pK_H2S :', -log10(zkh2s)))
print ('')
zknh4 = AK_AMMO_1_YAMI95(t_k, s, p_bar)
print ('')
print ('K_NH4 -- Yao and Millero (1995) -- pH_SWS')
print ('=========================================')
print ('')
print (c(' K_NH4 :', zknh4))
print (c(' ln(K_NH4) :', log(zknh4)))
print (c(' pK_NH4 :', -log10(zknh4)))
print ('')
return()
}
Try the SolveSAPHE package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.