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 (29.061 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.468 seconds).
[ Info: 2.497 minutes in 53.353 seconds with Δt = 1.999 seconds
[ Info: 6.969 minutes in 1.226 minutes with Δt = 3.946 seconds
[ Info: 15.471 minutes in 1.562 minutes with Δt = 5.843 seconds
[ Info: 27.369 minutes in 1.893 minutes with Δt = 9.077 seconds
[ Info: 42.191 minutes in 2.225 minutes with Δt = 9.513 seconds
[ Info: 56.632 minutes in 2.561 minutes with Δt = 8.378 seconds
[ Info: 1.158 hours in 2.897 minutes with Δt = 6.548 seconds
[ Info: 1.335 hours in 3.231 minutes with Δt = 6.208 seconds
[ Info: 1.496 hours in 3.563 minutes with Δt = 5.758 seconds
[ Info: 1.650 hours in 3.891 minutes with Δt = 5.423 seconds
[ Info: 1.791 hours in 4.225 minutes with Δt = 5.094 seconds
[ Info: 1.922 hours in 4.556 minutes with Δt = 4.612 seconds
[ Info: Simulation is stopping after running for 4.757 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.