How to simulate only CarbonChemistry and GasExchange

[lgloege]

I would like to run some idealized simulations with just carbon chemistry and air-sea gas exchange. Ideally, I would like to be able to do something like this:

using OceanBioME
using Oceananigans

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200));

biogeochemistry = CarbonChemistry(; grid) # also want GasExchange

CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()
boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),)
tracers = (:T, :S, :DIC, :ALK)

model = NonhydrostaticModel(; grid, biogeochemistry , boundary_conditions, tracers)

This obviously doesn't work. A somewhat hacky solution would be to to use the NPZD model and set P and Z to zero. Curious if there is something more elegant than this so I don't have to carry around P, Z and D tracers unnecessarily.

Please let me know if there is a way to achieve this. Thanks

[jagoosw]

Hi @lgloege,

The easiest way would probably be to manually define the tracers like you have done already, and just omit the biogeochemistry since the carbon chemistry lives in the gas exchange boundary condition anyway. This seems to work for me:

using OceanBioME
using Oceananigans

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200));

CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()
boundary_conditions = (; DIC = FieldBoundaryConditions(top = CO₂_flux))
tracers = (:T, :S, :DIC, :Alk) # note Alk not ALK

model = NonhydrostaticModel(; grid, boundary_conditions, tracers)

set!(model, DIC = 2200, Alk = 2000)

time_step!(model, 1000.0)

model.tracers.DIC
3×3×30 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 3×3×30 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: Flux, immersed: ZeroFlux
└── data: 9×9×36 OffsetArray(::Array{Float64, 3}, -2:6, -2:6, -2:33) with eltype Float64 with indices -2:6×-2:6×-2:33
    └── max=2200.0, min=2199.44, mean=2199.98

A next step if you wanted to specify something about what the DIC and alkalinity do could be to define a new bgc model with just DIC and Alkalinity then you won't have to specify the tracers too.

[lgloege]

Thank! @jagoosw ! I see that CarbonChemistry() is called in the CarbonDioxideGasExchangeBoundaryCondition(). Not obvious to me that carbon chem is being calculated at every grid cell because of this setup. Maybe i'm missing something, but it feels like you should setup the carbon chemistry in the model definition? This is why I thought defining it in biogeochemistry sort of made sense (although there is no bio if you're just defining the chemistry).

I'm getting closer to what I am trying to do. Now my issue is returning pH. My naive solution was to pass CarbonChemistry() to CarbonDioxideGasExchangeBoundaryCondition(), but I don't think that's the correct approach.

Nx, Nz = 256, 128
H, L = 2kilometers, 1000kilometers

grid = RectilinearGrid(size = (Nx, Nz), 
                       halo = (4, 4),
                       x = (-L, L), 
                       z = (-H, 0),
                       topology = (Bounded, Flat, Bounded))

coriolis = FPlane(latitude = -45)
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = 2e-4,
                                                                    haline_contraction = 8e-4))

# =====================================
# wind stress
# =====================================
u₁₀ = 2     # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3 # dimensionless drag coefficient
ρₐ = 1.225  # kg m⁻³, average density of air at sea-level
ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean

τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))


# =====================================
# CO2 gas transfer
# =====================================
carbon_chemistry = CarbonChemistry(return_pH = true)
CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition(carbon_chemistry, air_concentration = 425, wind_speed=u₁₀)
DIC_bcs = FieldBoundaryConditions(top = CO₂_flux)

# =====================================
# Alk forcing
# add an Alk Forcing() here
# to be implemented
# =====================================


tracers = (:T, :S, :DIC, :Alk, :pH) 

model = NonhydrostaticModel(; grid, buoyancy,
                            advection = UpwindBiased(order=5),
                            tracers,
                            coriolis,
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs, DIC=DIC_bcs))

set!(model, T = 20, S = 35, DIC = 2200, Alk = 2000)

simulation = Simulation(model, Δt=10.0, stop_time=40minutes)

# Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n",
                                iteration(sim), prettytime(sim), prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))

add_callback!(simulation, progress_message, IterationInterval(20))

