Merge branch 'main' into jsw/fix-fields-in-sinking-setup

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "bd1de93833dafc8fa6de30ae101c862f7fa5877e"
+project_hash = "37e11a3bdd396c973917441db34ded05b97faa85"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]

Project.toml

@@ -1,7 +1,7 @@
 name = "OceanBioME"
 uuid = "a49af516-9db8-4be4-be45-1dad61c5a376"
 authors = ["Jago Strong-Wright <[email protected]> and contributors"]
-version = "0.11.0"
+version = "0.11.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/Project.toml

@@ -19,6 +19,8 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 CairoMakie = "0.11"
 Documenter = "1"
 DocumenterCitations = "1"
-JLD2 = "0.4.38"
+EnsembleKalmanProcesses = "1"
+JLD2 = "0.4"
 Literate = "≥2.9.0"
-Oceananigans = "0.91.3"
+Oceananigans = "0.91"
+julia = "1.10"

src/BoxModel/timesteppers.jl

@@ -1,5 +1,6 @@
 using Oceananigans.Architectures: device
 using Oceananigans.Biogeochemistry: update_tendencies!, biogeochemical_auxiliary_fields
+using Oceananigans.Grids: nodes, Center
 using Oceananigans.TimeSteppers: rk3_substep_field!, store_field_tendencies!, RungeKutta3TimeStepper, QuasiAdamsBashforth2TimeStepper
 using Oceananigans.Utils: work_layout, launch!
 
@@ -39,7 +40,7 @@ function compute_tendencies!(model::BoxModel, callbacks)
     for tracer in required_biogeochemical_tracers(model.biogeochemistry)
         forcing = @inbounds model.forcing[tracer]
 
-        @inbounds Gⁿ[tracer][1, 1, 1] = tracer_tendency(Val(tracer), model.biogeochemistry, forcing, model.clock.time, model.field_values)
+        @inbounds Gⁿ[tracer][1, 1, 1] = tracer_tendency(Val(tracer), model.biogeochemistry, forcing, model.clock.time, model.field_values, model.grid)
     end
 
     for callback in callbacks
@@ -51,10 +52,14 @@ function compute_tendencies!(model::BoxModel, callbacks)
     return nothing
 end
 
-@inline tracer_tendency(val_name, biogeochemistry, forcing, time, model_fields) =
-    biogeochemistry(val_name, 0, 0, 0, time, model_fields...) + forcing(time, model_fields...)
+@inline boxmodel_xyz(nodes, grid) = map(n->boxmodel_coordinate(n, grid), nodes)
+@inline boxmodel_coordinate(::Nothing, grid) = zero(grid)
+@inline boxmodel_coordinate(nodes, grid) = @inbounds nodes[1]
 
-@inline tracer_tendency(::Val{:T}, biogeochemistry, forcing, time, model_fields) = 0
+@inline tracer_tendency(val_name, biogeochemistry, forcing, time, model_fields, grid) =
+    biogeochemistry(val_name, boxmodel_xyz(nodes(grid, Center(), Center(), Center()), grid)..., time, model_fields...) + forcing(time, model_fields...)
+
+@inline tracer_tendency(::Val{:T}, biogeochemistry, forcing, time, model_fields, grid) = 0
 
 function rk3_substep!(model::BoxModel, Δt, γⁿ, ζⁿ)
     model_fields = prognostic_fields(model)

src/Models/GasExchange/GasExchange.jl

@@ -40,6 +40,7 @@ include("gas_exchange.jl")
 include("carbon_dioxide_concentration.jl")
 include("schmidt_number.jl")
 include("gas_transfer_velocity.jl")
+include("gas_solubility.jl")
 
 using .ScaledGasTransferVelocity
 
