Merge pull request #211 from OceanBioME/jsw/fix-oxygen-gas-transfer

Fixes gas exchange for partially soluble gasses like oxygen

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/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_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