filename = "test_sim"

simulation.output_writers[:slices] =
    JLD2Writer(model, merge(model.velocities, model.tracers),
               filename = filename * ".jld2",
               indices = (:, :),
               schedule = TimeInterval(1minute),
               overwrite_existing = true)

run!(simulation)

[jagoosw]

Thank! @jagoosw ! I see that CarbonChemistry() is called in the CarbonDioxideGasExchangeBoundaryCondition(). Not obvious to me that carbon chem is being calculated at every grid cell because of this setup. Maybe i'm missing something, but it feels like you should setup the carbon chemistry in the model definition? This is why I thought defining it in biogeochemistry sort of made sense (although there is no bio if you're just defining the chemistry).

That makes sense and we're not set on the current architecture if it would be clearer to restructure a bit. The carbon chemistry isn't actually computed at every grid point in the normal configuration because all it does is diagnose pCO2 or pH, which is only usually used in the surface flux. The exception (in OceanBioME, but I'm sure there are others) is that you have to pass a carbon chemistry to PISCES as it uses it to compute the calcite saturation:

OceanBioME.jl/src/Models/AdvectedPopulations/PISCES/PISCES.jl

Lines 323 to 326 in 86bbfbc

	carbon_chemistry = CarbonChemistry(FT),
	calcite_saturation = CenterField(grid),

I'm getting closer to what I am trying to do. Now my issue is returning pH. My naive solution was to pass CarbonChemistry() to CarbonDioxideGasExchangeBoundaryCondition(), but I don't think that's the correct approach.

Since the carbon chemistry is diagnostic you can't make a pH tracer, a way you could make a field of the pH would be to write a kernel function operation. Using your code above:

...
tracers = (:T, :S, :DIC, :Alk, :pH) 

model = NonhydrostaticModel(; grid, buoyancy,
                            advection = UpwindBiased(order=5),
                            tracers,
                            coriolis,
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs, DIC=DIC_bcs))

set!(model, T = 20, S = 35, DIC = 2200, Alk = 2000)

DIC, Alk, T, S = model.tracers[(:DIC, :Alk, :T, :S)]

@inline function pH_kernel(i, j, k, grid, DIC, Alk, T, S, carbon_chemistry)
    @inbounds begin
        DIC = DIC[i, j, k]
        Alk = Alk[i, j, k]
        T = T[i, j, k]
        S = S[i, j, k]
    end

    return carbon_chemistry(; DIC, Alk, T, S, return_pH = true)
end

pH_op = KernelFunctionOperation{Center, Center, Center}(pH_kernel, grid, DIC, Alk, T, S, carbon_chemistry)

pH = Field(pH_op)

compute!(pH)

pH[1, 1, 1]
6.845136533754769

[jagoosw]

When you output fields with kernel function operations inside they get computed automatically

[lgloege]

This is maybe a separate discussion, but where i'm going with this is calculating uptake from an OAE release, very idealized right now. To do this I would need to difference two simulations. One with and one without additional alkalinity. I'm curious how challenging it would be calculate two carbonate systems in the same simulation. One that is a baseline "counterfactual" simulation and another that includes the additional alkalinity. Benefit is you only have to run the simulation once. I believe this is implemented in MARBL, but I don't see why it can't be generalized to the entire OceanBioME framework so it would work with all ecosystem models. I can open this up as a separate discussion if this is useful. I don't think my knowledge of Julia of or the codebase is there yet for me to make this contribution. Having that said, i'm happy to assist somebody is there is interest in this! I really want Oceanigans+OceanBioME to be the framework for ocean-based CDR modeling.

[jagoosw]

I guess you can open another discussion if/when you start implementing this but it seems like it would be relatively straight forward, you would just need to add a second DIC and alkalinity tracers which evolve the same (since nothing else depends on DIC and alkalinity in most models) and just do the OAE to one of the versions.

If you want to look into this, I can have a think about how best to restructure e.g. LOBSTER to allow this (since I think LOBSTER could do with tidying up anyway).

[lgloege]

