Currently we follow CO2SYS-MATLAB in using a finite difference scheme. Calculate directly instead!

The aim is to make all of the internal functions available to use without needing to run the full CO2SYS function. That means the internal functions need adding to the docs.

Modules to document

• assemble
• buffers
• concentrations
• constants
• convert
• equilibria
• meta
• original
• solubility
• solve

A paper by Honkanen et al. (2021) used a total boron parameterization by Kuliński et al. (2018). Is this implemented in PyCO2SYS? If not, is there interest in implementing it? This option may make the program more useful to researchers in brackish waters and/or the Baltic Sea.

In particular, the new NH3 and H2S inputs that were added.

At the moment, this calculation requires the k_constants["FugFac"] input. This makes it a bit bulky if you'd quickly like to convert pCO2 --> fCO2. All the components are there (i.e., gas.fugacityfactor), it would just be a matter of reorganising a little bit.

One possible way of doing this would be to call the current function convert._fCO2_to_pCO2 that gets called by the existing scripts in PyCO2SYS. The old name could then be used for a standalone function that uses gas.fugacityfactor the calculate the factor. This approach would require some refactoring. I had a quick look at the code and it could also be xarray friendly out of the box 👍

Alternatively, we can add it to the api.

For error propagation and for buffer factor calculations

See paper "Generalised expressions for the response of pH to changes in ocean chemistry" (Hagens & Middelburg, 2016).

The code of RGasConstant suggests that the function argument WhichR can be supplied as an array and that the returned variable RGas will be of the same shape as WhichR.

However, if WhichR is an array the isequal test (e.g. WhichR == 1) will always return false, leaving the default value in RGas, which is nan.
Examples

import PyCO2SYS as pyco2
import numpy as np

# As expected:
pyco2.constants.RGasConstant(1) # returns array(83.1451) as expected
pyco2.constants.RGasConstant(2) # returns array(83.14472) as expected
pyco2.constants.RGasConstant(np.array([1,1])) # returns array([83.1451, 83.1451]) as expected

# Not as expected:
pyco2.constants.RGasConstant([1]) # returns array([nan]), expected array([83.1451])
pyco2.constants.RGasConstant([1,2]) # returns array([nan, nan]), expected array([83.1451, 83.14472])
pyco2.constants.RGasConstant(None) # returns array(nan), expected array([83.1451, 83.14472])

# Expected a default value or an error:
pyco2.constants.RGasConstant(5) #  WhichR is out of range, this should error or return the default
pyco2.constants.RGasConstant() # The input argument should really be optional and the function should return the default value if none is supplied
pyco2.constants.RGasConstant(np.array([1.3, 2.345])) # The inputs are no valid integers

So input checks are necessary:

1. All WhichR inputs that are not scalar need to be numpy arrays of integers in the range from 1 to 3
2. For invalid inputs it might be sensible to return a default value

Please correct me, if this behaviour is intentional or I misunderstood anything.

Thanks for maintaining this awesome tool!

Pressure correction functions for equilibrium constants are currently in assemble.equilibria exactly as originally written in the MATLAB. To do:

• Extract these correction equations into their own subfunctions.
• Add the related citations to the docs references page.
• Add the pressure corrections directly into the functions in the equilibria module

Currently in the __init__.py

• Extract _CaSolubility and from the __init__.py into new module solubility.
• Break up _CaSolubility into sensible subfunctions.
• Add relevant citations to the docs references page.
• Make it all Autograd-able

Now they just give the conversion factors - better to make them do the actual conversions instead
All in module convert

Return logical arrays indicating whether calculated results fall in valid input ranges (e.g. of temperature, salinity and pressure) given the sets of constants used.

One way to do this:

• Add validity output to equilibrium constant functions
• Accumulate validity outputs through the different calculations
• Report validity results in main CO2SYS output

Another way would be to build a separate function that just checks validity, rather than adding extra outputs to the equilibrium constant functions and so on. This is less robust but less breaking?

I am attempting to do computations of the carbonate system in model output. I use xarray to load and analyze these data with lazy computations, meaning that data are only loaded into memory as a final step, and then only the data required for the computation. This makes parallelization of the computation easy when working with large datasets.

