oscarbranson / cbsyst Goto Github PK

Hello,
Firstly, thank you very much for this great work! I have quickly compared the outputs to pyCO2sys and not only is a near 1:1, it is also about 100 times faster! Amazing!

However, I ran into a small issue when trying to execute the given example in the README.
cbsyst.Csys and cbsyst.CBsys appear to run just fine, but when I try to run the following from you example:

  import cbsyst as cb
  import numpy as np

  pH = np.linspace(7,11,100)
  Bsw = cb.Bsys(pHtot=pH, BT=433., dBT=39.5)

It gives me the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.8/dist-packages/cbsyst/cbsyst.py", line 324, in Bsys
    ps.update(ABsys(pdict=ps))
  File "/usr/local/lib/python3.8/dist-packages/cbsyst/cbsyst.py", line 435, in ABsys
    ps.alphaB = alphaB_calc(ps.T)
AttributeError: 'Bunch' object has no attribute 'T'

So I went in the ABsys function and printed the Bunch object just before this line and found that there was indeed no attribute "T", but instead there was "T_in". So if I change the line from what it was to:

ps.alphaB = alphaB_calc(ps.T_in)

It works, but then the same kind of error arises when hitting this loop:

for k in ['ABO3', 'ABO4', 'ABT', 'Ca',
              'H', 'Mg', 'S', 'T', 'alphaB',
              'dBO3', 'dBO4', 'dBT', 'pHtot']:
        if not isinstance(ps[k], np.ndarray):
            ps[k] = np.array(ps[k], ndmin=1)

As in the current Bunch, there is no T and S, but T_in and S_in. So if you change them to give:

for k in ['ABO3', 'ABO4', 'ABT', 'Ca',
              'H', 'Mg', 'S_in', 'T_in', 'alphaB',
              'dBO3', 'dBO4', 'dBT', 'pHtot']:
        if not isinstance(ps[k], np.ndarray):
            ps[k] = np.array(ps[k], ndmin=1)

Then it works flawlessly.
I do not know if it wise to change it like this, because I do not know how you treat the "in" and "out" of the physical properties of water samples. That is why I posted it here.

By the way, I run python 3.8.5.

Thanks,
Douglas

Relatively simple, needs to happen.

Trivial, but not done yet.

Thought: when is the most efficient place to do this?

Implement Alternate [Mg] and [Ca] corrections

• Ben-Yaakov and Goldharber (1973)
• Tyrrell and Zeebe (2004)
• PHREEQC Pitzer model?

Would other constants be useful?

At present, there are no options for different constants in cbsyst. Do people want/need choices?

