diff --git a/Manifest.toml b/Manifest.toml
index ce037c528..7c83fb905 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.5"
+julia_version = "1.10.6"
 manifest_format = "2.0"
-project_hash = "b7b68895e3965b7930805d42785ca0f68052cac9"
+project_hash = "6711c2855dab6dc90d90db711a2fef6fa66fe405"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -20,9 +20,9 @@ version = "1.5.0"
 
 [[deps.Accessors]]
 deps = ["CompositionsBase", "ConstructionBase", "InverseFunctions", "LinearAlgebra", "MacroTools", "Markdown"]
-git-tree-sha1 = "b392ede862e506d451fc1616e79aa6f4c673dab8"
+git-tree-sha1 = "96bed9b1b57cf750cca50c311a197e306816a1cc"
 uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.38"
+version = "0.1.39"
 
     [deps.Accessors.extensions]
     AccessorsAxisKeysExt = "AxisKeys"
@@ -45,9 +45,9 @@ version = "0.1.38"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "d80af0733c99ea80575f612813fa6aa71022d33a"
+git-tree-sha1 = "50c3c56a52972d78e8be9fd135bfb91c9574c140"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.1.0"
+version = "4.1.1"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -116,9 +116,9 @@ version = "5.5.2"
 
 [[deps.CUDA_Driver_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ccd1e54610c222fadfd4737dac66bff786f63656"
+git-tree-sha1 = "14996d716a2eaaeccfc8d7bc854dd87fde720ac1"
 uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc"
-version = "0.10.3+0"
+version = "0.10.4+0"
 
 [[deps.CUDA_Runtime_Discovery]]
 deps = ["Libdl"]
@@ -128,27 +128,33 @@ version = "0.3.5"
 
 [[deps.CUDA_Runtime_jll]]
 deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "e43727b237b2879a34391eeb81887699a26f8f2f"
+git-tree-sha1 = "17f1536c600133f7c4113bae0a2d98dbf27c7ebc"
 uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
-version = "0.15.3+0"
+version = "0.15.5+0"
 
 [[deps.ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d"
+git-tree-sha1 = "c7acce7a7e1078a20a285211dd73cd3941a871d6"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.11.5"
+version = "0.12.0"
+
+    [deps.ColorTypes.extensions]
+    StyledStringsExt = "StyledStrings"
+
+    [deps.ColorTypes.weakdeps]
+    StyledStrings = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
 
 [[deps.Colors]]
 deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
-git-tree-sha1 = "362a287c3aa50601b0bc359053d5c2468f0e7ce0"
+git-tree-sha1 = "64e15186f0aa277e174aa81798f7eb8598e0157e"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.12.11"
+version = "0.13.0"
 
 [[deps.CommonDataModel]]
 deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf", "Statistics"]
-git-tree-sha1 = "d6fb5bf939a2753c74984b11434ea25d6c397a58"
+git-tree-sha1 = "98d64d5b9e5263884276656a43c45424b3a645c2"
 uuid = "1fbeeb36-5f17-413c-809b-666fb144f157"
-version = "0.3.6"
+version = "0.3.7"
 
 [[deps.CommonSolve]]
 git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
@@ -233,9 +239,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[deps.DiskArrays]]
 deps = ["LRUCache", "OffsetArrays"]
-git-tree-sha1 = "ef25c513cad08d7ebbed158c91768ae32f308336"
+git-tree-sha1 = "e0e89a60637a62d13aa2107f0acd169b9b9b77e7"
 uuid = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
-version = "0.3.23"
+version = "0.4.6"
 
 [[deps.Distances]]
 deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
@@ -285,9 +291,15 @@ version = "3.3.10+1"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "62ca0547a14c57e98154423419d8a342dca75ca9"
+git-tree-sha1 = "2dd20384bf8c6d411b5c7370865b1e9b26cb2ea3"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.4"
+version = "1.16.6"
+
+    [deps.FileIO.extensions]
+    HTTPExt = "HTTP"
+
+    [deps.FileIO.weakdeps]
+    HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
@@ -302,11 +314,6 @@ version = "0.8.5"
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
-[[deps.GMP_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
-version = "6.2.1+6"
-
 [[deps.GPUArrays]]
 deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]
 git-tree-sha1 = "62ee71528cca49be797076a76bdc654a170a523e"
@@ -342,12 +349,6 @@ git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
 uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
 version = "1.3.1"
 
-[[deps.GnuTLS_jll]]
-deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Nettle_jll", "P11Kit_jll", "Zlib_jll"]
-git-tree-sha1 = "383db7d3f900f4c1f47a8a04115b053c095e48d3"
-uuid = "0951126a-58fd-58f1-b5b3-b08c7c4a876d"
-version = "3.8.4+0"
-
 [[deps.HDF5_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
 git-tree-sha1 = "38c8874692d48d5440d5752d6c74b0c6b0b60739"
@@ -356,9 +357,9 @@ version = "1.14.2+1"
 
 [[deps.Hwloc_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "dd3b49277ec2bb2c6b94eb1604d4d0616016f7a6"
+git-tree-sha1 = "50aedf345a709ab75872f80a2779568dc0bb461b"
 uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
-version = "2.11.2+0"
+version = "2.11.2+1"
 
 [[deps.IncompleteLU]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -400,9 +401,9 @@ weakdeps = ["Dates", "Test"]
     InverseFunctionsTestExt = "Test"
 
 [[deps.InvertedIndices]]
-git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
+git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe"
 uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.IterativeSolvers]]
 deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
@@ -417,9 +418,9 @@ version = "1.0.0"
 
 [[deps.JLD2]]
 deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "Requires", "TranscodingStreams"]
-git-tree-sha1 = "783c1be5213a09609b23237a0c9e5dfd258ae6f2"
+git-tree-sha1 = "f1a1c1037af2a4541ea186b26b0c0e7eeaad232b"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.5.7"
+version = "0.5.10"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
@@ -434,10 +435,10 @@ uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e"
 version = "0.2.1+0"
 
 [[deps.KernelAbstractions]]
-deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
-git-tree-sha1 = "e73a077abc7fe798fe940deabe30ef6c66bdde52"
+deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs"]
+git-tree-sha1 = "b9a838cd3028785ac23822cded5126b3da394d1a"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.29"
+version = "0.9.31"
 
     [deps.KernelAbstractions.extensions]
     EnzymeExt = "EnzymeCore"
@@ -595,9 +596,9 @@ version = "2.28.2+1"
 
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f12a29c4400ba812841c6ace3f4efbb6dbb3ba01"
+git-tree-sha1 = "bc95bf4149bf535c09602e3acdf950d9b4376227"
 uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf"
-version = "10.1.4+2"
+version = "10.1.4+3"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -614,15 +615,15 @@ version = "2023.1.10"
 
 [[deps.NCDatasets]]
 deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "DiskArrays", "NetCDF_jll", "NetworkOptions", "Printf"]
-git-tree-sha1 = "77df6d3708ec0eb3441551e1f20f7503b37c2393"
+git-tree-sha1 = "2c9dc92001ac06d432f363f37ff5552954d9947c"
 uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
-version = "0.14.5"
+version = "0.14.6"
 
 [[deps.NVTX]]
 deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"]
-git-tree-sha1 = "53046f0483375e3ed78e49190f1154fa0a4083a1"
+git-tree-sha1 = "6a6f8bfaa91bb2e40ff562ab9f30dc827741daef"
 uuid = "5da4648a-3479-48b8-97b9-01cb529c0a1f"
-version = "0.3.4"
+version = "0.3.5"
 
 [[deps.NVTX_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -636,21 +637,15 @@ git-tree-sha1 = "a8af1798e4eb9ff768ce7fdefc0e957097793f15"
 uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3"
 version = "400.902.209+0"
 
-[[deps.Nettle_jll]]
-deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
-uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
-version = "3.7.2+0"
-
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
 
 [[deps.Oceananigans]]
 deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
-git-tree-sha1 = "cddb567fbab1c40b22914e24ad2fe917ff00d8c9"
+git-tree-sha1 = "6ea1bd01adf2ae95d96a049711b9db4e24959767"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.93.1"
+version = "0.95.5"
 
     [deps.Oceananigans.extensions]
     OceananigansEnzymeExt = "Enzyme"
@@ -662,9 +657,9 @@ version = "0.93.1"
     MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
+git-tree-sha1 = "5e1897147d1ff8d98883cda2be2187dcf57d8f0c"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.14.1"
+version = "1.15.0"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -677,9 +672,9 @@ version = "0.3.23+4"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"]
-git-tree-sha1 = "bfce6d523861a6c562721b262c0d1aaeead2647f"
+git-tree-sha1 = "2dace87e14256edb1dd0724ab7ba831c779b96bd"
 uuid = "fe0851c0-eecd-5654-98d4-656369965a5c"
-version = "5.0.5+0"
+version = "5.0.6+0"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -688,15 +683,9 @@ uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
 version = "3.0.15+1"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
+git-tree-sha1 = "12f1439c4f986bb868acda6ea33ebc78e19b95ad"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.3"
-
-[[deps.P11Kit_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "2cd396108e178f3ae8dedbd8e938a18726ab2fbf"
-uuid = "c2071276-7c44-58a7-b746-946036e04d0a"
-version = "0.24.1+0"
+version = "1.7.0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -788,9 +777,9 @@ version = "1.3.0"
 
 [[deps.Roots]]
 deps = ["Accessors", "CommonSolve", "Printf"]
-git-tree-sha1 = "3a7c7e5c3f015415637f5debdf8a674aa2c979c4"
+git-tree-sha1 = "8e3694d669323cdfb560e344dc872b984de23b71"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "2.2.1"
+version = "2.2.2"
 
     [deps.Roots.extensions]
     RootsChainRulesCoreExt = "ChainRulesCore"
@@ -833,9 +822,9 @@ version = "0.3.5"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "305becf8af67eae1dbc912ee9097f00aeeabb8d5"
+git-tree-sha1 = "712fb0231ee6f9120e005ccd56297abbc053e7e0"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.6"
+version = "1.4.8"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -856,9 +845,9 @@ version = "1.10.0"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "777657803913ffc7e8cc20f0fd04b634f871af8f"
+git-tree-sha1 = "7c01731da8ab6d3094c4d44c9057b00932459255"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.8"
+version = "1.9.9"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -950,9 +939,9 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[deps.TimerOutputs]]
 deps = ["ExprTools", "Printf"]
-git-tree-sha1 = "3a6f063d690135f5c1ba351412c82bae4d1402bf"
+git-tree-sha1 = "d7298ebdfa1654583468a487e8e83fae9d72dac3"
 uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
-version = "0.5.25"
+version = "0.5.26"
 
 [[deps.TranscodingStreams]]
 git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
@@ -971,17 +960,11 @@ git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278"
 uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
 version = "0.2.1"
 
-[[deps.UnsafeAtomicsLLVM]]
-deps = ["LLVM", "UnsafeAtomics"]
-git-tree-sha1 = "2d17fabcd17e67d7625ce9c531fb9f40b7c42ce4"
-uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249"
-version = "0.2.1"
-
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "6a451c6f33a176150f315726eba8b92fbfdb9ae7"
+git-tree-sha1 = "a2fccc6559132927d4c5dc183e3e01048c6dcbd6"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.13.4+0"
+version = "2.13.5+0"
 
 [[deps.XZ_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1018,10 +1001,10 @@ uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
 version = "5.11.0+0"
 
 [[deps.libzip_jll]]
-deps = ["Artifacts", "Bzip2_jll", "GnuTLS_jll", "JLLWrappers", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "668ac0297e6bd8f4d53dfdcd3ace71f2e00f4a35"
+deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "OpenSSL_jll", "XZ_jll", "Zlib_jll", "Zstd_jll"]
+git-tree-sha1 = "e797fa066eba69f4c0585ffbd81bc780b5118ce2"
 uuid = "337d8026-41b4-5cde-a456-74a10e5b31d1"
-version = "1.11.1+0"
+version = "1.11.2+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
diff --git a/Project.toml b/Project.toml
index ab8dbbfda..49b7ea19c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "OceanBioME"
 uuid = "a49af516-9db8-4be4-be45-1dad61c5a376"
 authors = ["Jago Strong-Wright <js2430@damtp.cam.ac.uk> and contributors"]
-version = "0.13.3"
+version = "0.13.4"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -24,7 +24,7 @@ EnsembleKalmanProcesses = "1"
 GibbsSeaWater = "0.1"
 JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
-Oceananigans = "^0.92.1, ^0.93.0"
+Oceananigans = "^0.95"
 OffsetArrays = "1.14"
 Roots = "2"
 SeawaterPolynomials = "0.3"
diff --git a/benchmark/benchmark_simulations.jl b/benchmark/benchmark_simulations.jl
index 1c45d4fed..961a7ed2e 100644
--- a/benchmark/benchmark_simulations.jl
+++ b/benchmark/benchmark_simulations.jl
@@ -1,28 +1,34 @@
-using OceanBioME, Oceananigans, BenchmarkTools
-using OceanBioME.Particles: infinitesimal_particle_field_coupling!
+using BenchmarkTools
+using Oceananigans
+using OceanBioME
+
 using Oceananigans.Units
+
+using Oceananigans: TendencyCallsite
+using OceanBioME.Particles: infinitesimal_particle_field_coupling!
+
 include("Benchmark.jl")
 
 function benchmark_LOBSTER(n)
-    grid = RectilinearGrid(size=(n, n, n), extent=(200, 200, 200)) 
-    PAR = Oceananigans.Fields.Field{Center, Center, Center}(grid)  
+    grid = RectilinearGrid(size=(n, n, n), extent=(200, 200, 200))
+    PAR = Oceananigans.Fields.Field{Center, Center, Center}(grid)
 
     @inline PAR⁰(t) = 10.0
 
-    bgc = Setup.Oceananigans(:LOBSTER, grid, LOBSTER.defaults) 
+    bgc = Setup.Oceananigans(:LOBSTER, grid, LOBSTER.defaults)
 
     model = NonhydrostaticModel(
         advection = WENO(;grid),
         timestepper = :RungeKutta3,
         grid = grid,
         tracers = bgc.tracers,
-        closure = ScalarDiffusivity(ν=1e-4, κ=1e-4), 
+        closure = ScalarDiffusivity(ν=1e-4, κ=1e-4),
         forcing = bgc.forcing,
         boundary_conditions = bgc.boundary_conditions,
         auxiliary_fields = (; PAR)
     )
     set!(model, P=0.03, Z=0.03, NO₃=11, NH₄=0.05)
-    simulation = Simulation(model, Δt=10minutes, stop_iteration=1) 
+    simulation = Simulation(model, Δt=10minutes, stop_iteration=1)
 
     simulation.callbacks[:update_par] = Callback(Light.twoBands.update!, IterationInterval(1), merge(merge(LOBSTER.defaults, Light.twoBands.defaults), (surface_PAR=PAR⁰,)), TimeStepCallsite());
     simulation.callbacks[:neg] = Callback(scale_negative_tracers!; parameters=(conserved_group=(:NO₃, :NH₄, :P, :Z, :D, :DD, :DOM), warn=false))
@@ -36,8 +42,8 @@ function benchmark_LOBSTER(n)
 end
 
 function benchmark_sediment(n)
-    grid = RectilinearGrid(size=(n, n, n), extent=(200, 200, 200)) 
-    PAR = Oceananigans.Fields.Field{Center, Center, Center}(grid)  
+    grid = RectilinearGrid(size=(n, n, n), extent=(200, 200, 200))
+    PAR = Oceananigans.Fields.Field{Center, Center, Center}(grid)
 
     @inline PAR⁰(t) = 10.0
     @inline t_function(x, y, z, t) = 13.0
@@ -49,7 +55,7 @@ function benchmark_sediment(n)
 
     # Have to prespecify fields for sediment model
     NO₃, NH₄, P, Z, D, DD, Dᶜ, DDᶜ, DOM, DIC, ALK, OXY = CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid), CenterField(grid)
-    
+
     #sediment bcs
     sediment=Boundaries.Sediments.Soetaert.setupsediment(grid, (;D, DD); POM_w=(D=LOBSTER.D_sinking, DD=LOBSTER.DD_sinking))
 
@@ -60,13 +66,13 @@ function benchmark_sediment(n)
         timestepper = :RungeKutta3,
         grid = grid,
         tracers = bgc.tracers,
-        closure = ScalarDiffusivity(ν=1e-4, κ=1e-4), 
+        closure = ScalarDiffusivity(ν=1e-4, κ=1e-4),
         forcing =  merge(bgc.forcing, sediment.forcing),
         boundary_conditions = bgc.boundary_conditions,
         auxiliary_fields = merge((; PAR), sediment.auxiliary_fields)
     )
     set!(model, P=0.03, Z=0.03, NO₃=11, NH₄=0.05)
-    simulation = Simulation(model, Δt=10minutes, stop_iteration=1) 
+    simulation = Simulation(model, Δt=10minutes, stop_iteration=1)
 
     simulation.callbacks[:update_par] = Callback(Light.twoBands.update!, IterationInterval(1), merge(merge(LOBSTER.defaults, Light.twoBands.defaults), (surface_PAR=PAR⁰,)), TimeStepCallsite());
     simulation.callbacks[:neg] = Callback(scale_negative_tracers!; parameters=(conserved_group=(:NO₃, :NH₄, :P, :Z, :D, :DD, :DOM), warn=false))
@@ -80,8 +86,8 @@ function benchmark_sediment(n)
 end
 
 function benchmark_SLatissima(n)
-    grid = RectilinearGrid(size=(n, n, n), extent=(200, 200, 200)) 
-    PAR = Oceananigans.Fields.Field{Center, Center, Center}(grid)  
+    grid = RectilinearGrid(size=(n, n, n), extent=(200, 200, 200))
+    PAR = Oceananigans.Fields.Field{Center, Center, Center}(grid)
 
     @inline PAR⁰(t) = 10.0
     @inline t_function(x, y, z, t) = 13.0
@@ -90,29 +96,29 @@ function benchmark_SLatissima(n)
     n_kelp = floor(Int, n/4) # number of kelp fronds
     z₀ = [-101+100/n_kelp:100/n_kelp:-1;]*1.0 # depth of kelp fronds
 
-    kelp_particles = SLatissima.setup(n_kelp, 200*rand(n_kelp), 200*rand(n_kelp), z₀, 
+    kelp_particles = SLatissima.setup(n_kelp, 200*rand(n_kelp), 200*rand(n_kelp), z₀,
                                                         30.0, 0.01, 0.1, 57.5;
-                                                        scalefactor = 5.0, 
-                                                        T = t_function, S = s_function, urel = 0.2, 
+                                                        scalefactor = 5.0,
+                                                        T = t_function, S = s_function, urel = 0.2,
                                                         optional_tracers = (:NH₄, :DIC, :DD, :DDᶜ, :OXY, :DOM))
 
     dic_bc = Boundaries.airseasetup(:CO₂, forcings=(T=t_function, S=s_function))
     oxy_bc = Boundaries.airseasetup(:O₂, forcings=(T=t_function, S=s_function))
-    bgc = Setup.Oceananigans(:LOBSTER, grid, LOBSTER.defaults, optional_sets=(:carbonates, :oxygen), topboundaries=(DIC=dic_bc, OXY=oxy_bc)) 
+    bgc = Setup.Oceananigans(:LOBSTER, grid, LOBSTER.defaults, optional_sets=(:carbonates, :oxygen), topboundaries=(DIC=dic_bc, OXY=oxy_bc))
 
     model = NonhydrostaticModel(
         advection = WENO(;grid),
         timestepper = :RungeKutta3,
         grid = grid,
         tracers = bgc.tracers,
-        closure = ScalarDiffusivity(ν=1e-4, κ=1e-4), 
+        closure = ScalarDiffusivity(ν=1e-4, κ=1e-4),
         forcing = bgc.forcing,
         boundary_conditions = bgc.boundary_conditions,
         auxiliary_fields = (; PAR),
         particles = kelp_particles
     )
     set!(model, P=0.03, Z=0.03, NO₃=11, NH₄=0.05, DIC=2200.0, ALK=2400.0, OXY=240.0)
-    simulation = Simulation(model, Δt=10minutes, stop_iteration=1) 
+    simulation = Simulation(model, Δt=10minutes, stop_iteration=1)
 
     simulation.callbacks[:update_par] = Callback(Light.twoBands.update!, IterationInterval(1), merge(merge(LOBSTER.defaults, Light.twoBands.defaults), (surface_PAR=PAR⁰,)), TimeStepCallsite());
     simulation.callbacks[:neg] = Callback(scale_negative_tracers!; parameters=(conserved_group=(:NO₃, :NH₄, :P, :Z, :D, :DD, :DOM), warn=false))
diff --git a/docs/Project.toml b/docs/Project.toml
index 742e5d79e..23d52f318 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -20,7 +20,7 @@ CairoMakie = "0.11, 0.12"
 Documenter = "1"
 DocumenterCitations = "1"
 EnsembleKalmanProcesses = "1"
-JLD2 = "0.4"
+JLD2 = "0.4, 0.5"
 Literate = "≥2.9.0"
-Oceananigans = "^0.92.1"
+Oceananigans = "^0.95"
 julia = "1.10"
diff --git a/docs/src/model_components/sediments/simple_multi_g.md b/docs/src/model_components/sediments/simple_multi_g.md
index 5c1a26696..08ed2a56d 100644
--- a/docs/src/model_components/sediments/simple_multi_g.md
+++ b/docs/src/model_components/sediments/simple_multi_g.md
@@ -14,16 +14,18 @@ sediment_model = SimpleMultiG(; grid)
 # output
 ┌ Warning: Sediment models are an experimental feature and have not yet been validated.
 └ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:104
-[ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
+┌ Warning: Sediment models currently do not pass tests and are probably broken!
+└ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:105
+[ Info: The SimpleMultiG sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
 Single-layer multi-G sediment model (Float64)
 ```
 
 You may optionally specify the model parameters. This can then be passed in the setup of a BGC model:
 
 ```julia
-biogeochemistry = LOBSTER(; grid, 
-                            carbonates = true, oxygen = true, variable_redfield = true, 
-                            open_bottom = true, 
+biogeochemistry = LOBSTER(; grid,
+                            carbonates = true, oxygen = true, variable_redfield = true,
+                            open_bottom = true,
                             sediment_model)
 ```
 
@@ -87,7 +89,7 @@ The original model of [Soetaert2000](@citet) also includes denitrification terms
 | Symbol                  | Variable name              | Units |
 |-------------------------|----------------------------|-------|
 | ``\lambda_\text{fast}`` | `fast_decay_rate`          | 1 / s |
-| ``\lambda_\text{slow}`` | `slow_decay_rate`          | 1 / s | 
+| ``\lambda_\text{slow}`` | `slow_decay_rate`          | 1 / s |
 | ``f_\text{fast}``       | `fast_fraction`            | -     |
 | ``f_\text{slow}``       | `slow_fraction`            | -     |
 | ``f_\text{ref}``        | `refactory_fraction`       | -     |
@@ -95,7 +97,7 @@ The original model of [Soetaert2000](@citet) also includes denitrification terms
 | ``a_i``                 | `anoxic_param`             | -     |
 | ``s_i``                 | `solid_dep_params`         | -     |
 
-All parameters are given in [Parameters](@ref parameters). 
+All parameters are given in [Parameters](@ref parameters).
 
 ## Model conservations
 
@@ -103,4 +105,4 @@ Nitrogen and carbon is conserved between the model domain and sediment, any nitr
 
 ## Model compatibility
 
-This model is currently only compatible with the [LOBSTER](@ref LOBSTER) biogeochemical model.
\ No newline at end of file
+This model is currently only compatible with the [LOBSTER](@ref LOBSTER) biogeochemical model.
diff --git a/src/Models/AdvectedPopulations/PISCES/common.jl b/src/Models/AdvectedPopulations/PISCES/common.jl
index dc96a45fc..0369ef416 100644
--- a/src/Models/AdvectedPopulations/PISCES/common.jl
+++ b/src/Models/AdvectedPopulations/PISCES/common.jl
@@ -1,12 +1,7 @@
 using KernelAbstractions: @kernel, @index
 
-using Oceananigans.Fields: flatten_node
-using Oceananigans.Grids: znode, zspacing, φnode
-
-import Oceananigans.Fields: flatten_node
-
-# TODO: move this to Oceananigans
-@inline flatten_node(::Nothing, ::Nothing, z) = tuple(z)
+using Oceananigans.Grids: znode, φnode
+using Oceananigans.Operators: zspacing
 
 @inline shear(z, zₘₓₗ, background_shear, mixed_layer_shear) = ifelse(z <= zₘₓₗ, background_shear, mixed_layer_shear) # Given as 1 in Aumont paper
 
@@ -37,7 +32,7 @@ end
                                  maximum_speed = 200/day,
                                  maximum_depth = 500)
 
-Returns sinking speed for particles which sink at `minimum_speed` in the 
+Returns sinking speed for particles which sink at `minimum_speed` in the
 surface ocean (the deepest of the mixed and euphotic layers), and accelerate
 to `maximum_speed` below that depth and `maximum_depth`.
 """
@@ -60,7 +55,7 @@ end
 end
 
 # don't actually use this version but might be useful since we can do it
-@inline (p::DepthDependantSinkingSpeed)(z, zₘₓₗ, zₑᵤ) = 
+@inline (p::DepthDependantSinkingSpeed)(z, zₘₓₗ, zₑᵤ) =
     ifelse(z < 0, - p.minimum_speed + (p.maximum_speed - p.minimum_speed) * min(0, z - min(zₘₓₗ, zₑᵤ)) / 5000, 0)
 
 @inline function anoxia_factor(bgc, O₂)
diff --git a/src/Models/AdvectedPopulations/PISCES/compute_calcite_saturation.jl b/src/Models/AdvectedPopulations/PISCES/compute_calcite_saturation.jl
index 112e97b77..25876fcd9 100644
--- a/src/Models/AdvectedPopulations/PISCES/compute_calcite_saturation.jl
+++ b/src/Models/AdvectedPopulations/PISCES/compute_calcite_saturation.jl
@@ -1,6 +1,6 @@
 using Oceananigans.Architectures: architecture
 using Oceananigans.BoundaryConditions: fill_halo_regions!
-using Oceananigans.BuoyancyModels: g_Earth
+using Oceananigans.BuoyancyFormulations: g_Earth
 using Oceananigans.Models: fields
 using Oceananigans.Utils: launch!
 
diff --git a/src/Models/Sediments/Sediments.jl b/src/Models/Sediments/Sediments.jl
index 78ee735a2..df317f086 100644
--- a/src/Models/Sediments/Sediments.jl
+++ b/src/Models/Sediments/Sediments.jl
@@ -13,8 +13,7 @@ using Oceananigans.Advection: advective_tracer_flux_z, FluxFormAdvection
 using Oceananigans.Units: day
 using Oceananigans.Fields: ConstantField
 using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity
-using Oceananigans.Grids: zspacing
-using Oceananigans.Operators: volume
+using Oceananigans.Operators: volume, zspacing
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, immersed_cell
 
 import Adapt: adapt_structure, adapt
diff --git a/src/Models/Sediments/coupled_timesteppers.jl b/src/Models/Sediments/coupled_timesteppers.jl
index 2888d8f6c..c7619b125 100644
--- a/src/Models/Sediments/coupled_timesteppers.jl
+++ b/src/Models/Sediments/coupled_timesteppers.jl
@@ -1,10 +1,11 @@
-using Oceananigans: NonhydrostaticModel, prognostic_fields, HydrostaticFreeSurfaceModel
-using OceanBioME.Models.Sediments: AbstractSediment
+using Oceananigans: NonhydrostaticModel, HydrostaticFreeSurfaceModel, prognostic_fields, location
 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.Models.HydrostaticFreeSurfaceModels: step_free_surface!, local_ab2_step!
 using Oceananigans.Architectures: AbstractArchitecture
+using Oceananigans.Utils: @apply_regionally
+using OceanBioME.Models.Sediments: AbstractSediment
 
 import Oceananigans.TimeSteppers: ab2_step!, rk3_substep!
 
@@ -58,7 +59,7 @@ end
     @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, χ)
+    step_free_surface!(model.free_surface, model, Δt, χ)
 
     sediment = model.biogeochemistry.sediment
     arch = model.architecture
diff --git a/src/Models/Sediments/instant_remineralization.jl b/src/Models/Sediments/instant_remineralization.jl
index 7f5a45008..8bc4776c1 100644
--- a/src/Models/Sediments/instant_remineralization.jl
+++ b/src/Models/Sediments/instant_remineralization.jl
@@ -2,7 +2,7 @@
     struct InstantRemineralisation
 
 Hold the parameters and fields the simplest benthic boundary layer where
-organic carbon is assumed to remineralise instantly with some portion 
+organic carbon is assumed to remineralise instantly with some portion
 becoming N, and a fraction being permanently buried.
 
 Burial efficiency from [RemineralisationFraction](@citet).
@@ -54,7 +54,8 @@ function InstantRemineralisation(; grid,
             burial_efficiency_constant2 = 0.53,
             burial_efficiency_half_saturation = 7.0)
 
-    @warn "Sediment models are an experimental feature and have not yet been validated"
+    @warn "Sediment models are an experimental feature and have not yet been validated."
+    @warn "Sediment models currently do not pass tests and are probably broken!"
 
     tracer_names = (:N_storage, )
 
@@ -73,14 +74,14 @@ function InstantRemineralisation(; grid,
                                    bottom_indices)
 end
 
-adapt_structure(to, sediment::InstantRemineralisation) = 
+adapt_structure(to, sediment::InstantRemineralisation) =
     InstantRemineralisation(adapt(to, sediment.burial_efficiency_constant1),
                             adapt(to, sediment.burial_efficiency_constant2),
                             adapt(to, sediment.burial_efficiency_half_saturation),
                             adapt(to, sediment.fields),
                             nothing,
                             adapt(to, sediment.bottom_indices))
-                  
+
 sediment_tracers(::InstantRemineralisation) = (:N_storage, )
 sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, )
 
@@ -105,4 +106,4 @@ function _calculate_sediment_tendencies!(i, j, sediment::InstantRemineralisation
 end
 
 summary(::InstantRemineralisation{FT}) where {FT} = string("Single-layer instant remineralisaiton ($FT)")
-show(io::IO, model::InstantRemineralisation) = print(io, summary(model))
\ No newline at end of file
+show(io::IO, model::InstantRemineralisation) = print(io, summary(model))
diff --git a/src/Models/Sediments/simple_multi_G.jl b/src/Models/Sediments/simple_multi_G.jl
index 138049c27..0c559dd88 100644
--- a/src/Models/Sediments/simple_multi_G.jl
+++ b/src/Models/Sediments/simple_multi_G.jl
@@ -25,11 +25,11 @@ struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B} <: FlatSediment
                       fields :: F
                   tendencies :: TE
               bottom_indices :: B
-    
+
     SimpleMultiG(fast_decay_rate::FT, slow_decay_rate::FT,
                  fast_redfield::FT,   slow_redfield::FT,
                  fast_fraction::FT,   slow_fraction::FT, refactory_fraction::FT,
-                 nitrate_oxidation_params::P1, 
+                 nitrate_oxidation_params::P1,
                  denitrification_params::P2,
                  anoxic_params::P3,
                  solid_dep_params::P4,
@@ -38,7 +38,7 @@ struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B} <: FlatSediment
         new{FT, P1, P2, P3, P4, F, TE, B}(fast_decay_rate, slow_decay_rate,
                                           fast_redfield,   slow_redfield,
                                           fast_fraction,   slow_fraction, refactory_fraction,
-                                          nitrate_oxidation_params, 
+                                          nitrate_oxidation_params,
                                           denitrification_params,
                                           anoxic_params,
                                           solid_dep_params,
@@ -84,7 +84,9 @@ julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200));
 julia> sediment_model = SimpleMultiG(; grid)
 ┌ Warning: Sediment models are an experimental feature and have not yet been validated.
 └ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:104
-[ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
+┌ Warning: Sediment models currently do not pass tests and are probably broken!
+└ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:105
+[ Info: The SimpleMultiG sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
 Single-layer multi-G sediment model (Float64)
 ```
 """
@@ -102,7 +104,8 @@ function SimpleMultiG(; grid,
                         solid_dep_params = on_architecture(architecture(grid), [0.233, 0.336, 982.0, - 1.548]))
 
     @warn "Sediment models are an experimental feature and have not yet been validated."
-    @info "This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC."
+    @warn "Sediment models currently do not pass tests and are probably broken!"
+    @info "The SimpleMultiG sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC."
 
     tracer_names = (:C_slow, :C_fast, :N_slow, :N_fast, :C_ref, :N_ref)
 
@@ -125,7 +128,7 @@ function SimpleMultiG(; grid,
                         bottom_indices)
 end
 
-adapt_structure(to, sediment::SimpleMultiG) = 
+adapt_structure(to, sediment::SimpleMultiG) =
     SimpleMultiG(adapt(to, sediment.fast_decay_rate),
                  adapt(to, sediment.slow_decay_rate),
                  adapt(to, sediment.fast_redfield),
@@ -140,7 +143,7 @@ adapt_structure(to, sediment::SimpleMultiG) =
                  adapt(to, sediment.fields),
                  nothing,
                  adapt(to, sediment.bottom_indices))
-                  
+
 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,
@@ -162,7 +165,7 @@ function _calculate_sediment_tendencies!(i, j, sediment::SimpleMultiG, bgc, grid
 
     @inbounds begin
         carbon_deposition = carbon_flux(i, j, k, grid, advection, bgc, tracers) * Δz
-                        
+
         nitrogen_deposition = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz
 
         # rates
@@ -174,7 +177,7 @@ function _calculate_sediment_tendencies!(i, j, sediment::SimpleMultiG, bgc, grid
 
         Cᵐⁱⁿ = C_min_slow + C_min_fast
         Nᵐⁱⁿ = N_min_slow + N_min_fast
-        
+
         reactivity = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1])
 
         # sediment evolution
@@ -190,7 +193,7 @@ function _calculate_sediment_tendencies!(i, j, sediment::SimpleMultiG, bgc, grid
         O₂  = tracers.O₂[i, j, k]
         NO₃ = tracers.NO₃[i, j, k]
         NH₄ = tracers.NH₄[i, j, k]
-        
+
         A, B, C, D, E, F = sediment.nitrate_oxidation_params
 
         pₙᵢₜ = exp(A +
@@ -208,9 +211,9 @@ function _calculate_sediment_tendencies!(i, j, sediment::SimpleMultiG, bgc, grid
                      sediment.denitrification_params.E * log(reactivity) ^ 2 +
                      sediment.denitrification_params.F * log(O₂) * log(reactivity)) / (Cᵐⁱⁿ * day)
         =#
-        
+
         A, B, C, D, E, F = sediment.anoxic_params
-        
+
         pₐₙₒₓ = exp(A +
                     B * log(Cᵐⁱⁿ * day) +
                     C * log(Cᵐⁱⁿ * day) ^ 2 +
diff --git a/src/Models/seawater_density.jl b/src/Models/seawater_density.jl
index 650638187..6539b859b 100644
--- a/src/Models/seawater_density.jl
+++ b/src/Models/seawater_density.jl
@@ -15,17 +15,17 @@ This function *will not* run on GPU so can not be used as the default.
     Pdbar = 10 * Pbar
 
     Sa = gsw_sa_from_sp(Sp, Pdbar, lon, lat)
-    
+
     return gsw_rho(Sa, T, Pdbar)
 end
 
 """
     teos10_polynomial_approximation(T, Sp, Pbar = 0, lon = 0, lat = 0)
 
-Returns sea water density computed by `SeawaterPolynomials.jl` TEOS10 
-polynomial approximation from `T` in degrees centigrade, `S` in PSU,  
+Returns sea water density computed by `SeawaterPolynomials.jl` TEOS10
+polynomial approximation from `T` in degrees centigrade, `S` in PSU,
 `P`ressure in bar. This approximation also requires us to approximate
-geopotential depth from pressure at 1 m / dbar, and absolute salinity 
+geopotential depth from pressure at 1 m / dbar, and absolute salinity
 as practial salinity.
 
 This method is less accurate than `teos10_density` but will run on GPU.
diff --git a/src/Utils/timestep.jl b/src/Utils/timestep.jl
index f6aa7db6e..3664c4e22 100644
--- a/src/Utils/timestep.jl
+++ b/src/Utils/timestep.jl
@@ -1,5 +1,6 @@
 using Oceananigans.Advection: cell_advection_timescale
-using Oceananigans.Grids: Center, znodes, zspacing, minimum_zspacing
+using Oceananigans.Grids: Center, znodes, minimum_zspacing
+using Oceananigans.Operators: zspacing
 
 @inline column_advection_timescale(model) = minimum_zspacing(model.grid) / maximum_sinking_velocity(model.biogeochemistry)
 
diff --git a/test/test_sediments.jl b/test/test_sediments.jl
index cda5478f7..b4b635b8f 100644
--- a/test/test_sediments.jl
+++ b/test/test_sediments.jl
@@ -1,14 +1,12 @@
 include("dependencies_for_runtests.jl")
 
-using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation
 using Oceananigans.Units
 
-using OceanBioME.Sediments: sediment_tracers, sediment_fields
-using Oceananigans: Field
+using Oceananigans: Field, TendencyCallsite
 using Oceananigans.Fields: TracerFields
-
 using Oceananigans.Operators: volume, Azᶠᶜᶜ
 
+using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation, sediment_tracers, sediment_fields
 using OceanBioME.Models.LOBSTERModel: VariableRedfieldLobster
 
 function intercept_tracer_tendencies!(model, intercepted_tendencies)
@@ -23,22 +21,22 @@ function set_defaults!(sediment::SimpleMultiG)
 
     set!(sediment.fields.C_fast, 0.5893)
     set!(sediment.fields.C_slow, 0.1677)
-end 
+end
 
 set_defaults!(::InstantRemineralisation) = nothing
 
 set_defaults!(::LOBSTER, model) =
-    set!(model, P = 0.4686, Z = 0.5363, 
-                NO₃ = 2.3103, NH₄ = 0.0010, 
+    set!(model, P = 0.4686, Z = 0.5363,
+                NO₃ = 2.3103, NH₄ = 0.0010,
                 DOM = 0.8115,
                 sPOM = 0.2299, bPOM = 0.0103)
 
 
 set_defaults!(::VariableRedfieldLobster, model) =
-    set!(model, P = 0.4686, Z = 0.5363, 
-                NO₃ = 2.3103, NH₄ = 0.0010, 
-                DIC = 2106.9, Alk = 2408.9, 
-                O₂ = 258.92, 
+    set!(model, P = 0.4686, Z = 0.5363,
+                NO₃ = 2.3103, NH₄ = 0.0010,
+                DIC = 2106.9, Alk = 2408.9,
+                O₂ = 258.92,
                 DOC = 5.3390, DON = 0.8115,
                 sPON = 0.2299, sPOC = 1.5080,
                 bPON = 0.0103, bPOC = 0.0781)
@@ -46,8 +44,8 @@ set_defaults!(::VariableRedfieldLobster, model) =
 
 set_defaults!(::NutrientPhytoplanktonZooplanktonDetritus, model) =  set!(model, N = 2.3, P = 0.4, Z = 0.5, D = 0.2)
 
-total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) + 
-                                    sum(sed.fields.N_slow) + 
+total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) +
+                                    sum(sed.fields.N_slow) +
                                     sum(sed.fields.N_ref)
 
 total_nitrogen(sed::InstantRemineralisation) = sum(sed.fields.N_storage)
@@ -67,20 +65,20 @@ total_nitrogen(::VariableRedfieldLobster, model) = sum(model.tracers.NO₃) +
                                                      sum(model.tracers.DON) +
                                                      sum(model.tracers.sPON) +
                                                      sum(model.tracers.bPON)
-                                                     
+
 total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(model.tracers.N) +
                                                                     sum(model.tracers.P) +
                                                                     sum(model.tracers.Z) +
                                                                     sum(model.tracers.D)
 
 function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAdamsBashforth2)
-    model = isa(model, NonhydrostaticModel) ? model(; grid, 
-                                                      biogeochemistry, 
+    model = isa(model, NonhydrostaticModel) ? model(; grid,
+                                                      biogeochemistry,
                                                       closure = nothing,
                                                       timestepper,
                                                       buoyancy = nothing) :
-                                              model(; grid, 
-                                                      biogeochemistry, 
+                                              model(; grid,
+                                                      biogeochemistry,
                                                       closure = nothing,
                                                       buoyancy = nothing,
                                                       tracers = nothing)
@@ -114,7 +112,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd
 
     # conservations
     rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7)
-    @test isapprox(N₀, N₁; rtol) 
+    @test isapprox(N₀, N₁; rtol)
 
     return nothing
 end
@@ -131,30 +129,82 @@ display_name(::ImmersedBoundaryGrid) = "Immersed boundary grid"
 bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill
 
 @testset "Sediment integration" begin
-    grids = [RectilinearGrid(architecture; size=(3, 3, 50), extent=(10, 10, 500)),
-             LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)),
-             ImmersedBoundaryGrid(
-                LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), 
-                GridFittedBottom(bottom_height))]
-    for grid in grids
-        for timestepper in (:QuasiAdamsBashforth2, :RungeKutta3),
-            sediment_model in (InstantRemineralisation(; grid), SimpleMultiG(; grid)),
-            model in (NonhydrostaticModel, HydrostaticFreeSurfaceModel)
-            for biogeochemistry in (NutrientPhytoplanktonZooplanktonDetritus(; grid, sediment_model),
-                                    LOBSTER(; grid,
-                                              carbonates = ifelse(isa(sediment_model, SimpleMultiG), true, false), 
-                                              oxygen = ifelse(isa(sediment_model, SimpleMultiG), true, false), 
-                                              variable_redfield = ifelse(isa(sediment_model, SimpleMultiG), true, false), 
-                                              sediment_model))
-                # get rid of incompatible combinations
-                run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) ||
-                             (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) ||
-                             (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true)
-                
-                if run
-                    @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))"
-                    @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" test_flat_sediment(grid, biogeochemistry, model; timestepper)
-                end
+    rectilinear_grid = RectilinearGrid(
+        architecture;
+        size = (3, 3, 50),
+        extent = (10, 10, 500)
+    )
+
+    latlon_grid = LatitudeLongitudeGrid(
+        architecture;
+        size = (3, 3, 16),
+        latitude = (0, 10),
+        longitude = (0, 10),
+        z = (-500, 0)
+    )
+
+    underlying_latlon_grid = LatitudeLongitudeGrid(
+        architecture;
+        size = (3, 3, 16),
+        latitude = (0, 10),
+        longitude = (0, 10),
+        z = (-500, 0)
+    )
+
+    immersed_latlon_grid = ImmersedBoundaryGrid(
+        underlying_latlon_grid,
+        GridFittedBottom(bottom_height)
+    )
+
+    grids = (rectilinear_grid, latlon_grid, underlying_latlon_grid)
+    timesteppers = (:QuasiAdamsBashforth2, :RungeKutta3)
+    sediment_models = (InstantRemineralisation(; grid), SimpleMultiG(; grid))
+    models = (NonhydrostaticModel, HydrostaticFreeSurfaceModel)
+
+    for grid in grids, timestepper in timesteppers, sediment_model in sediment_models, model in models
+        npzd_bgc_model = NutrientPhytoplanktonZooplanktonDetritus(; grid, sediment_model)
+
+        using_simple_multi_g = sediment_model isa SimpleMultiG
+
+        lobster_bgc_model = LOBSTER(;
+            grid,
+            sediment_model,
+            carbonates = using_simple_multi_g,
+            oxygen = using_simple_multi_g,
+            variable_redfield = using_simple_multi_g,
+        )
+
+        bgc_models = [npzd_bgc_model, lobster_bgc_model]
+
+        for biogeochemistry in bgc_models
+            nonhydrostatic = (model == NonhydrostaticModel)
+            hydrostatic = (model == HydrostaticFreeSurfaceModel)
+
+            grid_is_immersed = grid isa ImmersedBoundaryGrid
+            grid_is_latlon = grid isa LatitudeLongitudeGrid
+
+            rk3 = (timestepper == :RungeKutta3)
+            simple_multi_g = sediment_model isa SimpleMultiG
+            bgc_is_npzd = biogeochemistry.underlying_biogeochemistry isa NutrientPhytoplanktonZooplanktonDetritus
+
+            # Skip incompatible combinations
+            skip1 = nonhydrostatic && (grid_is_immersed || grid_is_latlon)
+            skip2 = hydrostatic && rk3
+            skip3 = simple_multi_g && bgc_is_npzd
+
+            if skip1 || skip2 || skip3
+                continue
+            end
+
+            arch_name = typeof(architecture)
+            sediment_name = display_name(sediment_model)
+            bgc_name = display_name(biogeochemistry.underlying_biogeochemistry)
+            grid_name = display_name(grid)
+
+            @info "Testing sediment on $arch_name with $timestepper and $sediment_name on $bgc_name with $grid_name"
+
+            @testset "$architecture, $timestepper, $sediment_name, $bgc_name, $grid_name" begin
+                @test_skip test_flat_sediment(grid, biogeochemistry, model; timestepper)
             end
         end
     end
diff --git a/test/test_utils.jl b/test/test_utils.jl
index b7c1c32b1..0815885a8 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -2,7 +2,7 @@ include("dependencies_for_runtests.jl")
 
 using OceanBioME: setup_velocity_fields, valid_sinking_velocity_locations
 
-using Oceananigans.Fields: AbstractField, CenterField, ConstantField, FunctionField, ZFaceField
+using Oceananigans.Fields: AbstractField, CenterField, ConstantField, FunctionField, ZFaceField, location
 
 function test_negative_scaling(arch)
     grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1))
@@ -57,7 +57,7 @@ sinking_speeds = merge(scalar_sinking_speeds, field_sinking_speeds)
 
     @test all(map(w -> isa(w, AbstractField) & (location(w) in valid_sinking_velocity_locations), values(sinking_velocities)))
 
-    sinking_velocities = 
+    sinking_velocities =
         @test_warn ("The sinking velocity provided for C is a field and therefore `open_bottom=false` can't be enforced automatically",
                     "The sinking velocity provided for D is a field and therefore `open_bottom=false` can't be enforced automatically",
                     "The sinking velocity provided for E is a field and therefore `open_bottom=false` can't be enforced automatically") setup_velocity_fields(sinking_speeds, grid, false)