Idealised ocean alkalinity enhancement (OAE)

In this example we setup a 3D model driven by a constant wind stress, with 2 instances of the carbon chemistry system (DIC and alkalinity). We then add alkalinity to one of the instances replicating an OAE release, and compare to the unperturbed counterfactual

Install dependencies

First we ensure we have the required dependencies installed

using Pkg
pkg "add OceanBioME, Oceananigans, CairoMakie"

Model setup

using Oceananigans, OceanBioME, Oceananigans.Units

using OceanBioME.Models.GasExchangeModel: CarbonDioxideConcentration

we define the grid and the biogeochemistry, adding in the carbonate_system with 2 replicates

grid = RectilinearGrid(size = (64, 64, 8),
                       extent = (500, 500, 15))

biogeochemistry = LOBSTER(grid;
                          carbonate_system = CarbonateSystem(2))

LOBSTER model (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC1, :DIC2, :Alk1, :Alk2) 
 Light attenuation: Two-band light attenuation model (Float64)
 Sediment: Nothing
 Particles: Nothing
 Modifiers: Nothing

Next we have to construct the boundary conditions for the DIC, changing the defaults to specify that the DIC and Alkalinity for the calculation should be DIC1 and Alk1/DIC2 and Alk2. Then we can put the boundary conditions together.

wind_speed = 5 # m/s
Cᴰ = 2e-3
ρₐ = 1.2
ρₒ = 1026

CO₂_flux1 =
    CarbonDioxideGasExchangeBoundaryCondition(;
        wind_speed,
        water_concentration = CarbonDioxideConcentration(; DIC = :DIC1,
                                                           Alk = :Alk1)
    )
CO₂_flux2 =
    CarbonDioxideGasExchangeBoundaryCondition(;
        wind_speed,
        water_concentration = CarbonDioxideConcentration(; DIC = :DIC2,
                                                           Alk = :Alk2)
    )

wind = FluxBoundaryCondition(-ρₐ/ρₒ * Cᴰ * wind_speed^2)

boundary_conditions = (; DIC1 = FieldBoundaryConditions(top = CO₂_flux1),
                         DIC2 = FieldBoundaryConditions(top = CO₂_flux2),
                         u    = FieldBoundaryConditions(top = wind))

(DIC1 = Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: DiscreteBoundaryFunction with GasExchange
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing), DIC2 = Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: DiscreteBoundaryFunction with GasExchange
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing), u = Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: -5.84795e-5
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing))

Next we define an alkalinity release in a circle in the center for 1 hour

@inline alkalinity_release(x, y, z, t, params) =
    ifelse((params.start_time <= t < params.start_time + params.duration) &
            (z > params.depth) &
            ((x-250)^2 + (y-250)^2 < params.radius^2),
            params.release_rate, 0)

total_release = 2.5e7 # 1t NaOH = 25000mol NaOH = 2.5e7 meq Alk
duration = 1hour
depth = -1
radius = 50
release_rate = total_release/duration/(π*radius^2*abs(depth))

oae = Forcing(alkalinity_release;
              parameters = (; start_time = 20minutes,
                              duration,
                              release_rate,
                              radius,
                              depth))

ContinuousForcing{@NamedTuple{start_time::Float64, duration::Float64, release_rate::Float64, radius::Int64, depth::Int64}}
├── func: alkalinity_release (generic function with 1 method)
├── parameters: (start_time = 1200.0, duration = 3600.0, release_rate = 0.8841941282883075, radius = 50, depth = -1)
└── field dependencies: ()

And then put the model together with Alk2 having the release while Alk1 evolves normally

model = NonhydrostaticModel(grid;
                            coriolis = FPlane(latitude = 45),
                            boundary_conditions,
                            biogeochemistry,
                            advection = WENO(),
                            tracers = (:T, :S),
                            forcing = (; Alk2 = oae))

set!(model, P = 1, NO₃ = 10,
            T = 15, S = 35,
            DIC1 = 2000, Alk1 = 2300,
            DIC2 = 2000, Alk2 = 2300,
            u = (x, y, z) -> randn()/2)

Then we can setup a simulation with a variable timestep and progress message and output the tracers

simulation = Simulation(model, Δt = 5, stop_time = 2hours)

conjure_time_step_wizard!(simulation)

prog(sim) =
    @info prettytime(sim)*" in "*prettytime(sim.run_wall_time)*" with Δt = "*prettytime(sim.Δt)

add_callback!(simulation, prog, IterationInterval(100))

