Sediment model fixes + simpler model

src/Boundaries/Sediments/Sediments.jl

@@ -1,6 +1,6 @@
 module Sediments
 
-export SimpleMultiG
+export SimpleMultiG, InstantRemineralisation
 
 using KernelAbstractions
 using OceanBioME: ContinuousFormBiogeochemistry
@@ -26,13 +26,14 @@ sediment_fields(::AbstractSediment) = ()
 function update_tendencies!(bgc, sediment::FlatSediment, model)
     arch = model.grid.architecture
 
-    for (i, tracer) in enumerate(sediment_tracers(sediment))    
-        launch!(arch, model.grid, :xy, store_flat_tendencies!, sediment.tendencies.Gⁿ[i], sediment.tendencies.G⁻[i])
+    for (i, tracer) in enumerate(sediment_tracers(sediment))
+        launch!(arch, model.grid, :xy, store_flat_tendencies!, sediment.tendencies.G⁻[i], sediment.tendencies.Gⁿ[i])
     end
 
     launch!(arch, model.grid, :xy,
             _calculate_tendencies!,
-            bgc.sediment_model, bgc, model.grid, model.tracers, model.timestepper)
+            bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper)
+
     return nothing
 end
 
@@ -41,7 +42,13 @@ end
     @inbounds G⁻[i, j, 1] = G⁰[i, j, 1]
 end
 
+
+@inline nitrogen_flux(grid, adveciton, bgc, tracers, i, j) = 0
+@inline carbon_flux(grid, adveciton, bgc, tracers, i, j) = 0
+@inline remineralizaiton_reciever(bgc, tendencies) = nothing
+
 include("coupled_timesteppers.jl")
 include("simple_multi_G.jl")
+include("instant_remineralization.jl")
 
 end # module

src/Boundaries/Sediments/coupled_timesteppers.jl

@@ -1,13 +1,15 @@
-using Oceananigans: NonhydrostaticModel, prognostic_fields
+using Oceananigans: NonhydrostaticModel, prognostic_fields, HydrostaticFreeSurfaceModel
 using OceanBioME: ContinuousFormBiogeochemistry
 using OceanBioME.Boundaries.Sediments: AbstractSediment
 using Oceananigans.TimeSteppers: ab2_step_field!, rk3_substep_field!, stage_Δt
 using Oceananigans.Utils: work_layout, launch!
 using Oceananigans.TurbulenceClosures: implicit_step!
+using Oceananigans.Models.HydrostaticFreeSurfaceModels: local_ab2_step!, ab2_step_free_surface!
+using Oceananigans.Architectures: AbstractArchitecture
 
 import Oceananigans.TimeSteppers: ab2_step!, rk3_substep!
 
-@inline function ab2_step!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ)
+@inline function ab2_step!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ) 
     workgroup, worksize = work_layout(model.grid, :xyz)
     arch = model.architecture
     step_field_kernel! = ab2_step_field!(device(arch), workgroup, worksize)
@@ -44,6 +46,25 @@ import Oceananigans.TimeSteppers: ab2_step!, rk3_substep!
     return nothing
 end
 
+@inline function ab2_step!(model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ)
+    # Step locally velocity and tracers
+    @apply_regionally local_ab2_step!(model, Δt, χ)
+
+    # blocking step for implicit free surface, non blocking for explicit
+    ab2_step_free_surface!(model.free_surface, model, Δt, χ)
+
+    sediment = model.biogeochemistry.sediment_model
+    arch = model.architecture
+
+    for (i, field) in enumerate(sediment_fields(sediment))
+        launch!(arch, model.grid, :xy, ab2_step_flat_field!, 
+                field, Δt, χ, 
+                sediment.tendencies.Gⁿ[i], 
+                sediment.tendencies.G⁻[i])
+
+    end
+end
+
 @kernel function ab2_step_flat_field!(u, Δt, χ, Gⁿ, G⁻)
     i, j = @index(Global, NTuple)
 
@@ -83,15 +104,14 @@ function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:A
         launch!(arch, model.grid, :xy, rk3_step_flat_field!, 
                 field, Δt, γⁿ, ζⁿ,
                 sediment.tendencies.Gⁿ[i], 
-                sediment.tendencies.G⁻[i])
+                sediment.tendencies.G⁻[i]) 
     end
 
     return nothing
 end
 
 @kernel function rk3_step_flat_field!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
     i, j = @index(Global, NTuple)
