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.603 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.407 seconds).
[ Info: 2.587 minutes in 52.827 seconds with Δt = 2.148 seconds
[ Info: 7.314 minutes in 1.219 minutes with Δt = 4.381 seconds
[ Info: 16.356 minutes in 1.560 minutes with Δt = 6.163 seconds
[ Info: 28.079 minutes in 1.897 minutes with Δt = 8.240 seconds
[ Info: 41.735 minutes in 2.232 minutes with Δt = 9.062 seconds
[ Info: 56.089 minutes in 2.571 minutes with Δt = 8.043 seconds
[ Info: 1.141 hours in 2.910 minutes with Δt = 7.227 seconds
[ Info: 1.323 hours in 3.246 minutes with Δt = 5.541 seconds
[ Info: 1.476 hours in 3.570 minutes with Δt = 6.128 seconds
[ Info: 1.636 hours in 3.873 minutes with Δt = 5.663 seconds
[ Info: 1.783 hours in 4.170 minutes with Δt = 5.161 seconds
[ Info: 1.916 hours in 4.461 minutes with Δt = 4.519 seconds
[ Info: Simulation is stopping after running for 4.659 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.