Thanks @jagoosw , so now the model is calculating the sqrt of a negative somewhere. This only happens after I turn on gas exchange. Here is the main part of the error message:

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (49.033 ms)
[ Info: Executing initial time step...
DomainError with -7.689605850147218e-10:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] sqrt
    @ ./math.jl:608 [inlined]
  [3] KB
    @ ~/.julia/packages/OceanBioME/aMtdF/src/Models/CarbonChemistry/equilibrium_constants.jl:243 [inlined]
  [4] #_#37
    @ ~/.julia/packages/OceanBioME/aMtdF/src/Models/CarbonChemistry/carbon_chemistry.jl:177 [inlined]
  [5] CarbonChemistry
    @ ~/.julia/packages/OceanBioME/aMtdF/src/Models/CarbonChemistry/carbon_chemistry.jl:145 [inlined]
  [6] surface_value
    @ ~/.julia/packages/OceanBioME/aMtdF/src/Models/GasExchange/carbon_dioxide_concentration.jl:42 [inlined]
  [7] GasExchange
    @ ~/.julia/packages/OceanBioME/aMtdF/src/Models/GasExchange/gas_exchange.jl:35 [inlined]
  [8] getbc
    @ ~/.julia/packages/Oceananigans/Xawvn/src/BoundaryConditions/discrete_boundary_function.jl:36 [inlined]
  [9] getbc
    @ ~/.julia/packages/Oceananigans/Xawvn/src/BoundaryConditions/boundary_condition.jl:112 [inlined]
 [10] apply_z_top_bc!
    @ ~/.julia/packages/Oceananigans/Xawvn/src/BoundaryConditions/apply_flux_bcs.jl:161 [inlined]
 [11] macro expansion
    @ ~/.julia/packages/Oceananigans/Xawvn/src/BoundaryConditions/apply_flux_bcs.jl:82 [inlined]
 [12] cpu__apply_z_bcs!
    @ ~/.julia/packages/KernelAbstractions/C3nYQ/src/macros.jl:304 [inlined]

and here is my code

using CairoMakie
using Dates
using NCDatasets
using Oceananigans
using OceanBioME
using Oceananigans.Units: meters, minute

# =======================================
# define architecture CPU or GPU
# =======================================
architecture = CPU()

# =======================================
# define model grid
# =======================================
x = (0, 15meters)
z = (-0.5meters, 0)
Δx, Δz = 0.01meters, 0.01meters
size = (Int(x[2]/Δx), Int(-z[1]/Δz))
topology = (Bounded, Flat, Bounded)

grid = RectilinearGrid(architecture; size, x, z, topology)

# =====================================
# wind stress
# =====================================
u₁₀ = 5     # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3 # dimensionless drag coefficient
ρₐ = 1.225  # kg m⁻³, average density of air at sea-level
ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean

τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))

# =====================================
# CO2 gas transfer
# =====================================
CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition(;air_concentration = 425, wind_speed=u₁₀)
DIC_bcs = FieldBoundaryConditions(top = CO₂_flux)

# =======================================
# Model paramters
# =======================================
# closure = ScalarDiffusivity(ν=1e-6, κ=1e-6)
closure = AnisotropicMinimumDissipation() 
tracers = (:T, :S, :c, :DIC, :Alk)
coriolis = FPlane(f=0.2)
advection= UpwindBiased(order=5) #WENO()

boundary_conditions = (u=u_bcs, DIC=DIC_bcs)

# ---------------------------------------
# equation of state and buoyancy
# b = g (- β S - γ C) # ignore temperature influence
# ---------------------------------------
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = 2e-4,
                                                                    haline_contraction = 8e-4))

model = NonhydrostaticModel(; grid, 
                              buoyancy,
                              tracers,
                              coriolis, 
                              advection, 
                              closure,
                              boundary_conditions)

# =======================================
# set intial salinity and tracers
# =======================================
# function to get index closest to value
closest_index(target, collection) = argmin(abs.(collection .- target))

# create x and z collections
xᶜ = collect(grid.xᶜᵃᵃ)
zᶜ = collect(grid.z.cᵃᵃᶜ)

