From 96620e5eb1cac4beaac82b8b8cb1a61b07507323 Mon Sep 17 00:00:00 2001 From: "Matthew P. Humphreys" Date: Mon, 10 Jun 2024 17:04:06 +0200 Subject: [PATCH] More updates --- PyCO2SYS/__init__.py | 7 +- PyCO2SYS/convert.py | 82 +++++++------ PyCO2SYS/engine/system.py | 145 +++++++++++++++------- PyCO2SYS/equilibria/__init__.py | 2 - PyCO2SYS/equilibria/p1atm.py | 210 ++++++++++++++++++++++++-------- newclass_v2.py | 21 +++- 6 files changed, 330 insertions(+), 137 deletions(-) diff --git a/PyCO2SYS/__init__.py b/PyCO2SYS/__init__.py index c2727d8..4514d6d 100644 --- a/PyCO2SYS/__init__.py +++ b/PyCO2SYS/__init__.py @@ -20,9 +20,9 @@ # bio, # buffers, # constants, - # convert, + convert, engine, - # equilibria, + equilibria, # gas, meta, # minimal, @@ -39,11 +39,14 @@ # Aliases for top-level access from .engine import CO2SYS, system + # from .engine.nd import CO2SYS as sys # from .engine.nd import assemble from .engine.system import CO2System + # from .api import CO2SYS_wrap, CO2SYS_MATLABv3 from .meta import hello # because history + # from .solve.get import speciation # from .api.ezio import ezio # from .uncertainty import all_OEDG18 as uncertainty_OEDG18 diff --git a/PyCO2SYS/convert.py b/PyCO2SYS/convert.py index 4bb1e13..66b0e46 100644 --- a/PyCO2SYS/convert.py +++ b/PyCO2SYS/convert.py @@ -3,7 +3,7 @@ """Convert units and calculate conversion factors.""" import copy -from autograd import numpy as np +from jax import numpy as np from . import constants from .equilibria import pressured @@ -58,44 +58,46 @@ def bar_to_decibar(Pbar): return Pbar * 10.0 -def pH_free_to_total(totals, k_constants): - """Free to Total pH scale conversion factor.""" - return 1.0 + totals["TSO4"] / k_constants["KSO4"] +def pH_free_to_total(total_sulfate, k_HSO4_free): + """Free to total pH scale conversion factor.""" + return 1.0 + 1e-6 * total_sulfate / k_HSO4_free -def pH_free_to_sws(totals, k_constants): - """Free to Seawater pH scale conversion factor.""" - return 1.0 + totals["TSO4"] / k_constants["KSO4"] + totals["TF"] / k_constants["KF"] +def pH_free_to_sws(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free): + """Free to seawater pH scale conversion factor.""" + return 1.0 + 1e-6 * total_sulfate / k_HSO4_free + 1e-6 * total_fluoride / k_HF_free -def pH_sws_to_free(totals, k_constants): - """Seawater to Free pH scale conversion factor.""" - return 1.0 / pH_free_to_sws(totals, k_constants) +def pH_sws_to_free(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free): + """Seawater to free pH scale conversion factor.""" + return 1.0 / pH_free_to_sws(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free) -def pH_sws_to_total(totals, k_constants): - """Seawater to Total pH scale conversion factor.""" - return pH_sws_to_free(totals, k_constants) * pH_free_to_total(totals, k_constants) +def pH_sws_to_total(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free): + """Seawater to total pH scale conversion factor.""" + return pH_sws_to_free( + total_fluoride, total_sulfate, k_HF_free, k_HSO4_free + ) * pH_free_to_total(total_sulfate, k_HSO4_free) -def pH_total_to_free(totals, k_constants): - """Total to Free pH scale conversion factor.""" - return 1.0 / pH_free_to_total(totals, k_constants) +def pH_total_to_free(total_sulfate, k_HSO4_free): + """Total to free pH scale conversion factor.""" + return 1.0 / pH_free_to_total(total_sulfate, k_HSO4_free) -def pH_total_to_sws(totals, k_constants): - """Total to Seawater pH scale conversion factor.""" - return 1.0 / pH_sws_to_total(totals, k_constants) +def pH_total_to_sws(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free): + """Total to seawater pH scale conversion factor.""" + return 1.0 / pH_sws_to_total(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free) -def pH_sws_to_nbs(totals, k_constants): +def pH_sws_to_nbs(fH): """Seawater to NBS pH scale conversion factor.""" - return k_constants["fH"] + return fH -def pH_nbs_to_sws(totals, k_constants): +def pH_nbs_to_sws(fH): """NBS to Seawater pH scale conversion factor.""" - return 1.0 / pH_sws_to_nbs(totals, k_constants) + return 1.0 / pH_sws_to_nbs(fH) def pH_total_to_nbs(totals, k_constants): @@ -118,22 +120,24 @@ def pH_nbs_to_free(totals, k_constants): return 1.0 / pH_free_to_nbs(totals, k_constants) -def fH_PTBO87(TempK, Sal): +def fH_PTBO87(temperature, salinity): """fH following PTBO87.""" # === CO2SYS.m comments: ======= # Peng et al, Tellus 39B:439-458, 1987: # They reference the GEOSECS report, but round the value # given there off so that it is about .008 (1#) lower. It # doesn't agree with the check value they give on p. 456. - return 1.29 - 0.00204 * TempK + (0.00046 - 0.00000148 * TempK) * Sal**2 + TempK = temperature + 273.15 + return 1.29 - 0.00204 * TempK + (0.00046 - 0.00000148 * TempK) * salinity**2 -def fH_TWB82(TempK, Sal): +def fH_TWB82(temperature, salinity): """fH following TWB82.""" # === CO2SYS.m comments: ======= # Takahashi et al, Chapter 3 in GEOSECS Pacific Expedition, # v. 3, 1982 (p. 80). - return 1.2948 - 0.002036 * TempK + (0.0004607 - 0.000001475 * TempK) * Sal**2 + TempK = temperature + 273.15 + return 1.2948 - 0.002036 * TempK + (0.0004607 - 0.000001475 * TempK) * salinity**2 def pH_to_all_scales(pH, pH_scale, totals, k_constants): @@ -177,7 +181,9 @@ def get_pHfactor_from_SWS(TempK, Sal, totals, k_constants, pHScale, WhichKs): ) # Total pHfactor = np.where(pHScale == 2, 1.0, pHfactor) # Seawater (SWS) pHfactor = np.where( - pHScale == 3, pH_sws_to_free(totals, k_constants), pHfactor + pHScale == 3, + pH_sws_to_free(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free), + pHfactor, ) # Free pHfactor = np.where( pHScale == 4, pH_sws_to_nbs(totals, k_constants), pHfactor @@ -197,7 +203,9 @@ def get_pHfactor_to_Free(TempK, Sal, totals, k_constants, pHScale, WhichKs): pHScale == 1, pH_total_to_free(totals, k_constants), pHfactor ) # Total pHfactor = np.where( - pHScale == 2, pH_sws_to_free(totals, k_constants), pHfactor + pHScale == 2, + pH_sws_to_free(total_fluoride, total_sulfate, k_HF_free, k_HSO4_free), + pHfactor, ) # Seawater (SWS) pHfactor = np.where(pHScale == 3, 1.0, pHfactor) # Free pHfactor = np.where( @@ -238,9 +246,11 @@ def _flattenfirst(args, dtype): npts = np.max(arglengths) return ( [ - np.full(npts, arg, dtype=dtype) - if np.size(arg) == 1 - else arg.ravel().astype(dtype) + ( + np.full(npts, arg, dtype=dtype) + if np.size(arg) == 1 + else arg.ravel().astype(dtype) + ) for arg in args ], npts, @@ -255,9 +265,11 @@ def _flattenafter(args, npts, dtype): ), "Inputs must all be the same length as each other or of length 1." # Make vectors of all inputs return [ - np.full(npts, arg, dtype=dtype) - if np.size(arg) == 1 - else arg.ravel().astype(dtype) + ( + np.full(npts, arg, dtype=dtype) + if np.size(arg) == 1 + else arg.ravel().astype(dtype) + ) for arg in args ] diff --git a/PyCO2SYS/engine/system.py b/PyCO2SYS/engine/system.py index e0c393c..cd1a788 100644 --- a/PyCO2SYS/engine/system.py +++ b/PyCO2SYS/engine/system.py @@ -4,16 +4,33 @@ import networkx as nx from jax import numpy as np from matplotlib import pyplot as plt -from . import salts +from . import constants, convert, equilibria, salts -# Define functions for calculations that depend neither on icase nor opts +# Define functions for calculations that depend neither on icase nor opts: get_funcs = { + # Total salt contents "ionic_strength": salts.ionic_strength_DOE94, "total_fluoride": salts.total_fluoride_R65, "total_sulfate": salts.total_sulfate_MR66, + # Equilibrium constants at 1 atm and on reported pH scale + "k_H2S_total_1atm": equilibria.p1atm.k_H2S_total_YM95, + # pH scale conversion factors at 1 atm + "free_to_sws_1atm": lambda total_fluoride, total_sulfate, k_HF_free_1atm, k_HSO4_free_1atm: convert.pH_free_to_sws( + total_fluoride, total_sulfate, k_HF_free_1atm, k_HSO4_free_1atm + ), + "nbs_to_sws": convert.pH_nbs_to_sws, + "total_to_sws_1atm": lambda total_fluoride, total_sulfate, k_HF_free_1atm, k_HSO4_free_1atm: convert.pH_total_to_sws( + total_fluoride, total_sulfate, k_HF_free_1atm, k_HSO4_free_1atm + ), + # Equilibrium constants at 1 atm and on the seawater pH scale + "k_H2S_sws_1atm": lambda k_H2S_total_1atm, total_to_sws_1atm: ( + k_H2S_total_1atm * total_to_sws_1atm + ), + # Pressure correction factors for equilibrium constants + # Equilibrium constants at pressure and on the seawater pH scale } -# Define functions for calculations that depend on icase +# Define functions for calculations that depend on icase: get_funcs_core = {} get_funcs_core[0] = {} # get_funcs_core[5] = { @@ -27,8 +44,22 @@ # "CO3": CO3_from_dic_pH, # } -# Define functions for calculations that depend on opts +# Define functions for calculations that depend on opts: +# (unlike in previous versions, each opt may only affect one parameter) get_funcs_opts = {} +get_funcs_opts["opt_fH"] = { + 1: dict(fH=convert.fH_TWB82), + 2: dict(fH=convert.fH_PTBO87), +} +get_funcs_opts["opt_k_HF"] = { + 1: dict(k_HF_free_1atm=equilibria.p1atm.k_HF_free_DR79), + 2: dict(k_HF_free_1atm=equilibria.p1atm.k_HF_free_PF87), +} +get_funcs_opts["opt_k_HSO4"] = { + 1: dict(k_HSO4_free_1atm=equilibria.p1atm.k_HSO4_free_D90a), + 2: dict(k_HSO4_free_1atm=equilibria.p1atm.k_HSO4_free_KRCB77), + 3: dict(k_HSO4_free_1atm=equilibria.p1atm.k_HSO4_free_WM13), +} get_funcs_opts["opt_total_borate"] = { 1: dict(total_borate=salts.total_borate_U74), 2: dict(total_borate=salts.total_borate_LKB10), @@ -42,32 +73,34 @@ # Automatically set up graph for calculations that depend neither on icase nor opts # based on the function names and signatures in get_funcs -links = nx.DiGraph() +graph = nx.DiGraph() for k, func in get_funcs.items(): - func_args = func.__code__.co_varnames + fcode = func.__code__ + func_args = fcode.co_varnames[: fcode.co_argcount] for f in func_args: - links.add_edge(f, k) + graph.add_edge(f, k) # Automatically set up graph for each icase based on the function names and signatures # in get_funcs_core -links_core = {} +graph_core = {} for icase, funcs in get_funcs_core.items(): - links_core[icase] = nx.DiGraph() + graph_core[icase] = nx.DiGraph() for t, func in get_funcs_core[icase].items(): for f in func.__name__.split("_")[2:]: - links_core[icase].add_edge(f, t) + graph_core[icase].add_edge(f, t) # Automatically set up graph for each opt based on the function names and signatures in # get_funcs_opts -links_opts = {} +graph_opts = {} for o, opts in get_funcs_opts.items(): - links_opts[o] = {} + graph_opts[o] = {} for opt, funcs in opts.items(): - links_opts[o][opt] = nx.DiGraph() + graph_opts[o][opt] = nx.DiGraph() for k, func in funcs.items(): - func_args = func.__code__.co_varnames + fcode = func.__code__ + func_args = fcode.co_varnames[: fcode.co_argcount] for f in func_args: - links_opts[o][opt].add_edge(f, k) + graph_opts[o][opt].add_edge(f, k) parameters_core = [ "alkalinity", @@ -94,7 +127,9 @@ } default_opts = { - # "opt_k_carbonic": 10, + "opt_fH": 1, + "opt_k_HF": 1, + "opt_k_HSO4": 1, "opt_total_borate": 1, "opt_Ca": 1, } @@ -123,11 +158,11 @@ def __init__(self, values=None, opts=None, use_default_values=True): ), "{} is not allowed for {}!".format(v, k) self.opts.update(opts) # Assemble graph and functions - self.links = nx.compose(links, links_core[self.icase]) + self.graph = nx.compose(graph, graph_core[self.icase]) self.get_funcs = get_funcs.copy() self.get_funcs.update(get_funcs_core[self.icase]) for opt, v in self.opts.items(): - self.links = nx.compose(self.links, links_opts[opt][v]) + self.graph = nx.compose(self.graph, graph_opts[opt][v]) self.get_funcs.update(get_funcs_opts[opt][v]) # Assign default values, if requested if use_default_values: @@ -135,14 +170,14 @@ def __init__(self, values=None, opts=None, use_default_values=True): for k, v in default_values.items(): if k not in values: values[k] = v - self.links.add_node(k) + self.graph.add_node(k) # Save arguments self.values = {} for k, v in values.items(): if k != "self" and v is not None: self.values[k] = v # state 1 means that the value was provided as an argument - nx.set_node_attributes(self.links, {k: 1}, name="state") + nx.set_node_attributes(self.graph, {k: 1}, name="state") def get(self, parameters, save_steps=True): """Calculate and return parameter(s) and (optionally) save them internally. @@ -165,11 +200,11 @@ def get(self, parameters, save_steps=True): parameters = [parameters] parameters = set(parameters) # get rid of duplicates # needs: which intermediate parameters we need to get the requested parameters - links_unknown = self.links.copy() - links_unknown.remove_nodes_from([v for v in self.values if v not in parameters]) + graph_unknown = self.graph.copy() + graph_unknown.remove_nodes_from([v for v in self.values if v not in parameters]) needs = parameters.copy() for p in parameters: - needs = needs | nx.ancestors(links_unknown, p) + needs = needs | nx.ancestors(graph_unknown, p) # The got counter increments each time we successfully get a value, either from # the arguments, already-calculated values, or by calculating it. # The loop stops once got reaches the number of parameters in `needs`, because @@ -189,39 +224,63 @@ def get(self, parameters, save_steps=True): got += 1 print("{} is available!".format(p)) else: - priors = self.links.pred[p] + priors = self.graph.pred[p] if len(priors) == 0 or all([r in self_values for r in priors]): print("Calculating {}".format(p)) self_values[p] = self.get_funcs[p]( *[ self_values[r] - for r in self.get_funcs[p].__code__.co_varnames + for r in self.get_funcs[p].__code__.co_varnames[ + : self.get_funcs[p].__code__.co_argcount + ] ] ) # state 2 means that the value was calculated internally if save_steps: - nx.set_node_attributes(self.links, {p: 2}, name="state") - for f in self.get_funcs[p].__code__.co_varnames: + nx.set_node_attributes(self.graph, {p: 2}, name="state") + for f in self.get_funcs[p].__code__.co_varnames[ + : self.get_funcs[p].__code__.co_argcount + ]: nx.set_edge_attributes( - self.links, {(f, p): 2}, name="state" + self.graph, {(f, p): 2}, name="state" ) results[p] = self_values[p] got += 1 print("Got", got, "of", len(set(needs))) + # Get rid of jax overhead on results + for k, v in results.items(): + try: + results[k] = v.item() + except: + pass if save_steps: + for k, v in self_values.items(): + try: + self_values[k] = v.item() + except: + pass self.values.update(self_values) return results - def plot_links(self, ax=None, show_tsp=True, show_missing=True, show_isolated=True): + def plot_graph( + self, + ax=None, + prog_graphviz="neato", + show_tsp=True, + show_unknown=True, + show_isolated=True, + ): """Draw a graph showing the relationships between the different parameters. Parameters ---------- ax : matplotlib axes, optional The axes, by default None, in which case new axes are generated. + prog_graphviz : str, optional + Name of Graphviz layout program, by default "neato". show_tsp : bool, optional Whether to show temperature, salinity and pressure nodes, by default False. - show_missing : bool, optional + show_unknown : bool, optional Whether to show nodes for parameters that have not (yet) been calculated, by default True. show_isolated : bool, optional @@ -235,26 +294,28 @@ def plot_links(self, ax=None, show_tsp=True, show_missing=True, show_isolated=Tr """ if ax is None: ax = plt.subplots(dpi=300, figsize=(8, 7))[1] - self_links = self.links.copy() - node_states = nx.get_node_attributes(self_links, "state", default=0) - edge_states = nx.get_edge_attributes(self_links, "state", default=0) + self_graph = self.graph.copy() + node_states = nx.get_node_attributes(self_graph, "state", default=0) + edge_states = nx.get_edge_attributes(self_graph, "state", default=0) if not show_tsp: - self_links.remove_nodes_from(["pressure", "salinity", "temperature"]) - if not show_missing: - self_links.remove_nodes_from([n for n, s in node_states.items() if s == 0]) + self_graph.remove_nodes_from(["pressure", "salinity", "temperature"]) + if not show_unknown: + self_graph.remove_nodes_from([n for n, s in node_states.items() if s == 0]) if not show_isolated: - self_links.remove_nodes_from( - [n for n, d in dict(self_links.degree).items() if d == 0] + self_graph.remove_nodes_from( + [n for n, d in dict(self_graph.degree).items() if d == 0] ) state_colours = {0: "xkcd:grey", 1: "xkcd:grass", 2: "xkcd:azure"} - node_colour = [state_colours[node_states[n]] for n in nx.nodes(self_links)] - edge_colour = [state_colours[edge_states[e]] for e in nx.edges(self_links)] - nx.draw_planar( - self_links, + node_colour = [state_colours[node_states[n]] for n in nx.nodes(self_graph)] + edge_colour = [state_colours[edge_states[e]] for e in nx.edges(self_graph)] + pos = nx.nx_agraph.graphviz_layout(self.graph, prog=prog_graphviz) + nx.draw_networkx( + self_graph, ax=ax, clip_on=False, with_labels=True, node_color=node_colour, edge_color=edge_colour, + pos=pos, ) return ax diff --git a/PyCO2SYS/equilibria/__init__.py b/PyCO2SYS/equilibria/__init__.py index 3dc3991..2a25fb9 100644 --- a/PyCO2SYS/equilibria/__init__.py +++ b/PyCO2SYS/equilibria/__init__.py @@ -6,8 +6,6 @@ from . import p1atm, pcx, pressured from .. import constants, convert, gas, solubility -__all__ = ["p1atm", "pcx", "pressured"] - def prepare(temperature, pressure, equilibria): """Initialise equilibria dict if needed and convert temperature/pressure units to diff --git a/PyCO2SYS/equilibria/p1atm.py b/PyCO2SYS/equilibria/p1atm.py index df0bacf..9032f5a 100644 --- a/PyCO2SYS/equilibria/p1atm.py +++ b/PyCO2SYS/equilibria/p1atm.py @@ -1,17 +1,39 @@ # PyCO2SYS: marine carbonate system calculations in Python. # Copyright (C) 2020--2024 Matthew P. Humphreys et al. (GNU GPLv3) -"""Estimate stoichiometric equilibrium constants at atmospheric pressure.""" - -from autograd import numpy as np -from .. import salts - - -def kCO2_W74(TempK, Sal): +""" +PyCO2SYS.equilibria.p1atm +========================= +Calculate stoichiometric equilibrium constants under standard atmospheric pressure. + +Functions +--------- +k_H2S_total_YM95 + Hydrogen sulfide dissociation constant on the total scale following YM95. +k_HF_free_DR79 + Hydrogen fluoride dissociation constant on the free scale following DR79. + Used when opt_k_HF = 1. +k_HF_free_PF87 + Hydrogen fluoride dissociation constant on the free scale following PF87. + Used when opt_k_HF = 2. +k_HSO4_free_D90a + Bisulfate dissociation constant in mol/kg-sw on the free scale following D90a. + Used when opt_k_HSO4 = 1. +k_HSO4_free_KRCB77 + Bisulfate dissociation constant in mol/kg-sw on the free scale following KRCB77. + Used when opt_k_HSO4 = 2. +k_HSO4_free_WM13 + Bisulfate dissociation constant in mol/kg-sw on the free scale following WM13, + with the corrections of WMW14. Used when opt_k_HSO4 = 3. +""" +from jax import numpy as np + + +def kCO2_W74(temperature, salinity): """Henry's constant for CO2 solubility in mol/kg-sw/atm following W74.""" # === CO2SYS.m comments: ======= # Weiss, R. F., Marine Chemistry 2:203-215, 1974. # This is in mol/kg-SW/atm. - TempK100 = TempK / 100 + TempK100 = (temperature + 273.15) / 100 lnK0 = ( -60.2409 + 93.4517 / TempK100 @@ -21,8 +43,24 @@ def kCO2_W74(TempK, Sal): return np.exp(lnK0) -def kHSO4_FREE_D90a(TempK, Sal): - """Bisulfate dissociation constant following D90a.""" +def k_HSO4_free_D90a(temperature, salinity, ionic_strength): + """Bisulfate dissociation constant in mol/kg-sw on the free scale following D90a. + Used when opt_k_HSO4 = 1. + + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + ionic_strength : float + Ionic strength in mol/kg-sw. + + Returns + ------- + float + HSO4 dissociation constant. + """ # === CO2SYS.m comments: ======= # Dickson, A. G., J. Chemical Thermodynamics, 22:113-127, 1990 # The goodness of fit is .021. @@ -30,22 +68,38 @@ def kHSO4_FREE_D90a(TempK, Sal): # TYPO on p. 121: the constant e9 should be e8. # Output KS is on the free pH scale in mol/kg-sw. # This is from eqs 22 and 23 on p. 123, and Table 4 on p 121: + TempK = temperature + 273.15 logTempK = np.log(TempK) - IonS = salts.ionic_strength_DOE94(Sal) - lnKSO4 = ( + lnk_HSO4 = ( -4276.1 / TempK + 141.328 - 23.093 * logTempK - + (-13856 / TempK + 324.57 - 47.986 * logTempK) * np.sqrt(IonS) - + (35474 / TempK - 771.54 + 114.723 * logTempK) * IonS - + (-2698 / TempK) * np.sqrt(IonS) * IonS - + (1776 / TempK) * IonS**2 + + (-13856 / TempK + 324.57 - 47.986 * logTempK) * np.sqrt(ionic_strength) + + (35474 / TempK - 771.54 + 114.723 * logTempK) * ionic_strength + + (-2698 / TempK) * np.sqrt(ionic_strength) * ionic_strength + + (1776 / TempK) * ionic_strength**2 ) - return np.exp(lnKSO4) * (1 - 0.001005 * Sal) - - -def kHSO4_FREE_KRCB77(TempK, Sal): - """Bisulfate dissociation constant following KRCB77.""" + return np.exp(lnk_HSO4) * (1 - 0.001005 * salinity) + + +def k_HSO4_free_KRCB77(temperature, salinity, ionic_strength): + """Bisulfate dissociation constant in mol/kg-sw on the free scale following KRCB77. + Used when opt_k_HSO4 = 2. + + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + ionic_strength : float + Ionic strength in mol/kg-sw. + + Returns + ------- + float + HSO4 dissociation constant. + """ # === CO2SYS.m comments: ======= # Khoo, et al, Analytical Chemistry, 49(1):29-34, 1977 # KS was found by titrations with a hydrogen electrode @@ -60,13 +114,30 @@ def kHSO4_FREE_KRCB77(TempK, Sal): # The rms error is .0021 in pKS, or about .5% in KS. # This is equation 20 on p. 33: # Output KS is on the free pH scale in mol/kg-sw. - IonS = salts.ionic_strength_DOE94(Sal) - pKSO4 = 647.59 / TempK - 6.3451 + 0.019085 * TempK - 0.5208 * np.sqrt(IonS) - return 10.0**-pKSO4 * (1 - 0.001005 * Sal) + TempK = temperature + 273.15 + pk_HSO4 = ( + 647.59 / TempK - 6.3451 + 0.019085 * TempK - 0.5208 * np.sqrt(ionic_strength) + ) + return 10.0**-pk_HSO4 * (1 - 0.001005 * salinity) + +def k_HSO4_free_WM13(temperature, salinity): + """Bisulfate dissociation constant in mol/kg-sw on the free scale following WM13, + with the corrections of WMW14. Used when opt_k_HSO4 = 3. -def kHSO4_FREE_WM13(TempK, Sal): - """Bisulfate dissociation constant following WM13/WMW14.""" + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + + Returns + ------- + float + HSO4 dissociation constant. + """ + TempK = temperature + 273.15 logKS0 = ( 562.69486 - 102.5154 * np.log(TempK) @@ -81,27 +152,55 @@ def kHSO4_FREE_WM13(TempK, Sal): + 0.0267059 * TempK * np.log(TempK) - 0.000042128 * TempK**2 ) - * Sal**0.5 - + (0.2542181 - 0.00509534 * TempK + 0.00071589 * TempK * np.log(TempK)) * Sal - + (-0.00291179 + 0.0000209968 * TempK) * Sal**1.5 - + -0.0000403724 * Sal**2 + * salinity**0.5 + + (0.2542181 - 0.00509534 * TempK + 0.00071589 * TempK * np.log(TempK)) + * salinity + + (-0.00291179 + 0.0000209968 * TempK) * salinity**1.5 + + -0.0000403724 * salinity**2 ) - kSO4 = (1 - 0.001005 * Sal) * 10.0 ** (logKSK0 + logKS0) - return kSO4 + k_HSO4 = (1 - 0.001005 * salinity) * 10.0 ** (logKSK0 + logKS0) + return k_HSO4 + + +def k_HF_free_DR79(temperature, salinity, ionic_strength): + """Hydrogen fluoride dissociation constant on the free scale following DR79. + Used when opt_k_HF = 1. + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. -def kHF_FREE_DR79(TempK, Sal): - """Hydrogen fluoride dissociation constant following DR79.""" + Returns + ------- + float + HF dissociation constant. + """ # === CO2SYS.m comments: ======= # Dickson, A. G. and Riley, J. P., Marine Chemistry 7:89-99, 1979: # this is on the free pH scale in mol/kg-sw - IonS = salts.ionic_strength_DOE94(Sal) - lnKF = 1590.2 / TempK - 12.641 + 1.525 * IonS**0.5 - return np.exp(lnKF) * (1 - 0.001005 * Sal) + lnKF = 1590.2 / (temperature + 273.15) - 12.641 + 1.525 * ionic_strength**0.5 + return np.exp(lnKF) * (1 - 0.001005 * salinity) + +def k_HF_free_PF87(temperature, salinity): + """Hydrogen fluoride dissociation constant on the free scale following PF87. + Used when opt_k_HF = 2. -def kHF_FREE_PF87(TempK, Sal): - """Hydrogen fluoride dissociation constant following PF87.""" + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + + Returns + ------- + float + HF dissociation constant. + """ # Note that this is not currently used or an option in CO2SYS, # despite the equations below appearing in CO2SYS.m (commented out). # === CO2SYS.m comments: ======= @@ -110,7 +209,7 @@ def kHF_FREE_PF87(TempK, Sal): # Nonetheless, P&F87 might actually be better than the fit of D&R79 above, # which is based on only three Salinities: [0 26.7 34.6] # Output is on the free pH scale in mol/kg-SW. - lnKF = 874 / TempK - 9.68 + 0.111 * Sal**0.5 + lnKF = 874 / (temperature + 273.15) - 9.68 + 0.111 * salinity**0.5 return np.exp(lnKF) @@ -745,8 +844,21 @@ def kH2CO3_TOT_PLR18(TempK, Sal): return K1, K2 -def kH2S_TOT_YM95(TempK, Sal): - """Hydrogen sulfide dissociation constant following YM95.""" +def k_H2S_total_YM95(temperature, salinity): + """Hydrogen sulfide dissociation constant on the total scale following YM95. + + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + + Returns + ------- + float + H2S dissociation constant. + """ # === CO2SYS_v1_21.m comments: ======= # H2S Millero et. al.( 1988) Limnol. Oceanogr. 33,269-274. # Yao and Millero, Aquatic Geochemistry 1:53-88, 1995. Total Scale. @@ -754,12 +866,13 @@ def kH2S_TOT_YM95(TempK, Sal): # they agree with Millero 1988 which are on Total Scale. # Also, calculations agree at high H2S with AquaEnv when assuming it is on # Total Scale. + TempK = temperature + 273.15 lnkH2S = ( 225.838 - 13275.3 / TempK - 34.6435 * np.log(TempK) - + 0.3449 * np.sqrt(Sal) - - 0.0274 * Sal + + 0.3449 * np.sqrt(salinity) + - 0.0274 * salinity ) return np.exp(lnkH2S) @@ -790,14 +903,9 @@ def kNH3_TOT_CW95(TempK, Sal): -13.6416 + 1.176949 * TempK**0.5 - 0.02860785 * TempK + 545.4834 / TempK ) * Sal**0.5 PKNH3expCW += ( - -0.1462507 - + 0.0090226468 * TempK**0.5 - - 0.0001471361 * TempK - + 10.5425 / TempK + -0.1462507 + 0.0090226468 * TempK**0.5 - 0.0001471361 * TempK + 10.5425 / TempK ) * Sal**1.5 - PKNH3expCW += ( - 0.004669309 - 0.0001691742 * TempK**0.5 - 0.5677934 / TempK - ) * Sal**2 + PKNH3expCW += (0.004669309 - 0.0001691742 * TempK**0.5 - 0.5677934 / TempK) * Sal**2 PKNH3expCW += (-2.354039e-05 + 0.009698623 / TempK) * Sal**2.5 KNH3 = 10.0**-PKNH3expCW # this is on the total pH scale in mol/kg-H2O KNH3 = KNH3 * (1 - 0.001005 * Sal) # convert to mol/kg-SW diff --git a/newclass_v2.py b/newclass_v2.py index 5b92b35..796fa70 100644 --- a/newclass_v2.py +++ b/newclass_v2.py @@ -1,15 +1,26 @@ -from PyCO2SYS import CO2System +from PyCO2SYS import CO2System, system import networkx as nx sys = CO2System( # dict(salinity=32), - # opts=dict(opt_total_borate=2), + opts=dict(opt_k_HF=2), # use_default_values=False, ) -sys.get(("total_borate", "total_sulfate")) -sys.plot_links( - # show_missing=False, +sys.get( + ( + # "total_sulfate", + # "k_HSO4_free_1atm", + # "ionic_strength", + # "k_HF_free_1atm", + # "total_to_sws_1atm", + # "nbs_to_sws", + "k_H2S_sws_1atm", + ) +) +sys.plot_graph( + # show_unknown=False, # show_isolated=False, + # prog_graphviz='neato', ) print(sys.values)