Request constants in the comments below. If they get traction (i.e. lots of +1's), we'll implement them.

At present:

Constants as a function of Temperature and Salinity:

• K₀, K₁, K₂, K_B, K_W, K_SO4 from Chapter 5 of Dickson, Sabine and Christian (2007), corrected for [Mg] and [Ca] after Hain et al (2015)
• K_sp(arag) and K_sp(cal) from Zeebe and Wolf-Gladrow (2001)
• K_P1, K_P2, K_P3, K_Si from from Chapter 5 of Dickson, Sabine and Christian (2007) not corrected for [Mg] and [Ca].
• K_F from Dickson and Riley (1979)

Pressure Corrections
All pressure correction parameters from Millero (2007), with modifications:

• TYPO: Third parameter for K_B was positive, should have been negative
• TYPO: All 'b' parameters were missing a factor of 10^-3
• TYPO: K_W parameters in original were for fresh water - replaced with seawater values from Millero (1979)

Conserved Seawater Compositions:
As a function of salinity:

• TS of Morris and Riley (1966)
• TB of Lee et al (2010) Uppstrom, (1974) IF no value provided. You can specify BT at S=35, and it scales with salinity.
• TF of Riley (1965)

Implement it... using uncertainties package?

This is low priority but would be nice.
Current incompatibility stems from repeated use of 'Extended Tuple Unpacking' (using *) in numerous function inputs throughout.
Can be solved with a generator wrapper (as described here).

• Replace all instances of Extended Tuple Unpacking with generators.
• Make import structure python 2 & 3 compatible.

Note: To run tests in python 2.x, you need to remove the python 2 incompatability warning in setup.py

cbsyst errors out when some inputs are set to None due to a bug in the current version of Kgen.

Steps to reproduce:
Install cbsyst

python3.10 -m venv .environment
source .environment/bin/activate
pip install cbsyst

Run a test

python
from cbsyst import Csys

result = Csys(pHtot=8.0,DIC=2000/1e6,unit="mol") // Produces error in kgen

Currently, input and output conditions are the same.
It would be good to be able to separate input/output T, S, P conditions.

This should be relatively simple and inexpensive, and involve a second call to MyAMI_get_Ks and re-calculating the final speciations.
Suggest wrapping this in a separate function which takes Ks, and outputs all conditions from a single parameter combination (e.g. DIC and H for Csys).

A lot of the solvers in carbon_fns.py would be sped up by vectorising.
At the moment, most serially apply zero_finder functions to solve for H.

The following functions in cbsyst/carbon_fns.py need to be vecotrised:

• CO2SYS-style zero finder for HCO3_TA (case 11)
• CO2SYS-style zero finder for CO3_TA (case 13)
• Vectorise CO2_HCO3 (case 2)
• Vectorise CO2_CO3 (case 3)
• Vectorise CO2_DIC (case 5)
• Vectorise HCO3_CO3 (case 10)
• Vectorise HCO3_DIC (case 12)
• Vectorise CO3_DIC (case 14)

To fix deviation from GLODAPv2 at depth.

Clarification

cbsyst should not be used to calculate C system at different [Mg] and [Ca] conditions until this issue is resolved. Modern ocean calculations are unaffected by this problem.

Issue

Output of MyAMI_K_calc is incorrect at non-ambient [Ca] and [Mg] (c.v. sensitivies listed in Hain et al, 2015).

This is an underlying issue in MyAMI. Hain's original code produces different K values when using PyMyAMI.m compared to MyAMI_V1.py.

Main difference in code is that PyMyAMI.m calculates correction factor for experimental condition directly, while MyAMI_V1.py calculates a grid of corrected Ks at different (T,S) conditions, fits a K_conditional function to it, and then returns parameters to be used for K calculation across a temperature range. Uncertain why these produce different results.

Resolution

• Contact Hain
• Either: Implement direct experimental correction method (as perPyMyAMI.m), or fix grid calculation method in MyAMI_K_calc function.

Dear developers,

thank you for your help and your work!
I have following issue, when I am running the script(s).

partially initialized module 'cbsyst' has no attribute 'Bsys' (most likely due to a circular import)

Do you have any clou, how to solve that problem?

best regards, Matthias

Currently using CO2SYS method. Consider others:

• As in LOSCAR (Follows et al, 2006).
• SolveSAPHE

After Orr et al (2015).

This is shelved for the moment - revisit later if it becomes an issue

At the moment, the Xsys functions can only take single combinations of input parameters.
It might be useful to allow for the calculation of multiple combinations of input combinations at once (a la CO2SYS.m).
However, algorithms currently identify parameters to be calculated by whether they are 'None' or not, so this would require a major re-tooling.

I reckon the effort involved here would be much greater than the inconvenience for the user having to run the function twice for different sets of input parameters.
How often do people actually use the ability of CO2SYS.m to take more than one combination of input parameters?

In the list of input from the Bsys or CBsys functions, the alphaB parametre can be manually typed. However, this parametre is never really used. The ps Bunch object is automatically updated using the alphaB_calc function (dependant on temperature, based on the Hoenisch book) and simply ignores the alphaB parametre (i.e. cbsyst.py, line 692 and 1037).
Would it be possible to make a simple statement checking if the alphaB is given, and if not calculate it?

CO2 can become negative under high pH, high TA:

cb.Csys(pHtot=10.5, TA=2300).CO2

# array([-9.01237887e-06])

Csys hangs if pH is too high:

cb.Csys(pHtot=11, TA=2300)

# does not run

Suspect that slight negative trend in DIC predicted from TA and pH could be because of incorrect K temperature sensitivity.

Check this!

Current vectorised fsolve calls in carbon_fns produce incorrect values when multiple parameters are passed. Identified from Luecker et al (2000) test data. Unclear why.

Solution:

• Cast all inputs into ndarray and iterate over rows of inputs. This gives correct results, but is slower.
• When to do this?
  • At [N]sys level, after Ks calculation? Brings problems with identification of None params?
  • But these problems need to be fixed anyway, because at the moment doesn't allow for simultaneous calculation of different input calculations.

• Total
• SWS
• Free
• NBS

import cbsyst as cb

cb.Csys(pHtot=5., pCO2=300, T_in=26, S_in=33)

Yields a TA of -8.92526313.