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 inorganic_carbon with 2 replicates

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

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

NutrientsPlanktonDetritus{Float64} with (:NO₃, :NH₄, :P, :Z, :DOM, :sPOM, :bPOM, :DIC1, :DIC2, :Alk1, :Alk2) 
 Light attenuation: TwoBandPhotosyntheticallyActiveRadiation{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, DOM, sPOM, bPOM, 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.295 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.211 seconds).
[ Info: 2.428 minutes in 1.007 minutes with Δt = 1.866 seconds
[ Info: 6.912 minutes in 1.478 minutes with Δt = 3.413 seconds
[ Info: 14.976 minutes in 1.949 minutes with Δt = 6.745 seconds
[ Info: 27.486 minutes in 2.422 minutes with Δt = 8.613 seconds
[ Info: 41.930 minutes in 2.892 minutes with Δt = 8.759 seconds
[ Info: 55.259 minutes in 3.364 minutes with Δt = 8.320 seconds
[ Info: 1.132 hours in 3.824 minutes with Δt = 7.293 seconds
[ Info: 1.321 hours in 4.288 minutes with Δt = 6.267 seconds
[ Info: 1.487 hours in 4.751 minutes with Δt = 5.473 seconds
[ Info: 1.625 hours in 5.212 minutes with Δt = 5.081 seconds
[ Info: 1.766 hours in 5.663 minutes with Δt = 4.731 seconds
[ Info: 1.897 hours in 6.092 minutes with Δt = 4.742 seconds
[ Info: Simulation is stopping after running for 6.457 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.