-
     @inbounds U[i, j, 1] += Δt * (γⁿ * Gⁿ[i, j, 1] + ζⁿ * G⁻[i, j, 1])
 end
 

src/Boundaries/Sediments/instant_remineralization.jl

@@ -0,0 +1,97 @@
+"""
+    struct InstantRemineralisation
+
+Hold the parameters and fields the simplest benthic boundary layer where
+organic carbon is assumed to remineralise instantly with some portion 
+becoming N, and a fraction being perminantly burried.
+
+Burial efficiency by 10.1029/2006GB002907, 2007
+"""
+struct InstantRemineralisation{FT, F, TE} <: FlatSediment
+          burial_efficiency_constant1 :: FT
+          burial_efficiency_constant2 :: FT
+    burial_efficiency_half_saturaiton :: FT
+
+                               fields :: F
+                           tendencies :: TE
+
+    function InstantRemineralisation(burial_efficiency_constant1::FT,
+                                     burial_efficiency_constant2::FT,
+                                     burial_efficiency_half_saturaiton::FT,
+                                     fields::F,
+                                     tendencies::TE) where {FT, F, TE}
+
+        return new{FT, F, TE}(burial_efficiency_constant1,
+                              burial_efficiency_constant2,
+                              burial_efficiency_half_saturaiton,
+                              fields,
+                              tendencies)
+    end
+end
+
+"""
+    InstantRemineralisation(grid; 
+        burial_efficiency_constant1::FT = 0.013,
+        burial_efficiency_constant2::FT = 0.53,
+        burial_efficiency_half_saturaiton::FT = 7) where FT
+
+Returns a single-layer instant remineralisaiton model for NPZD bgc models.
+
+Example
+=======
+
+```@example
+using OceanBioME, Oceananigans, OceanBioME.Sediments
+
+grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))
+
+sediment_model = InstantRemineralisation(grid)
+```
+"""
+function InstantRemineralisation(; grid,
+            burial_efficiency_constant1 = 0.013,
+            burial_efficiency_constant2 = 0.53,
+            burial_efficiency_half_saturaiton = 7.0)
+
+    @warn "Sediment models are an experimental feature and have not yet been validated"
+
+    tracer_names = (:N_storage, )
+
+    # add field slicing back ( indices = (:, :, 1)) when output writer can cope
+    fields = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names))
+    tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)),
+                  G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)))
+
+    return InstantRemineralisation(burial_efficiency_constant1,
+                                   burial_efficiency_constant2,
+                                   burial_efficiency_half_saturaiton,
+                                   fields,
+                                   tendencies)
+end
+
+adapt_structure(to, sediment::InstantRemineralisation) = 
+    SimpleMultiG(sediment.burial_efficiency_constant1,
+                 sediment.burial_efficiency_constant2,
+                 sediment.burial_efficiency_half_saturaiton,
+                 adapt(to, sediment.fields),
+                 adapt(to, sediment.tendencies))
+
+sediment_tracers(::InstantRemineralisation) = (:N_storage, )
+sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, )
+
+@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, timestepper)
+    i, j = @index(Global, NTuple)
+
+    @inbounds begin                        
+        Δz = zspacing(i, j, 1, grid, Center(), Center(), Center())
+
+        flux = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz
+
+        burial_efficiency = sediment.burial_efficiency_constant1 + sediment.burial_efficiency_constant2 * ((flux / 6.56) / (7 + flux / 6.56)) ^ 2
+
+        # sediment evolution
+        sediment.tendencies.Gⁿ.N_storage[i, j, 1] = burial_efficiency * flux
+
+        remineralizaiton_reciever(bgc, timestepper.Gⁿ)[i, j, 1] += flux * (1 - burial_efficiency) / Δz
+    end
+end

src/Boundaries/Sediments/simple_multi_G.jl

