mvdh7 / PyCO2SYS

Marine carbonate system calculations in Python
https://PyCO2SYS.readthedocs.io
GNU General Public License v3.0
50 stars 9 forks source link

Pressure Correction #68

Closed uliw closed 8 months ago

uliw commented 8 months ago

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")
d-sandborn commented 8 months ago

I obtain the following result when running the code above (using PyCO2SYS 1.8.2):

pk1 = 5.8563 vs 5.8563 @ 0 bar pk2 = 8.9249 vs 8.9249 @ 0 bar pk1 = 5.8445 vs 5.7392 @ 300 bar pk2 = 8.9162 vs 8.4746 @ 300 bar

The pressure term in pyco2.sys() should be given in dbar, as given in the docs. Additionally, the check values for 300 bar given in the code above seem to disagree with those in the text. In my copy of CO2 in Seawater, table A.11.2 on page 268 lists the following check values:

image

Revising the code with pressure = 3000 and substituting the text values yields:

pk1 = 5.7389 vs 5.7397 @ 300 bar pk2 = 8.8401 vs 8.8409 @ 300 bar

These values are in considerably better agreement, differing by 0.14 ‰ and 0.090 ‰.

uliw commented 8 months ago

as usual, RTFM .... Thank you so much, and my apologies!

mvdh7 commented 8 months ago

Thanks @d-sandborn !