Illustration of using a carbon system solver
This example illustrates how to use ClimaOceanBiogeochemistry's UniversalRobustCarbonSystem model in a 0-d context.
using ClimaOceanBiogeochemistry.CarbonSystemSolvers.UniversalRobustCarbonSolver: UniversalRobustCarbonSystem
using ClimaOceanBiogeochemistry.CarbonSystemSolvers: CarbonSystemParameters, CarbonSolverParameters, CarbonCoefficientParameters
using Oceananigans
using Oceananigans.Units
using Printf
using CairoMakieWe are going to simulate two bottles of soda, one opened and left in the fridge the other opened and left on the counter to go flat.
Model setup
The 0-d grid represents a 10cm bottle of soda
grid = RectilinearGrid(size = 1,
z = (-1meter/10, 0),
topology = (Flat, Flat, Bounded))1×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×1 halo
├── Flat x
├── Flat y
└── Bounded z ∈ [-0.1, 0.0] regularly spaced with Δz=0.1For CO₂ exchange, we use a FluxBoundaryCondition for the "top" of the dissolved inorganic carbon (DIC) tracer. We'll write a callback later to calculate the flux every time step.
co₂_flux = Field{Center, Center, Nothing}(grid)
dic_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(co₂_flux))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: 1×1×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Field arrays to store pCO₂ values filled in compute_co₂_flux!
soda_pCO₂ = Field{Center, Center, Nothing}(grid)
atmos_pCO₂ = Field{Center, Center, Nothing}(grid)1×1×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU
├── grid: 1×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×1 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 1×1×1 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, 1:1) with eltype Float64 with indices 1:1×1:1×1:1
└── max=0.0, min=0.0, mean=0.0Build function for CO₂ flux calculation. Dissolved CO₂ in the soda will exchange with the overlying atmosphere These are some coefficients and constants that we'll use in the calculation
Base.@kwdef struct CO₂_flux_parameters{FT}
surface_wind_speed :: FT = 10. # ms⁻¹
atmospheric_pCO₂ :: FT = 280e-6 # atm
exchange_coefficient :: FT = 0.337 # cm hr⁻¹
salinity :: FT = 0.0 # psu
alkalinity :: FT = 2.35e-3 # mol kg⁻¹
silicate :: FT = 0e-6 # mol kg⁻¹
phosphate :: FT = 0e-6 # mol kg⁻¹
initial_pH_guess :: FT = 8.0
reference_density :: FT = 1024.5 # kg m⁻³
end
"""
compute_co₂_flux!(simulation)
Returns the tendency due to CO₂ flux using the piston velocity
formulation of Wanninkhof (1992) and the solubility/activity of
CO₂ that depends on temperature, etc.
"""
@inline function compute_co₂_flux!(simulation; solver_params = ())
Nz = size(simulation.model.grid, 3)
# Get coefficients from CO₂_flux_parameters struct
# I really want the option to take these from the model
(; surface_wind_speed,
atmospheric_pCO₂,
exchange_coefficient,
salinity,
alkalinity,
silicate,
phosphate,
initial_pH_guess,
reference_density,
) = CO₂_flux_parameters()
U₁₀ = surface_wind_speed
pCO₂ᵃᵗᵐ = atmospheric_pCO₂
Kʷₐᵥₑ = exchange_coefficient
Sᴬ = salinity
Aᵀ = alkalinity
Siᵀ = silicate
Pᵀ = phosphate
pH = initial_pH_guess
ρʳᵉᶠ = reference_density
cmhr⁻¹_per_ms⁻¹ = 1 / 3.6e5 # conversion factor from cm/hr to m/s
co₂_flux = simulation.model.tracers.DIC.boundary_conditions.top.condition
Θᶜ = simulation.model.tracers.T[1,1,Nz]
Cᵀ = simulation.model.tracers.DIC[1,1,Nz]/ρʳᵉᶠ
# applied pressure in atm (i.e. Pˢᵘʳᶠ-Pᵃᵗᵐ)
# Positive when the can is sealed, then zero when the can is opens
# On average, the 12 ounce soda cans sold in the US tend to have a pressure of roughly 120 kPa when canned at 4 °C
if iteration(simulation) <= 1
Δpᵦₐᵣ = 0.2
else
# *pssshhhht* the bottle is opened
Δpᵦₐᵣ = 0.0
end
# compute soda pCO₂ using the UniversalRobustCarbonSystem solver
# Returns soda pCO₂ (in atm) and atmosphere/soda solubility coefficients (mol kg⁻¹ atm⁻¹)
(; pCO₂ᵒᶜᵉ, Pᵈⁱᶜₖₛₒₗₐ, Pᵈⁱᶜₖₛₒₗₒ) = UniversalRobustCarbonSystem(
pH = pH,
pCO₂ᵃᵗᵐ = pCO₂ᵃᵗᵐ,
Θᶜ = Θᶜ,
Sᴬ = Sᴬ,
Δpᵦₐᵣ = Δpᵦₐᵣ,
Cᵀ = Cᵀ,
Aᵀ = Aᵀ,
Pᵀ = Pᵀ,
Siᵀ = Siᵀ,
solver_params...,
)
# store the soda and atmospheric CO₂ concentrations into Fields
soda_pCO₂[1,1,Nz] = (pCO₂ᵒᶜᵉ * Pᵈⁱᶜₖₛₒₗₒ ) * ρʳᵉᶠ # Convert mol kg⁻¹ m s⁻¹ to mol m⁻² s⁻¹
atmos_pCO₂[1,1,Nz] = (pCO₂ᵃᵗᵐ * Pᵈⁱᶜₖₛₒₗₐ) * ρʳᵉᶠ # Convert mol kg⁻¹ to mol m⁻³
# compute schmidt number for DIC
kˢᶜ = CarbonCoefficientParameters(
a₀ = 2116.8,
a₁ = 136.25,
a₂ = 4.7353,
a₃ = 9.2307e-2,
a₄ = 7.555e-4,
a₅ = 660.0,
)
Cˢᶜᵈⁱᶜ = ( kˢᶜ.a₀ -
kˢᶜ.a₁ * Θᶜ +
kˢᶜ.a₂ * Θᶜ^2 -
kˢᶜ.a₃ * Θᶜ^3 +
kˢᶜ.a₄ * Θᶜ^4
)/kˢᶜ.a₅
# compute gas exchange coefficient/piston velocity and correct with Schmidt number
Kʷ = (
(Kʷₐᵥₑ * cmhr⁻¹_per_ms⁻¹) * U₁₀^2
) / sqrt(Cˢᶜᵈⁱᶜ)
# compute co₂ flux (-ve for uptake, +ve for outgassing since convention is +ve upwards in the soda)
co₂_flux[1,1,Nz] = - Kʷ * (
pCO₂ᵃᵗᵐ * Pᵈⁱᶜₖₛₒₗₐ -
pCO₂ᵒᶜᵉ * Pᵈⁱᶜₖₛₒₗₒ
) * ρʳᵉᶠ # Convert mol kg⁻¹ m s⁻¹ to mol m⁻² s⁻¹
return nothing
endMain.compute_co₂_flux!Simulation of a soda outgassing CO₂ in the fridge
model_open_in_the_fridge = NonhydrostaticModel(;
grid,
velocities = nothing,
buoyancy = nothing,
closure = nothing,
tracers = (:T, :DIC),
boundary_conditions = (; DIC=dic_bcs),
)NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×1 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── advection scheme: Centered reconstruction order 2
├── tracers: (T, DIC)
├── closure: Nothing
├── buoyancy: Nothing
└── coriolis: NothingInitial conditions for the refridgerated soda
Tᵢ = 4 # °C
DICᵢ = 2.4 # mol/m³
set!(model_open_in_the_fridge, T = Tᵢ, DIC = DICᵢ)
simulation = Simulation(model_open_in_the_fridge, Δt=10minutes, stop_time=24hours)Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 1 day
├── Stop iteration : Inf
├── Wall time limit: Inf
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entriesAdd an output writer...
fnm_fridge = "soda_in_the_fridge.jld2"
outputs = (; co₂_flux,
soda_pCO₂,
atmos_pCO₂,
model_open_in_the_fridge.tracers.T,
model_open_in_the_fridge.tracers.DIC,
)
simulation.output_writers[:jld2] = JLD2OutputWriter(
model_open_in_the_fridge, outputs;
filename = fnm_fridge,
schedule = TimeInterval(10minutes),
overwrite_existing = true,
)JLD2OutputWriter scheduled on TimeInterval(10 minutes):
├── filepath: ./soda_in_the_fridge.jld2
├── 5 outputs: (co₂_flux, soda_pCO₂, atmos_pCO₂, T, DIC)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0.0 B... and don't forget to add a callback to compute the CO₂ flux
add_callback!(simulation, compute_co₂_flux!)Run the simulation
run!(simulation)[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (1.088 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.743 seconds).
[ Info: Simulation is stopping after running for 5.385 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.A simulation of a soda outgassing CO₂ on the counter
We simulate soda Warming up on the counter using a forcing function linear increase from 4°C to 30°C over 12 hours then stops warming
temperature_increase(z, t, p) = ifelse(t <= 12hours, p.∂T∂t * p.Δt , 0.0)
warming = Forcing(
temperature_increase,
parameters=(∂T∂t=1e-6, Δt=10minutes),
)ContinuousForcing{@NamedTuple{∂T∂t::Float64, Δt::Float64}}
├── func: temperature_increase (generic function with 1 method)
├── parameters: (∂T∂t = 1.0e-6, Δt = 600.0)
└── field dependencies: ()Build the second model for the soda on the counter
model_open_on_the_counter = NonhydrostaticModel(;
grid,
velocities = nothing,
buoyancy = nothing,
closure = nothing,
forcing = (; T=warming),
tracers = (:T, :DIC),
boundary_conditions = (; DIC=dic_bcs),
)
set!(model_open_on_the_counter, T = Tᵢ, DIC = DICᵢ)
simulation = Simulation(model_open_on_the_counter, Δt=10minutes, stop_time=24hours)Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 1 day
├── Stop iteration : Inf
├── Wall time limit: Inf
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entriesAdd an output writer...
fnm_counter = "soda_on_the_counter.jld2"
outputs = (; co₂_flux,
soda_pCO₂,
atmos_pCO₂,
model_open_on_the_counter.tracers.T,
model_open_on_the_counter.tracers.DIC,
)
simulation.output_writers[:jld2] = JLD2OutputWriter(
model_open_on_the_counter, outputs;
filename = fnm_counter,
schedule = TimeInterval(10minutes),
overwrite_existing = true)JLD2OutputWriter scheduled on TimeInterval(10 minutes):
├── filepath: ./soda_on_the_counter.jld2
├── 5 outputs: (co₂_flux, soda_pCO₂, atmos_pCO₂, T, DIC)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0.0 B..and don't forget to add a callback to compute the CO₂ flux
add_callback!(simulation, compute_co₂_flux!)Run the simulation
run!(simulation)[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (995.146 ms)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (2.019 seconds).
[ Info: Simulation is stopping after running for 3.570 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.Visualization
All that's left is to visualize the results.
fridge_soda_pCO₂ = FieldTimeSeries(fnm_fridge, "soda_pCO₂")
fridge_atmo_co₂ = FieldTimeSeries(fnm_fridge, "atmos_pCO₂")
fridge_temp = FieldTimeSeries(fnm_fridge, "T")
counter_soda_pCO₂ = FieldTimeSeries(fnm_counter, "soda_pCO₂")
counter_atmo_co₂ = FieldTimeSeries(fnm_counter, "atmos_pCO₂")
counter_temp = FieldTimeSeries(fnm_counter, "T")
t = fridge_soda_pCO₂.times
nt = length(t)
fig = Figure(size=(1200, 900))
ax = Axis(fig[1,1], xlabel="Time", ylabel="CO₂ conc [mmol m⁻³]")
lines!(t/(3600),
interior(fridge_soda_pCO₂, 1, 1, 1, :)*1e3;
linestyle = :dash,
label = "fridge soda",
)
lines!(t/(3600),
interior(counter_soda_pCO₂, 1, 1, 1, :)*1e3;
linestyle = :solid,
label = "counter soda",
)
lines!(t/(3600),
interior(fridge_atmo_co₂, 1, 1, 1, :)*1e3;
linestyle = :dash,
label = "fridge saturated")
lines!(t/(3600),
interior(counter_atmo_co₂, 1, 1, 1, :)*1e3;
linestyle = :solid,
label = "counter saturated")
axislegend()
ax = Axis(fig[2,1], xlabel="Time", ylabel="Temp (°C)")
lines!(t/(3600),
interior(fridge_temp, 1, 1, 1, :),
linestyle = :dash,
label = "fridge soda temperature",
)
lines!(t/(3600),
interior(counter_temp, 1, 1, 1, :),
linestyle = :solid,
label = "counter soda temperature",
)
axislegend()Makie.Legend()The cool soda's CO₂ concentration approaches equilibrium with the atmosphere (the saturated CO₂ concentration) quickly. The warming soda continues to outgas, since the solubility of CO₂ decreases with temperature. It'll taste flatter because of the lower CO₂ concentration.
#save("soda_outgassing_0d.png", fig)
fig
This page was generated using Literate.jl.