@@ -130,7 +131,7 @@ specified by the the `OxygenConcentration` in the base model, and `air_concentra
 """
 OxygenGasExchangeBoundaryCondition(; transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = OxygenPolynomialSchmidtNumber()),
                                      water_concentration = OxygenConcentration(),
-                                     air_concentration = 9352.7, # mmolO₂/m³
+                                     air_concentration = PartiallySolubleGas(9352.7, OxygenSolubility()), # mmolO₂/m³
                                      wind_speed = 2,
                                      kwargs...) = 
     GasExchangeBoundaryCondition(; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...)

src/Models/GasExchange/gas_exchange.jl

@@ -30,7 +30,7 @@ end
 
     k = g.transfer_velocity(u₁₀, T)
 
-    air_concentration = surface_value(g.air_concentration, i, j, grid, clock)
+    air_concentration = surface_value(g.air_concentration, i, j, grid, clock, model_fields)
 
     water_concentration = surface_value(g.water_concentration, i, j, grid, clock, model_fields)
 

src/Models/GasExchange/gas_solubility.jl

@@ -0,0 +1,42 @@
+"""
+    PartiallySolubleGas(; air_concentration, solubility)
+
+Parameterises the available concentration of a gas dissolving in water in the form ``\\alpha C_a``
+where ``\alpha`` is the Ostwald solubility coeffieient and ``C_a`` is the concentration in the air.
+"""
+struct PartiallySolubleGas{AC, S}
+    air_concentration :: AC
+           solubility :: S
+end
+
+
+@inline surface_value(gs::PartiallySolubleGas, i, j, grid, clock, model_fields) = 
+    surface_value(gs.air_concentration, i, j, grid, clock) * surface_value(gs.solubility, i, j, grid, clock, model_fields)
+
+"""
+    Wanninkhof92Solubility
+
+Parameterises the Ostwald solubility coefficient as given in Wanninkhof, 1992.
+"""
+struct Wanninkhof92Solubility{FT}
+    A1 :: FT
+    A2 :: FT
+    A3 :: FT
+    B1 :: FT
+    B2 :: FT
+    B3 :: FT
+end
+
+function surface_value(sol::Wanninkhof92Solubility, i, j, grid, clock, model_fields)
+    Tk = @inbounds model_fields.T[i, j, grid.Nz] + 273.15
+    S = @inbounds model_fields.S[i, j, grid.Nz]
+
+    Tk_100 = Tk / 100
+
+    β = exp(sol.A1 + sol.A2  /Tk_100 + sol.A3 * log(Tk_100) + S * (sol.B1 + sol.B2 * Tk_100 + sol.B2 * Tk_100^2))
+
+    return β / Tk
+end
+
+OxygenSolubility(; A1 = -58.3877, A2 = 85.8079, A3 = 23.8439, B1 = -0.034892, B2 = 0.015578, B3 = -0.0019387) =
+    Wanninkhof92Solubility(A1, A2, A3, B1, B2, B3)

src/Models/GasExchange/surface_values.jl

@@ -13,7 +13,7 @@ end
     x = xnode(i, j, grid.Nz, grid, Center(), Center(), Center())
     y = ynode(i, j, grid.Nz, grid, Center(), Center(), Center())
 
-    return f.func(x, y, t, args...)
+    return f.func(x, y, t)
 end
 
 struct DiscreteSurfaceFuncton{F}

test/test_boxmodel.jl

@@ -1,14 +1,18 @@
 #include("dependencies_for_runtests.jl")
 
+using OceanBioME.BoxModels: boxmodel_xyz
+
 using Oceananigans: AbstractGrid, AbstractModel, AbstractOutputWriter
 using Oceananigans.Fields: FunctionField
+using Oceananigans.Grids: nodes
 
 PAR(t) = (60 * (1 - cos((t + 15days) * 2π / 365days)) * (1 / (1 + 0.2 * exp(-((mod(t, 365days) - 200days) / 50days)^2))) + 2) * exp(-2)
 
 defaults = (NO₃ = 10.0, NH₄ = 0.1, P = 0.1, Z = 0.01)
 
-function simple_box_model()
-    grid = BoxModelGrid()
+function simple_box_model(; x = nothing, y = nothing, z = nothing)
+    grid = BoxModelGrid(; x, y, z)
+
     clock = Clock(time = zero(grid))
 
     light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation(FunctionField{Center, Center, Center}(PAR, grid; clock))
@@ -21,11 +25,14 @@ function simple_box_model()
 end
 
 @testset "Construct `BoxModel`" begin
-    model = simple_box_model()
+    for z in (nothing, -100)
+        model = simple_box_model(; z)
 
-    @test grid isa AbstractGrid
-    @test model isa AbstractModel
-    @test set!(model; defaults...) isa Nothing
+        @test model.grid isa AbstractGrid
+        @test model isa AbstractModel
+        @test set!(model; defaults...) isa Nothing
+        @test boxmodel_xyz(nodes(model.grid, Center(), Center(), Center()), model.grid) == (0, 0, ifelse(isnothing(z), 0, z))
+    end
 end
 
 @testset "Timestep `BoxModel`" begin
@@ -38,8 +45,9 @@ end
     end
 
     # values are being updated
-    for (name, field) in enumerate(model.fields)
+    for (name, field) in pairs(model.fields)
         default = name in keys(defaults) ? defaults[name] : 0
+
         @test field[1, 1, 1] != default
     end
 end

test/test_gasexchange_carbon_chem.jl

@@ -132,6 +132,9 @@ end
     S = ConstantField(35)
     DIC = ConstantField(2136.242890518708)
     Alk = ConstantField(2500)
+
     # value from Dickson et. al, 2007
     @test ≈(surface_value(CO₂_exchange.water_concentration, 1, 1, BoxModelGrid(), Clock(; time = 0), (; T, S, DIC, Alk)), 350, atol = 0.1)
+
+    @test ≈(surface_value(O₂_exchange.air_concentration, 1, 1, BoxModelGrid(), Clock(; time = 0), (; T, S)), 200, atol = 50) # ball park correct
 end