Refactors gas exchange carbon chemistry (+ adds some complexity)

[jagoosw]

This PR refactors the carbon chemistry calculations improving the clarity and supersedes #179.

It also moves both the carbon chemistry and gas exchange functions to Models since they are models, and separates them.

This PR also implements the full carbon chemistry standards of the best practice (https://www.ncei.noaa.gov/access/ocean-carbon-acidification-data-system/oceans/Handbook_2007/Guide_all_in_one.pdf) so it can depend on Si and PO4, and approximated fluoride and sulfate. It also introduces proper (TEOS10 and polynomial approximations) density calculations.

Also, all of the equilibrium and pressure correction constants are also tested now.

I also introduce a function calcite_saturation that can compute calcite_concentration and saturation for use in other models.

This should also be faster since it removed a branching operation.

Finally, I changed the structure of GasExchange to make it more flexible (and correct, there was a bug in the computation of x(CO2) from fCO2).

Also, closes https://github.com/OceanBioME/OceanBioME.jl/issues/197, #200

[jagoosw]

@bethcandish you could do some testing with the "best practice" parameters with this branch and see if it makes any difference, and if so change them here.

[jagoosw]

@johnryantaylor do you think it would make sense to move carbon chemistry to be a property of Biogeochemistry?

Previously we had only used it in the CO2 flux but PISCES needs the calcite saturation for various things, and it makes sense to me that we separate it so it can be used for any model.

It wouldn't cause many changes and will just be there to be called, e.g. by the boundary condition.

[jagoosw]

Also, I will update the docs page at some point but the code now includes the pressure correction for the rate constants. I have added some tests which ensure we have the constants right (when using the "best practice" parameters) like how Kgen does, so I think the only difference between the functionality here and cbsyst on the carbon front is the Mg and Ca bits.

[jagoosw]

Final comment for now, I'm pretty sure that this PR isn't breaking since OCMPI_default wasn't part of the public API, but it's not clear to me if it counts as breaking because lots of things have moved around. I don't really want to make a breaking release because it was a lot of effort to bump the compact of all my packages that rely on it like the giant kelp stuff.

[jagoosw]

I downloaded all of the latest GLODAP data (~200000 complete records) and get pH error to be 0.00022 ± 0.01223 with the "best practice" parameters (vs 0.006 ± 0.012 with the parameters from the NEMO source code).

Also, in the latest GLODAP, there are 9 cruises which report 10059 directly measured pCO2 values. With the "best practice" parameters I get an error of -3 ± 28.8 (vs 11 ± 28 for the NEMO values).

It also looks to me like the error variance is larger at lower temperatures and greater error at lower pH.

Given these results, I think we should change the defaults to the "best practice" parameters.

(for both pH and pCO2 I excluded the "calculated" values with flag 0 since it seems silly to compare to calculated values which presumably use a very similar algorithm and cause a huge spike a histogram of the error)

@johnryantaylor @oscarbranson

[jagoosw]

Oh dear GibbsSeaWater/ccall doesn't work on GPU

Edit: alternative using SeawaterPolynomials sorted now

[jagoosw]

Docs preview is done now and includes much better details of whats going on: https://oceanbiome.github.io/OceanBioME.jl/previews/PR186/model_components/carbon-chemistry/

(+/- typos which I think I've fixed now)

[jagoosw]

@vtamsitt the refactoring of the gas exchange model here should make adding different gas transfer parameterizations (https://github.com/OceanBioME/OceanBioME.jl/issues/168) a lot easier because you can put anything in the transfer_velocity that you like as long as it has a method with arguments u, T.

I also changed the default to Ho et al., 2006 (which is the same as Wanninkhof, 2014 - https://x.com/_david_ho_/status/1414754481074589697) but it is important that the coefficients are wind product specific so you should use different values.

[vtamsitt]

Great thanks, I'll go back and take another look at this. I'm not sure if it would be useful to have a way for users to input which wind product is used and then have the correct coefficient chosen from there, or whether to just include all the wind product-specific parameterizations and leave it up to the user.

@vtamsitt the refactoring of the gas exchange model here should make adding different gas transfer parameterizations (#168) a lot easier because you can put anything in the transfer_velocity that you like as long as it has a method with arguments u, T.

I also changed the default to Ho et al., 2006 (which is the same as Wanninkhof, 2014 - https://x.com/_david_ho_/status/1414754481074589697) but it is important that the coefficients are wind product specific so you should use different values.

[jagoosw]

@vtamsitt do you know if there is a list somewhere of the coefficients for each wind product? I can make a function that builds a ScaledTransferVelocity with the appropriate coefficient if the user specifies a wind product.

[vtamsitt]

@vtamsitt do you know if there is a list somewhere of the coefficients for each wind product? I can make a function that builds a ScaledTransferVelocity with the appropriate coefficient if the user specifies a wind product.

The PySeaFlux data product includes the coefficients they have calculated for "ERA5", "JRA55", "NCEP1", "CCMP2" in SeaFlux.v2023.02_alpha_1982-2022.nc. These will presumably be recalculated each time the product is updated.

But these are for Wanninkhof 1992! The code includes a bunch of other gas transfer velocity parameterizations including Ho et al. 2006 and others but as far as I can tell haven't included coefficients for different wind products for the other parameterizations.

[jagoosw]

If I understand correctly, the data product is a climatology of $k_w$ based on the Wanninkhof 1992 parameterization right? But I downloaded it and it doesn't include the parameter value they used so I assume they used the original values maybe.

In their code it says which wind products each parameterization is suitable for (e.g. Ho 2006 scale factor says "The gas transfer velocity is for the QuickSCAT satellite wind product").

Since all of the parameterizations are some form of:$kw = k{660}(u) (Sc/660)^{-1/2}$ where $f(u)$ is some polynomial (possibly piecewise e.g. k_Li86) parameterization I think I'll rename ScaledTransferVelocity to SchmidtScaledTransferVelocity which has Sc and f properties, and then I'll add in the default configurations.

[vtamsitt]

If I understand correctly, the data product is a climatology of kw based on the Wanninkhof 1992 parameterization right? But I downloaded it and it doesn't include the parameter value they used so I assume they used the original values maybe.

In their code it says which wind products each parameterization is suitable for (e.g. Ho 2006 scale factor says "The gas transfer velocity is for the QuickSCAT satellite wind product").

There's a separate file SeaFlux.v2023.02_alpha_1982-2022.nc that contains the coefficient of gas transfer for the four wind products:

netcdf SeaFlux.v2023.02_alpha_1982-2022 {
dimensions:
        wind = 4 ;
variables:
         float alpha(wind) ;
                 alpha:_FillValue = NaNf ;
                 alpha:units = "(cm/hr) (m/s)^-2" ;
                 alpha:long_name = "coefficient of gas transfer" ;
                 string alpha:description = "The coefficient of gas transfer is scaled so that the global average kw is 16.5 cm/hr as recommended by Naegler (2009) for the period 1990-2019 using bomb 14-C estimates. Note that there is 20% uncertainty around the global average estimate of kw (16.5 ± 3.2 cm/hr).Note that during scaling, we calculate the area weighted average for kw over only the ice free ocean (ice frac < 5%) where pixels contain only ocean (sea frac > 95%). We also exclude the Black and Caspian sea from the globally averaged kw. " ;
         string wind(wind) ;
 data:

  alpha = 0.270875, 0.2601975, 0.2866424, 0.256789 ;

  wind = "ERA5", "JRA55", "NCEP1", "CCMP2" ;
 }
 (END) 

Since all of the parameterizations are some form of:$kw = k{660}(u) (Sc/660)^{-1/2}$ where f(u) is some polynomial (possibly piecewise e.g. k_Li86) parameterization I think I'll rename ScaledTransferVelocity to SchmidtScaledTransferVelocity which has Sc and f properties, and then I'll add in the default configurations.

This sounds good, thanks!

[jagoosw]

Ah I see thank you, I've added the different parameterizations now and interestingly they're all very similar except Wanninkhof99 and McGillis01, and you can that Wanninkhof09 is quite differently parameterised from the others but is near by.

[jagoosw]

alpha = 0.270875, 0.2601975, 0.2866424, 0.256789 ;

wind = "ERA5", "JRA55", "NCEP1", "CCMP2" ;

It is interesting that these are all within the error estimates of Ho 2006 (and I can't find errors for other ones).

I've added them now too.

[vtamsitt]

Ah I see thank you, I've added the different parameterizations now and interestingly they're all very similar except Wanninkhof99 and McGillis01, and you can that Wanninkhof09 is quite differently parameterised from the others but is near by.

Thanks for adding all these. This is really interesting, thanks!

[jagoosw]

TODO: check model_implemntation page when rendered

update: it was broken