Carbonate Chemistry Algorithm

[j-ostiguy]

Hello,

I've done a comparison between the carbonate chemistry solver implemented here, the algorithm presented in this paper, and the CO2SYS algorithm (CO2System.jl). For a range of DIC/Alkalinity or temperature, the calculated pCO2 can differ by up to ~5 to 15 uatm which is shown in the figure below.

I've also run Oceananigans simulations with both the OceanBioME carbonate solver and the Follows (2006) algorithm and found that the choice of carbonate solver results in DIC fields that are quite different. I'm wondering if you have any thoughts on this or if this kind of comparison has already been done. I'm also curious to know if the algorithm chosen here has some advantage over the others. 

Thank you!

[johnryantaylor]

Thanks for the comment. I don't immediately know the answer. In order to look into this further, can you let us know which version of OceanBioME you are using and what value of salinity you are using for this test?

[jagoosw]

Hi @j-ostiguy,

This is interesting thank you. Firstly for the simulations with DIC fields, I found a bug in the air/sea flux code that I'm fixing in #303 at the moment is now fixed in main.

As for the carbon chemistry, it is meant to follow "Dickson, A.G., Sabine, C.L. and Christian, J.R. (2007), Guide to Best Practices for Ocean CO 2 Measurements. PICES Special Publication 3, 191 pp." which looks similar to the Follows paper you linked. Could you share the exact setup you used to calculate them?

[jagoosw]

Could this be because CarbonChemistry returns fCO2 not pCO2, we convert it here for the air/sea exchange:

OceanBioME.jl/src/Models/GasExchange/carbon_dioxide_concentration.jl

Lines 31 to 59 in 3c7465a

	@inline function surface_value(cc::CarbonDioxideConcentration, i, j, grid, clock, model_fields)
	DIC = @inbounds model_fields.DIC[i, j, grid.Nz]
	Alk = @inbounds model_fields.Alk[i, j, grid.Nz]
	T = @inbounds model_fields.T[i, j, grid.Nz]
	S = @inbounds model_fields.S[i, j, grid.Nz]
	FT = eltype(grid)
	silicate, phosphate = silicate_and_phosphate(cc.silicate_and_phosphate_names, model_fields)
	fCO₂ = cc.carbon_chemistry(; DIC, Alk, T, S, silicate, phosphate)
	P = surface_value(cc.air_pressure, i, j, grid, clock) * convert(FT, ATM)
	Tk = T + convert(FT, 273.15)
	B = cc.first_virial_coefficient(Tk)
	δ = cc.cross_virial_coefficient(Tk)
	xCO₂ = fCO₂ / P
	# Experimentally this converged xCO₂ to machine precision
	for n = 1:3
	φ = P * exp((B + convert(FT, 2) * (one(FT) - xCO₂)^convert(FT, 2) * δ) * P / (convert(FT, GAS_CONSTANT) * Tk))
	xCO₂ = fCO₂ / φ
	end
	return xCO₂ * P # mol/mol to ppm
	end

Perhaps when I implemented this I miss understood something about the algorithm.

[jagoosw]

I checked with PyCO2SYS and the value I get from cc(; DIC = 2000.0, Alk = 2100.0, T = 15.0, S = 35.0) is very close to pyco2.sys(par1=1949.317738791423, par2=2046.7836257309941, par1_type=2, par2_type=1, temperature = 15, salinity = 35)["fCO2"] (724.5881474791435 vs 724.454609155267)

[johnryantaylor]

Are par1 and par2 the values of DIC and Alk? Where do the numbers in the co2sys call come from (e.g. par1=1949.3177...)?

[jagoosw]

Yes sorry should have clarified, in co2sys you specify what par1 and par2 are by the par1/2_type which specifies that they're DIC and Alkalinity, and co2sys uses $\mu$ mol/kg so the values are *10^3/1026 (so that might be where the small difference comes from)

[jagoosw]

When I use oceanbiome to convert first I get 724.5290381997396 from co2sys

[johnryantaylor]

Ok, that makes sense, thanks. How different are the values of fCO2 and pCO2 from CO2SYS? I would have thought that they would be very close to each other. Usually fCO2 would be preferred, but we could make it an option in cc to output pCO2 instead (similar to how we optionally output pH).

[jagoosw]

For the same inputs pCO2 is 727.1548437614236 in co2sys.

Okay I'll make a PR to add them all to be options, perhaps pCO2 should be the default since thats whats used in the air/sea flux calculation?

[jagoosw]

Perhaps I should add a function to output pCO2, or change the default output of cc(; DIC...) to be pCO2

[johnryantaylor]

I think fCO2 should be the default. This includes corrections for non-ideal behavior of the gas. pCO2 could be useful for comparing with other models (as in this case).

[j-ostiguy]

The difference is because of the units - I was using the same input Alk/DIC for all three algorithms but didn't realize that carbon_chemistry takes in mmol/m^3. It could be helpful to clarify this more explicitly in the documentation (independently of LOBSTER, which I just found lists the units). Thank you for looking into this so thoroughly!

[jagoosw]

No problem! If you remake the figure it would be interesting to see.

I'm going to add a note to the docstrings about the units

[j-ostiguy]

Here is the figure, for a wider range of T/Alk/DIC (and S = 35) with results for fCO2. For this range, the fCO2 differs by less than 0.5% between the different algorithms.

[jagoosw]

Thank you thats really helpful to see!