# find index range for x ∈ [0, 0.15] 
i₁ = closest_index(0.0,   xᶜ)
i₂ = closest_index(0.15,  xᶜ)

# find index range for z ∈ [-0.15, 0]
k₁ = closest_index(-0.15, zᶜ)
k₂ = closest_index(0.0,   zᶜ)

# create zero array on center grid
S_initial = CenterField(grid)
c_initial = CenterField(grid)

# set values within the index range
S_initial[i₁:i₂, :, k₁:k₂] .= 25.0;
c_initial[i₁:i₂, :, k₁:k₂] .= 0.0078;

# now set the values in the model
set!(model, T=25, S=S_initial, c = c_initial, DIC = 2200, Alk = 2000)

simulation = Simulation(model; Δt = 0.005, stop_iteration=10)

# --------------------------------------------------
# output writer
# --------------------------------------------------
#filename = "test-release_$(string(today()))"
filename = "./test_sim"

fields = Dict("T" => model.tracers.T, "S" => model.tracers.S, "c" => model.tracers.c, "DIC" => model.tracers.DIC, "Alk" => model.tracers.Alk)

simulation.output_writers[:tracers] = JLD2Writer(model, fields; 
                                                 filename = filename * ".jld2",
                                                 schedule = IterationInterval(20),
                                                 overwrite_existing = true)

run!(simulation)

simulation.stop_iteration = 100
simulation.Δt = 0.05
run!(simulation)

i'm at a loss what the issue could be, any ideas on how to diagnose the issue?

[jagoosw]

Errors like this tend to arise when the model gets a value for T/S/DIC/Alk significantly out of the normal range. Looking at the specific equilibrium constant, it looks like it is only S that gets rooted:

OceanBioME.jl/src/Models/CarbonChemistry/equilibrium_constants.jl

Lines 243 to 250 in 5dafe8f

	@inline (c::KB)(T::FT, S; P = nothing) where FT =
	c.pressure_correction(T, P) *
	exp(c.constant
	+ (c.inverse_T + c.inverse_T_sqrt_S * √S + c.inverse_T_S * S + c.inverse_T_sqrt_S³ * S^convert(FT, 1.5) + c.inverse_T_S² * S^convert(FT, 2)) / T
	+ c.sqrt_S * √S
	+ c.S * S
	+ (c.log_T + c.log_T_sqrt_S * √S + c.S_log_T * S ) * log(T)
	+ c.T_sqrt_S * √S * T)

I suspect that the timestep might not be small enough for the physics (to keep the CFL below 0.5 with this timestep the speed has to be less than 0.1m/s which seems unlikely with a 5m/s wind), so perhaps just trying a smaller timestep would work?

[lgloege]

Thanks @jagoosw , this helped a lot. I was confused where the sqrt was being applied. Time step wasn't the issue. I am simulating a freshwater plume and the model did not like having a salinity of zero in a small region and values of 25 PSU everywhere else. My solution was to set the salinity to. small value in that region:

# set values within the index range
S_initial[:, :, :] .= 0.1;
S_initial[i₁:i₂, :, k₁:k₂] .= 25;
c_initial[i₁:i₂, :, k₁:k₂] .= 0.0078;

# now set the values in the model
set!(model, T=25, S=S_initial, c = c_initial, DIC = 2200, Alk = 2000)

But I can have salinity be 0 everywhere and the model does not complain. Thank you for helping solve the issue. Any idea why I need to set salinity to a small value? Why can't I have it be zero in a small region?

[jagoosw]

I guess maybe it's because when it's exactly zero, the numerical error in the advection/diffusion can make it go slightly below zero (e.g. 10^-10 above), which the carbon chemistry can't cope with. I don't recall if progress has been made on this front but there might be positivity preserving advection and diffusion in Oceananigans now which would prevent this, alternatively a hack that we have implemented (and is very common in operational bgc models) is to either clip values below zero (ZeroNegativeTracers) or for groups of conserved tracers scale them (ScaleNegativeTracers, although that wouldn't be appropriate in this case).