simulation.output_writers[:tracers] = JLD2Writer(model, model.tracers;
                                                filename = "oae.jld2",
                                                schedule = AveragedTimeInterval(5minutes),
                                                overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(5 minutes):
├── filepath: oae.jld2
├── 13 outputs: (T, S, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC1, DIC2, Alk1, Alk2) averaged on AveragedTimeInterval(window=5 minutes, stride=1, interval=5 minutes)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

We can also create a BoundaryConditionOperation which records the flux through the top boundary for both the DIC fields which we can save

qCO₂1 = BoundaryConditionOperation(model.tracers.DIC1, :top, model)
qCO₂2 = BoundaryConditionOperation(model.tracers.DIC2, :top, model)


simulation.output_writers[:carbon_flux] = JLD2Writer(model, (; qCO₂1, qCO₂2),
                                                     indices = (:, :, grid.Nz),
                                                     filename = "oae_surface_flux.jld2",
                                                     schedule = TimeInterval(5minutes),
                                                     overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(5 minutes):
├── filepath: oae_surface_flux.jld2
├── 2 outputs: (qCO₂1, qCO₂2)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

and then run the simulation

run!(simulation)

[ Info: Initializing simulation...
[ Info: 0 seconds in 0 seconds with Δt = 2.500 seconds
[ Info:     ... simulation initialization complete (28.763 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.420 seconds).
[ Info: 2.422 minutes in 52.825 seconds with Δt = 1.996 seconds
[ Info: 6.364 minutes in 1.207 minutes with Δt = 3.321 seconds
[ Info: 14.070 minutes in 1.525 minutes with Δt = 5.931 seconds
[ Info: 25.132 minutes in 1.850 minutes with Δt = 8.104 seconds
[ Info: 39.373 minutes in 2.180 minutes with Δt = 8.102 seconds
[ Info: 52.107 minutes in 2.503 minutes with Δt = 7.188 seconds
[ Info: 1.069 hours in 2.834 minutes with Δt = 7.398 seconds
[ Info: 1.250 hours in 3.164 minutes with Δt = 6.384 seconds
[ Info: 1.425 hours in 3.469 minutes with Δt = 5.863 seconds
[ Info: 1.575 hours in 3.755 minutes with Δt = 4.710 seconds
[ Info: 1.717 hours in 4.041 minutes with Δt = 5.286 seconds
[ Info: 1.852 hours in 4.327 minutes with Δt = 4.733 seconds
[ Info: 1.980 hours in 4.613 minutes with Δt = 4.461 seconds
[ Info: Simulation is stopping after running for 4.662 minutes.
[ Info: Simulation time 2 hours equals or exceeds stop time 2 hours.

Finally we load the data and create some plots

using CairoMakie, Statistics

fds = FieldDataset("oae.jld2")
fds_surface = FieldDataset("oae_surface_flux.jld2")

n = Observable(1)

N = length(fds["P"])

Δ_Alk = @lift (mean(interior(fds["Alk2"][$n]), dims=3)[:, :, 1] .-
               mean(interior(fds["Alk1"][$n]), dims=3)[:, :, 1]).* 15
Δ_qCO₂ = @lift (interior(fds_surface["qCO₂2"][$n], :, :, 1) .-
                interior(fds_surface["qCO₂1"][$n], :, :, 1)) .* (12+2*16)*1e-3*day

Δ_Alk_range = maximum(abs, mean(interior(fds["Alk2"]), dims = 3) .-
                           mean(interior(fds["Alk1"]), dims = 3)) .* 15
Δ_qCO₂_range = maximum(abs, interior(fds_surface["qCO₂2"]) .-
                            interior(fds_surface["qCO₂1"])) .* (12+2*16)*1e-3*day

xc = xnodes(fds["P"])
yc = ynodes(fds["P"])

fig = Figure(size=(1000, 500))

ax = Axis(fig[1, 1], aspect = DataAspect())
ax2 = Axis(fig[1, 2], aspect = DataAspect())

hm = heatmap!(ax, xc, yc, Δ_Alk, colormap = :balance, colorrange = (-1, 1) .* Δ_Alk_range)
hm2 = heatmap!(ax2, xc, yc, Δ_qCO₂, colormap = :balance, colorrange = (-1, 1) .* Δ_qCO₂_range)

Colorbar(fig[2, 1], hm, vertical = false,
                        label = "Alkalinity pertubation (mmol / m³)",
                        flip_vertical_label = true)
Colorbar(fig[2, 2], hm2, vertical = false,
                         label = "Surface CO₂ flux difference (gCO₂/m²/day)",
                         flip_vertical_label = true)

supertitle = Label(fig[0, :], (@lift prettytime(fds["P"].times[$n])))

CairoMakie.record(fig, "oae.mp4", 1:N, framerate=10) do i
    @info "$i"
    n[] = i
end

"oae.mp4"

fig = Figure();

ax = Axis(fig[1, 1], title = "Surface flux (gC/day)")

surface_area = 500*500
mmol_to_g = 12/1000

total_carbon_flux1 =  mean(interior(fds_surface["qCO₂1"], :, :, 1, :), dims = (1, 2))[1, 1, :]
total_carbon_flux2 =  mean(interior(fds_surface["qCO₂2"], :, :, 1, :), dims = (1, 2))[1, 1, :]

lines!(ax, fds_surface["qCO₂1"].times./hours, total_carbon_flux1 .* mmol_to_g * surface_area * day)
lines!(ax, fds_surface["qCO₂2"].times./hours, total_carbon_flux2 .* mmol_to_g * surface_area * day)

ax2 = Axis(fig[2, 1], xlabel = "Time (hours)", title = "Cumulative difference (gC)")

Δt = diff(fds_surface["qCO₂1"].times./hours)[1]

change =  mean(interior(fds_surface["qCO₂1"], :, :, 1, :), dims = (1, 2))[1, 1, :] .-
          mean(interior(fds_surface["qCO₂2"], :, :, 1, :), dims = (1, 2))[1, 1, :]
change .*= mmol_to_g * surface_area

Δt_save = diff(fds_surface["qCO₂1"].times)[1]

cumulative_change = cumsum(Δt_save .* change)

lines!(ax2, fds_surface["qCO₂1"].times./hours, cumulative_change)

fig

This page was generated using Literate.jl.