@@ -64,40 +64,41 @@ grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))
 sediment_model = SimpleMultiG(grid)
 ```
 """
-function SimpleMultiG(grid; 
-                      fast_decay_rate::FT = 2/day,
-                      slow_decay_rate::FT = 0.2/day,
-                      fast_redfield::FT = 0.1509,
-                      slow_redfield::FT = 0.13,
-                      fast_fraction::FT = 0.74,
-                      slow_fraction::FT = 0.26,
-                      refactory_fraction::FT = 0.1,
-                      nitrate_oxidation_params::P1 = (A = - 1.9785, 
-                                                      B = 0.2261, 
-                                                      C = -0.0615, 
-                                                      D = -0.0289, 
-                                                      E = - 0.36109, 
-                                                      F = - 0.0232),
-                      denitrifcaiton_params::P2 = (A = - 3.0790, 
-                                                   B = 1.7509, 
-                                                   C = 0.0593, 
-                                                   D = - 0.1923, 
-                                                   E = 0.0604, 
-                                                   F = 0.0662),
-                      anoxic_params::P3 = (A = - 3.9476, 
-                                           B = 2.6269, 
-                                           C = - 0.2426, 
-                                           D = -1.3349, 
-                                           E = 0.1826, 
-                                           F = - 0.0143),
-                      depth = abs(znodes(grid, Face())[1]),
-                      solid_dep_params::P4 = (A = 0.233, 
-                                              B = 0.336, 
-                                              C = 982, 
-                                              D = - 1.548, 
-                                              depth = depth)) where {FT, P1, P2, P3, P4}
-
+function SimpleMultiG(; grid,
+                       fast_decay_rate::FT = 2/day,
+                       slow_decay_rate::FT = 0.2/day,
+                       fast_redfield::FT = 0.1509,
+                       slow_redfield::FT = 0.13,
+                       fast_fraction::FT = 0.74,
+                       slow_fraction::FT = 0.26,
+                       refactory_fraction::FT = 0.1,
+                       nitrate_oxidation_params::P1 = (A = - 1.9785, 
+                                                       B = 0.2261, 
+                                                       C = -0.0615, 
+                                                       D = -0.0289, 
+                                                       E = - 0.36109, 
+                                                       F = - 0.0232),
+                       denitrifcaiton_params::P2 = (A = - 3.0790, 
+                                                    B = 1.7509, 
+                                                    C = 0.0593, 
+                                                    D = - 0.1923, 
+                                                    E = 0.0604, 
+                                                    F = 0.0662),
+                       anoxic_params::P3 = (A  = - 3.9476, 
+                                            B = 2.6269, 
+                                            C = - 0.2426, 
+                                            D = -1.3349, 
+                                            E = 0.1826, 
+                                            F = - 0.0143),
+                       depth = abs(znodes(grid, Face())[1]),
+                       solid_dep_params::P4 = (A = 0.233, 
+                                               B = 0.336, 
+                                               C = 982, 
+                                               D = - 1.548, 
+                                               depth = depth)) where {FT, P1, P2, P3, P4}
+ 
     @warn "Sediment models are an experimental feature and have not yet been validated"
+    @info "This sediment model is only compatible with models providing NH₄, NO₃, O₂ and DIC"
 
     tracer_names = (:C_slow, :C_fast, :N_slow, :N_fast, :C_ref, :N_ref)
 
@@ -138,30 +139,37 @@ adapt_structure(to, sediment::SimpleMultiG) =
 sediment_tracers(::SimpleMultiG) = (:C_slow, :C_fast, :C_ref, :N_slow, :N_fast, :N_ref)
 sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, C_fast = model.fields.C_fast, N_slow = model.fields.N_slow, N_fast = model.fields.N_fast, C_ref = model.fields.C_ref, N_ref = model.fields.N_ref)
 
-@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, tracers, timestepper)
+@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, timestepper)
     i, j = @index(Global, NTuple)
 
+    Δz = zspacing(i, j, 1, grid, Center(), Center(), Center())
+
     @inbounds begin
-        carbon_deposition = -(biogeochemical_drift_velocity(bgc, Val(:sPOC)).w[i, j, 1] * tracers.sPOC[i, j, 1] +
-                              biogeochemical_drift_velocity(bgc, Val(:bPOC)).w[i, j, 1] * tracers.bPOC[i, j, 1])
+
+        carbon_deposition = carbon_flux(grid, advection, bgc, tracers, i, j) * Δz
 
-        nitrogen_deposition = -(biogeochemical_drift_velocity(bgc, Val(:sPON)).w[i, j, 1] * tracers.sPON[i, j, 1] +
-                                biogeochemical_drift_velocity(bgc, Val(:bPON)).w[i, j, 1] * tracers.bPON[i, j, 1])
+        nitrogen_deposition = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz
 
         # rates
-        Cᵐⁱⁿ = sediment.fields.C_slow[i, j, 1] * sediment.slow_decay_rate + sediment.fields.C_fast[i, j, 1] * sediment.fast_decay_rate
-        Nᵐⁱⁿ = sediment.fields.N_slow[i, j, 1] * sediment.slow_decay_rate + sediment.fields.N_fast[i, j, 1] * sediment.fast_decay_rate
+        C_min_slow = sediment.fields.C_slow[i, j, 1] * sediment.slow_decay_rate
+        C_min_fast = sediment.fields.C_fast[i, j, 1] * sediment.fast_decay_rate
+
+        N_min_slow = sediment.fields.N_slow[i, j, 1] * sediment.slow_decay_rate
+        N_min_fast = sediment.fields.N_fast[i, j, 1] * sediment.fast_decay_rate
+
+        Cᵐⁱⁿ = C_min_slow + C_min_fast
+        Nᵐⁱⁿ = N_min_slow + N_min_fast
 
         k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1])
 
         # sediment evolution
-        sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - sediment.slow_decay_rate * sediment.fields.C_slow[i, j, 1]
-        sediment.tendencies.Gⁿ.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - sediment.slow_decay_rate * sediment.fields.C_fast[i, j, 1]
-        sediment.tendencies.Gⁿ.C_ref[i, j, 1] = max(0.0, sediment.refactory_fraction * carbon_deposition)
+        sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow
+        sediment.tendencies.Gⁿ.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast
+        sediment.tendencies.Gⁿ.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition
 
-        sediment.tendencies.Gⁿ.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - sediment.slow_decay_rate * sediment.fields.N_slow[i, j, 1]
-        sediment.tendencies.Gⁿ.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - sediment.slow_decay_rate * sediment.fields.N_fast[i, j, 1]
-        sediment.tendencies.Gⁿ.N_ref[i, j, 1] = max(0.0, sediment.refactory_fraction * nitrogen_deposition)
+        sediment.tendencies.Gⁿ.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow
+        sediment.tendencies.Gⁿ.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast
+        sediment.tendencies.Gⁿ.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition
 
         # efflux/influx
         O₂  = tracers.O₂[i, j, 1]
@@ -191,15 +199,15 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, C_fast = m
 
         if isnan(pₐₙₒₓ)
             println("$(Cᵐⁱⁿ), $(k), $(O₂), $(NO₃)")
-            error("fucn")
+            error("Sediment anoxia has caused model failure")
         end
 
         pₛₒₗᵢ = sediment.solid_dep_params.A * (sediment.solid_dep_params.C * sediment.solid_dep_params.depth ^ sediment.solid_dep_params.D) ^ sediment.solid_dep_params.B
 
         Δz = grid.Δzᵃᵃᶜ[1]
 
-        timestepper.Gⁿ.NH₄[i, j, 1] += (Nᵐⁱⁿ * (1 - pₙᵢₜ)) / Δz
-        timestepper.Gⁿ.NO₃[i, j, 1] += (Nᵐⁱⁿ * pₙᵢₜ - Cᵐⁱⁿ * pᵈᵉⁿⁱᵗ * 4//5) / Δz
+        timestepper.Gⁿ.NH₄[i, j, 1] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz
+        timestepper.Gⁿ.NO₃[i, j, 1] += Nᵐⁱⁿ * pₙᵢₜ / Δz
         timestepper.Gⁿ.DIC[i, j, 1] += Cᵐⁱⁿ / Δz
         timestepper.Gⁿ.O₂[i, j, 1]  -= max(0, (1 - pᵈᵉⁿⁱᵗ - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ / Δz) # this seems dodge but this model doesn't cope with anoxia properly
     end