Corrected calcite solubility parameterisation

[jagoosw]

As caught in #367

[jagoosw]

Conversation from #367

jagoosw:
Do you know where this is in the paper? I can't see any log10s?

jagoosw:

jagoosw:
I think there is a problem with the current calcite saturation, I think it is off by a factor of 1000, but I thought currently that it might be a unit problem with the Ca+ ions?

johnryantaylor:
It is in Mucci 1983 page 792. It is confusing because it just says "log", but below equation 10, the text says that log refers to log_10. It is worth double checking with another source, though. Presumably CO2SYS has this formula in it?

johnryantaylor:
The image that you posted does have ln(). Is this from Millero?

jagoosw:
Yeah this is from Millero, I guess that this is a mistake in Millero because pyco2sys has 10^

jagoosw:
I'll make a separate PR to correct this

johnryantaylor:
This is what CO2SYS has:
% CalciteSolubility:
% ' Mucci, Alphonso, Amer. J. of Science 283:781-799, 1983.
logKCa = -171.9065 - 0.077993.*TempK + 2839.319./TempK;
logKCa = logKCa + 71.595.*logTempK./log(10);
logKCa = logKCa + (-0.77712 + 0.0028426.*TempK + 178.34./TempK).*sqrSal;
logKCa = logKCa - 0.07711.*Sal + 0.0041249.*sqrSal.*Sal;
% ' sd fit = .01 (for Sal part, not part independent of Sal)
KCa = 10.^(logKCa);% ' this is in (mol/kg-SW)^2

johnryantaylor:
That is the Matlab implementation of CO2SYS.

johnryantaylor:
While I'm looking at it, the pressure correction term for aragonite in CO2SYS looks very different:
% PressureCorrectionForAragonite:
% ' Millero, Geochemica et Cosmochemica Acta 43:1651-1661, 1979,
% ' same as Millero, GCA 1995 except for typos (-.5304, -.3692,
% ' and 10^3 for Kappa factor)
deltaVKAr = deltaVKCa + 2.8;
KappaKAr = KappaKCa;
lnKArfac = (-deltaVKAr + 0.5.*KappaKAr.*Pbar).*Pbar./RT;
KAr = KAr.*exp(lnKArfac);

The pressure correction for calcite is the same as what we have. Is it possible that this needs to be updated too?

johnryantaylor:
Sorry, scratch that last comment - I hadn't spotted that CO2SYS is just modifying the calcite values. This looks fine

jagoosw:
I found how this slipped through, the book I used for the check values also just lists "log", and I realise now they use "ln" for the natural log

[johnryantaylor]

lgtm, thanks!

[johnryantaylor]

lgtm, thanks!

[jagoosw]

I was just looking at the MARBL source code and I think they're making the same mistake here too:
https://github.com/marbl-ecosys/MARBL/blob/6b2a5e41dc5b0faaf9b600c96432e2487515e878/src/marbl_co2calc_mod.F90#L1127-L1136