This functionality doesn't seem to play nicely with PyCO2SYS (or vice versa!), which instead loads all the data into memory during the application of the sys command. I'm not sure how to put together a good minimum working example to show the issue here, because it relies on loading a dataset. I can try to put something together if the issue is not clear.

I'm wondering if there is a plan or interest to integrate the functionality of PyCO2SYS with that of xarray? That is, to allow the calculations inherent in PyCO2SYS to be performed lazily. I would be happy to help implement this, although I am no expert. On the other hand, I would be happy to hear if there is a workaround that I am missing.

Tagging @ognancy4life and @jbusecke as folk who might be interested in this problem.

P.S. Thanks so much for implementing this package - it's a great resource!

How does pyco2sys parametrize the pressure correction, and is there a way to select what correction to use? I am comparing the output against (1995), but the values are quite different. Reference values below are from Zeebe & Gladrow 2001 CO2 in Seawater: Equilibrium, Kinetics, Isotopes (page 268), who use Roy 1993 for pk1 & pk2. The values for 0 bar match, but those for 300 bar are quite different.

import PyCO2SYS as pyco2
from math import log, log10, exp

# defaults for S & T are 35 & 25
r = pyco2.sys(
    pressure=1,
    par1_type=1,  # "1" =  "alkalinity"
    par1=2399, # umol/kg 
    par2_type=2,  # "1" = dic
    par2=2291, # umol/kg 
    opt_k_carbonic=1, # 1: Roy et al 
    opt_pH_scale=1, # 1:total, 3:free
    opt_buffers_mode=2, # 
)

print(f"pk1 = {-log10(r['k_carbonic_1']):.4f} vs 5.8563 @ 0 bar")
print(f"pk2 = {-log10(r['k_carbonic_2']):.4f} vs 8.9249 @ 0 bar")

and for 300 bar

r = pyco2.sys(
    pressure=301,
    par1_type=1,  # "1" =  "alkalinity"
    par1=2399, # umol/kg 
    par2_type=2,  # "1" = dic
    par2=2291, # umol/kg 
    opt_k_carbonic=1, # 1: Roy et al 
    opt_pH_scale=1, # 1:total, 3:free
    opt_buffers_mode=2, # 
)

print(f"pk1 = {-log10(r['k_carbonic_1']):.4f} vs 5.7392 @ 300 bar")
print(f"pk2 = {-log10(r['k_carbonic_2']):.4f} vs 8.4746 @ 300 bar")

• __init__
• assemble
• buffers
• concentrations
• convert
• equilibria
• meta
• original
• solubility
• solve

Into convert module?

I am building a package that depends on some pyCO2sys routines. I would like to release this package via conda. However, that will only work if all dependencies (including pyCO2sys) can be satisfied by the conda package manager.

I see that adding pyCO2sys to conda has been attempted before: conda-forge/staged-recipes#21796

But apparently never succeeded. A quick check suggests that all of pyCO2sys's requirements exist on conda, so this should be straightforward unless there is OS-dependent code somewhere.

Are there plans to try this again? Also, I am happy to help with this.

Make a 'black box' function that allows a user to input a CSV/spreadsheet file and receive all the CO2SYS outputs in the same format.

Seacarb (R)

https://cran.r-project.org/web/packages/seacarb/index.html

mocsy (Fortran)

https://github.com/jamesorr/mocsy

CO2SYSv2 (MATLAB)

https://github.com/jamesorr/CO2SYS-MATLAB

CO2SYS (DOS)

Add the factors derived by Egleston et al. (2010) as additional outputs in CO2dict.

Note the typo corrections given by Richier et al. (2018) - see their supplementary info files.

Can modify from MATLAB version: https://github.com/mvdh7/oceancarb-constants/blob/master/ESM10_CO2SYS.m

Add isocapnic quotient of Humphreys et al. (2018) to CO2dict output.

Can use MATLAB scripts as a guide: https://github.com/mvdh7/biogeochem-phi.

assemble has become too random and unwieldy. Break up into more sensible chunks (e.g. engine)