Non-zero CO2 surface flux for an ocean in equilibrium (with non-zero average wind speed)

[ali-ramadhan]

Thank you everyone for the super nice package (especially @jagoosw)! I know everything is still in development with lots of changes happening.

I've been playing around with LOBSTER (well, really just the carbonate chemistry component) and the air-sea gas exchange model. I'm just initializing an ocean at equilibrium with constant T, S, DIC, and Alk. I set the atmospheric pCO2 equal to the ocean's pCO2. I expect a surface CO2 flux of zero. I see this but only when average_wind_speed = 0.

When average_wind_speed = 10.0 (the default) I believe the flux should still be zero since pCO2 is equal across the air-ocean interface but it's slightly negative after one time step. It's zero at time step 0.

Is this behavior expected or am I misunderstanding the gas exchange model?

Thank you!

---

Here's a minimal working example to showcase what I'm seeing:

using Oceananigans
using OceanBioME

using OceanBioME.Boundaries: OCMIP_default

const pCO₂_OCMIP_default = OCMIP_default

grid = RectilinearGrid(
    CPU(),
    topology = (Periodic, Flat, Bounded),
    size = (256, 128),
    x = (0, 200),
    z = (-100, 0),
    halo = (4, 4)
)

ambient_Alk = 2300.0
ambient_DIC = 2100.0
SST = 20.0
SSS = 35.0

# Let's set the atmospheric pCO₂ to equal the oceanic pCO₂.
pCO₂_atmosphere = pCO₂_OCMIP_default(ambient_DIC, ambient_Alk, SST, SSS)

CO₂_flux = GasExchange(;
    gas = :CO₂,
    air_concentration = pCO₂_atmosphere,
    average_wind_speed = 10.0
)

boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), )

model = NonhydrostaticModel(;
    grid,
    buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState()),
    biogeochemistry = LOBSTER(; grid, carbonates=true),
    tracers = (:T, :S, :NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk),
    boundary_conditions = boundary_conditions
)

set!(model, T=SST, S=SSS, Alk=ambient_Alk, DIC=ambient_DIC)

function surface_CO2_flux(model)
    grid = model.grid
    T_surface = interior(model.tracers.T, :, 1, grid.Nz)
    S_surface = interior(model.tracers.S, :, 1, grid.Nz)
    DIC_surface = interior(model.tracers.DIC, :, 1, grid.Nz)
    Alk_surface = interior(model.tracers.Alk, :, 1, grid.Nz)
    flux_func = model.tracers.DIC.boundary_conditions.top.condition.func
    return flux_func.(0, 0, 0, DIC_surface, Alk_surface, T_surface, S_surface)
end

Time step 0:

julia> surface_CO2_flux(model)
256-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0

Time step 1:

julia> time_step!(model, 1)

julia> surface_CO2_flux(model)
256-element Vector{Float64}:
 -8.941545089802037e-8
 -8.941545089802037e-8
 -8.941545089802037e-8
 -8.941545089802037e-8
 -8.941545089802037e-8
  ⋮
 -8.941545089802037e-8
 -8.941545089802037e-8
 -8.941545089802037e-8
 -8.941545089802037e-8
 -8.941545089802037e-8

---

I'm using Julia 1.10.4 with

  [a49af516] OceanBioME v0.10.5
  [9e8cae18] Oceananigans v0.91.4

[jagoosw]

Hi @ali-ramadhan, thank you for raising this issue!

I've had a look and am pretty sure this is because the grid is flat, so when Oceananigans calls getbc the arguments are x, t, DIC, Alk, T, S, but when I wrote this it would have passed x, 0, t, DIC, Alk, T, S, so its not computing the pCO2 in the water correctly first. I.e. its using the second method not the first one here:

https://github.com/OceanBioME/OceanBioME.jl/blob/4bcaf29797abcb7c7293a299c2e1c549e3d3d7eb/src/Boundaries/gasexchange.jl#L223-L232

I tried running your code with y periodic and it appears to be correct.

Currently, I'm doing a lot of work on the gas exchange model because I've discovered a few issues with how I originally implemented it, and following that PR (#186) it will be more sophisticated and better tested. This situation would also cause an issue with the refactored code so I'll work out how to fix it there.

For now, the best thing is probably to only use 3D grids I think.

[ali-ramadhan]

Thanks for having a look @jagoosw! Tricky dispatch haha. It does work with a 3D grid and I was able to hack in a boundary condition wrapper for 2D grids.

gas_exchange = GasExchange(;
    gas = :CO₂,
    air_concentration = pCO₂_atmosphere
)

@inline CO₂_flux(x, t, DIC, Alk, T, S, ge) = ge.condition.func(x, 0, t, DIC, Alk, T, S)

DIC_top_bc = FluxBoundaryCondition(CO₂_flux, field_dependencies = (:DIC, :Alk, :T, :S), parameters = gas_exchange)

Out of curiousity, do you know why it works correctly when average_wind_speed = 0.0?

[glwagner]

Just to chime in --- I recommend using the discrete form in packages that build out functionality for Oceananigans. That should avoid the issues that occur when trying to write generic code for Flat or 3D grids.

[jagoosw]

Out of curiousity, do you know why it works correctly when average_wind_speed = 0.0?

The flux is like $k(u{wind})(pCO{2, \text{air}} - pCO{2, \text{water}})$ and the common form of $k(u{wind})$ is $k0u{wind}^2(660/Sc)^{1/2}$ so it also goes to zero when the wind speed is zero. The new implementation allows the actual wind speed to be put in, as well as having different transfer velocities since a few have a constant factor was well.