diff --git a/Manifest.toml b/Manifest.toml
index 8babc3a20..f7baaf529 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -17,9 +17,9 @@ weakdeps = ["ChainRulesCore", "Test"]
 
 [[deps.Accessors]]
 deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Markdown", "Test"]
-git-tree-sha1 = "c0d491ef0b135fd7d63cbc6404286bc633329425"
+git-tree-sha1 = "f61b15be1d76846c0ce31d3fcfac5380ae53db6a"
 uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.36"
+version = "0.1.37"
 
     [deps.Accessors.extensions]
     AccessorsAxisKeysExt = "AxisKeys"
@@ -52,9 +52,9 @@ version = "1.1.1"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "ed2ec3c9b483842ae59cd273834e5b46206d6dda"
+git-tree-sha1 = "5c9b74c973181571deb6442d41e5c902e6b9f38e"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.11.0"
+version = "7.12.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -137,9 +137,9 @@ version = "5.4.2"
 
 [[deps.CUDA_Driver_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
-git-tree-sha1 = "c48f9da18efd43b6b7adb7ee1f93fe5f2926c339"
+git-tree-sha1 = "97df9d4d6be8ac6270cb8fd3b8fc413690820cbd"
 uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc"
-version = "0.9.0+0"
+version = "0.9.1+1"
 
 [[deps.CUDA_Runtime_Discovery]]
 deps = ["Libdl"]
@@ -149,9 +149,9 @@ version = "0.3.4"
 
 [[deps.CUDA_Runtime_jll]]
 deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "bcba305388e16aa5c879e896726db9e71b4942c6"
+git-tree-sha1 = "afea94249b821dc754a8ca6695d3daed851e1f5a"
 uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
-version = "0.14.0+1"
+version = "0.14.1+0"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
@@ -217,9 +217,9 @@ weakdeps = ["InverseFunctions"]
 
 [[deps.ConstructionBase]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "260fd2400ed2dab602a7c15cf10c1933c59930a2"
+git-tree-sha1 = "d8a9c0b6ac2d9081bf76324b39c78ca3ce4f0c98"
 uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.5.5"
+version = "1.5.6"
 
     [deps.ConstructionBase.extensions]
     ConstructionBaseIntervalSetsExt = "IntervalSets"
@@ -346,9 +346,9 @@ version = "6.2.1+6"
 
 [[deps.GPUArrays]]
 deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]
-git-tree-sha1 = "5c9de6d5af87acd2cf719e214ed7d51e14017b7a"
+git-tree-sha1 = "a74c3f1cf56a3dfcdef0605f8cdb7015926aae30"
 uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
-version = "10.2.2"
+version = "10.3.0"
 
 [[deps.GPUArraysCore]]
 deps = ["Adapt"]
@@ -381,9 +381,9 @@ version = "1.14.2+1"
 
 [[deps.Hwloc_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1d334207121865ac8c1c97eb7f42d0339e4635bf"
+git-tree-sha1 = "5e19e1e4fa3e71b774ce746274364aef0234634e"
 uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
-version = "2.11.0+0"
+version = "2.11.1+0"
 
 [[deps.IfElse]]
 git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
@@ -397,22 +397,23 @@ uuid = "40713840-3770-5561-ab4c-a76e7d0d7895"
 version = "0.2.1"
 
 [[deps.InlineStrings]]
-deps = ["Parsers"]
-git-tree-sha1 = "86356004f30f8e737eff143d57d41bd580e437aa"
+git-tree-sha1 = "45521d31238e87ee9f9732561bfee12d4eebd52d"
 uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
-version = "1.4.1"
+version = "1.4.2"
 
     [deps.InlineStrings.extensions]
     ArrowTypesExt = "ArrowTypes"
+    ParsersExt = "Parsers"
 
     [deps.InlineStrings.weakdeps]
     ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd"
+    Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 
 [[deps.IntelOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "be50fe8df3acbffa0274a744f1a99d29c45a57f4"
+git-tree-sha1 = "14eb2b542e748570b56446f4c50fbfb2306ebc45"
 uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2024.1.0+0"
+version = "2024.2.0+0"
 
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
@@ -420,9 +421,9 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.InverseFunctions]]
 deps = ["Test"]
-git-tree-sha1 = "e7cbed5032c4c397a6ac23d1493f3289e01231c4"
+git-tree-sha1 = "18c59411ece4838b18cd7f537e56cf5e41ce5bfd"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.14"
+version = "0.1.15"
 weakdeps = ["Dates"]
 
     [deps.InverseFunctions.extensions]
@@ -600,9 +601,9 @@ version = "1.9.4+0"
 
 [[deps.MKL_jll]]
 deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"]
-git-tree-sha1 = "80b2833b56d466b3858d565adcd16a4a05f2089b"
+git-tree-sha1 = "f046ccd0c6db2832a9f639e2c669c6fe867e5f4f"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2024.1.0+0"
+version = "2024.2.0+0"
 
 [[deps.MPI]]
 deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"]
@@ -620,9 +621,9 @@ version = "0.20.16"
 
 [[deps.MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "4099bb6809ac109bfc17d521dad33763bcf026b7"
+git-tree-sha1 = "19d4bd098928a3263693991500d05d74dbdc2004"
 uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
-version = "4.2.1+1"
+version = "4.2.2+0"
 
 [[deps.MPIPreferences]]
 deps = ["Libdl", "Preferences"]
@@ -717,9 +718,9 @@ version = "0.91.3"
     Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e"
+git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.14.0"
+version = "1.14.1"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -737,9 +738,9 @@ version = "0.8.1+2"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"]
-git-tree-sha1 = "a9de2f1fc98b92f8856c640bf4aec1ac9b2a0d86"
+git-tree-sha1 = "2f0a1d8c79bc385ec3fcda12830c9d0e72b30e71"
 uuid = "fe0851c0-eecd-5654-98d4-656369965a5c"
-version = "5.0.3+0"
+version = "5.0.4+0"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -937,9 +938,9 @@ version = "0.3.4"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "90b4f68892337554d31cdcdbe19e48989f26c7e6"
+git-tree-sha1 = "ff11acffdb082493657550959d4feb4b6149e73a"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.3"
+version = "1.4.5"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -970,9 +971,9 @@ weakdeps = ["ChainRulesCore"]
 
 [[deps.Static]]
 deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"]
-git-tree-sha1 = "0bbff21027dd8a107551847528127b62a35f7594"
+git-tree-sha1 = "87d51a3ee9a4b0d2fe054bdd3fc2436258db2603"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "1.1.0"
+version = "1.1.1"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
@@ -987,9 +988,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"]
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "20833c5b7f7edf0e5026f23db7f268e4f23ec577"
+git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.6"
+version = "1.9.7"
 weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
@@ -1079,10 +1080,10 @@ uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
 version = "1.0.1"
 
 [[deps.Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits"]
-git-tree-sha1 = "cb76cf677714c095e535e3501ac7954732aeea2d"
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"]
+git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.11.1"
+version = "1.12.0"
 
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
diff --git a/Project.toml b/Project.toml
index 7781e1e0b..b1951ea73 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.10.4"
+version = "0.10.5"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/benchmark/box_model.jl b/benchmark/box_model.jl
index 21dfebc22..3be049801 100644
--- a/benchmark/box_model.jl
+++ b/benchmark/box_model.jl
@@ -17,6 +17,7 @@
 # ```
 
 using OceanBioME, Oceananigans.Units, Oceananigans, BenchmarkTools
+using Oceananigans.Fields: FunctionField
 
 const year = years = 365day
 
@@ -25,13 +26,20 @@ const year = years = 365day
 @inline PAR⁰(t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
 
 z = -10 # nominal depth of the box for the PAR profile
-@inline PAR(t)::Float64 = PAR⁰(t) * exp(0.2z) # Modify the PAR based on the nominal depth and exponential decay 
+@inline PAR_func(t)::Float64 = PAR⁰(t) * exp(0.2z) # Modify the PAR based on the nominal depth and exponential decay 
 
 function run_box_simulation()
-    biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid,
-                                                                 light_attenuation_model = nothing)
 
-    model = BoxModel(; biogeochemistry, forcing = (; PAR))
+    clock = Clock(time = Float64(0))
+
+    grid = BoxModelGrid()
+
+    PAR = FunctionField{Center, Center, Center}(PAR_func, grid; clock)
+
+    biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid,
+                                                                 light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation(PAR))
+
+    model = BoxModel(; biogeochemistry, clock)
 
     set!(model, N = 10.0, P = 0.1, Z = 0.01)
 
@@ -41,25 +49,12 @@ function run_box_simulation()
 
     fast_output = SpeedyOutput("box_benchmarking.jld2")
 
-    simulation.callbacks[:output] = Callback(fast_output, IterationInterval(20);)
+    simulation.callbacks[:output] = Callback(fast_output, IterationInterval(20))
 
     @info "Running the model..."
     run!(simulation)
 end
 
-function fast_output(sim, fname)
-    file = jldopen(fname, "w+")
-
-    model = sim.model
-
-    t = time(sim)
-
-    file["fields/$t"] = model.field_values
-
-    close(file)
-
-    return nothing
-end
 #####
 ##### results
 #####
diff --git a/docs/make.jl b/docs/make.jl
index 8374327e6..66bf1f6dd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -52,12 +52,12 @@ example_pages = [ title => "generated/$(filename).md" for (title, filename) in e
 
 if !isdir(OUTPUT_DIR) mkdir(OUTPUT_DIR) end
 
-model_parameters = (LOBSTER(; grid = BoxModelGrid, light_attenuation_model = nothing).underlying_biogeochemistry,
-                    NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid, light_attenuation_model = nothing).underlying_biogeochemistry,
+model_parameters = (LOBSTER(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry,
+                    NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry,
                     SLatissima(),
                     TwoBandPhotosyntheticallyActiveRadiation(; grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))),
-                    SimpleMultiG(; grid = BoxModelGrid),
-                    InstantRemineralisation(; grid = BoxModelGrid),
+                    SimpleMultiG(; grid = BoxModelGrid()),
+                    InstantRemineralisation(; grid = BoxModelGrid()),
                     OCMIP_default,
                     GasExchange(; gas = :CO₂).condition.func,
                     GasExchange(; gas = :O₂).condition.func)
diff --git a/examples/box.jl b/examples/box.jl
index 4574c0641..85f7fe7a1 100644
--- a/examples/box.jl
+++ b/examples/box.jl
@@ -17,7 +17,7 @@ using Oceananigans.Fields: FunctionField
 
 const year = years = 365day
 
-grid = BoxModelGrid
+grid = BoxModelGrid()
 clock = Clock(time = zero(grid))
 
 nothing #hide
diff --git a/examples/data_assimilation.jl b/examples/data_assimilation.jl
index cf62b6bb7..c30212ea0 100644
--- a/examples/data_assimilation.jl
+++ b/examples/data_assimilation.jl
@@ -42,7 +42,7 @@ function run_box_simulation(initial_photosynthetic_slope,
                             phyto_base_mortality_rate,
                             j)
 
-    biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid,
+    biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid(),
                                                                  initial_photosynthetic_slope,
                                                                  base_maximum_growth,
                                                                  nutrient_half_saturation,
diff --git a/src/BoxModel/boxmodel.jl b/src/BoxModel/boxmodel.jl
index 1a6546461..48abc5c8f 100644
--- a/src/BoxModel/boxmodel.jl
+++ b/src/BoxModel/boxmodel.jl
@@ -3,7 +3,7 @@ Integrate biogeochemical models on a single point
 """
 module BoxModels
 
-export BoxModel, run!, set!, SpeedyOutput
+export BoxModel, run!, set!, SpeedyOutput, load_output
 
 using Oceananigans: Clock, prettytime
 using Oceananigans.Biogeochemistry: 
@@ -48,36 +48,37 @@ end
                clock = Clock(; time = 0.0),
                prescribed_tracers = (:T, ))
 
-Constructs a box model of a `biogeochemistry` model. Once this has been constructed you can set initial condiitons by `set!(model, X=1.0...)` and then `run!(model)`.
+Constructs a box model of a `biogeochemistry` model. Once this has been constructed you can set initial condiitons by `set!(model, X=1.0...)`.
 
 Keyword Arguments 
 =================
 
-- `biogeochemistry`: (required) an OceanBioME biogeochemical model, most models must be passed a `grid` which can be set to `BoxModelGrid` for box models
+- `biogeochemistry`: (required) an OceanBioME biogeochemical model, most models must be passed a `grid` which can be set to a `BoxModelGrid` for box models
 - `forcing`: NamedTuple of additional forcing functions for the biogeochemical tracers to be integrated
 - `timestepper`: Timestepper to integrate model
 - `clock`: Oceananigans clock to keep track of time
 - `prescribed_tracers`: Tuple of fields names (Symbols) which are not integrated but provided in `forcing` as a function of time with signature `f(t)`
 """
 function BoxModel(; biogeochemistry::B,
+                    grid = BoxModelGrid(),
                     forcing = NamedTuple(),
                     timestepper = :RungeKutta3,
                     clock::C = Clock(; time = 0.0),
                     prescribed_tracers::PT = (T = (t) -> 0, )) where {B, C, PT}
 
     variables = required_biogeochemical_tracers(biogeochemistry)
-    fields = NamedTuple{variables}([CenterField(BoxModelGrid) for var in eachindex(variables)])
+    fields = NamedTuple{variables}([CenterField(grid) for var in eachindex(variables)])
     field_values = zeros(length(variables)+length(required_biogeochemical_auxiliary_fields(biogeochemistry)))
     forcing = NamedTuple{variables}([variable in keys(forcing) ? forcing[variable] : no_func for variable in variables])
 
-    timestepper = BoxModelTimeStepper(timestepper, BoxModelGrid, keys(fields))
+    timestepper = BoxModelTimeStepper(timestepper, grid, keys(fields))
 
     F = typeof(fields)
     FV = typeof(field_values)
     FO = typeof(forcing)
     TS = typeof(timestepper)
 
-    return BoxModel{typeof(BoxModelGrid), B, F, FV, FO, TS, C, PT}(BoxModelGrid, biogeochemistry, fields, field_values, forcing, timestepper, clock, prescribed_tracers)
+    return BoxModel{typeof(grid), B, F, FV, FO, TS, C, PT}(grid, biogeochemistry, fields, field_values, forcing, timestepper, clock, prescribed_tracers)
 end
 
 function update_state!(model::BoxModel, callbacks=[]; compute_tendencies = true)
@@ -101,7 +102,7 @@ function update_state!(model::BoxModel, callbacks=[]; compute_tendencies = true)
     return nothing
 end
 
-architecture(::BoxModel) = architecture(BoxModelGrid)
+architecture(model::BoxModel) = architecture(model.grid) # this might be the default fallback
 default_nan_checker(::BoxModel) = nothing
 iteration(model::BoxModel) = model.clock.iteration
 prognostic_fields(model::BoxModel) = @inbounds model.fields[required_biogeochemical_tracers(model.biogeochemistry)]
@@ -139,13 +140,13 @@ default_included_properties(::BoxModel) = [:grid]
 include("timesteppers.jl")
 include("output_writer.jl")
 
-summary(::BoxModel{B, V, F, TS, C}) where {B, V, F, TS, C} = string("Biogeochemical box model")
-show(io::IO, model::BoxModel{B, V, F, TS, C}) where {B, V, F, TS, C} = 
+summary(::BoxModel)  = string("Biogeochemical box model")
+show(io::IO, model::BoxModel{G, B, F, FV, FO, TS, C, PT}) where {G, B, F, FV, FO, TS, C, PT} = 
        print(io, summary(model), "\n",
                 "  Biogeochemical model: ", "\n",
                 "    └── ", summary(model.biogeochemistry), "\n",
                 "  Time-stepper:", "\n", 
-                "    └── ", summary(model.timestepper), "\n",
+                "    └── ", nameof(typeof(model.timestepper)), "\n",
                 "  Time:", "\n",
                 "    └── $(prettytime(model.clock.time))")
 
diff --git a/src/BoxModel/output_writer.jl b/src/BoxModel/output_writer.jl
index 519faccf9..dc7e26470 100644
--- a/src/BoxModel/output_writer.jl
+++ b/src/BoxModel/output_writer.jl
@@ -3,6 +3,10 @@ struct SpeedyOutput{FN, B} # I will make this an `AbstractOutputWriter` at some
     overwrite_existing :: B
 
     function SpeedyOutput(filename::FN; overwrite_existing::B = true) where {FN, B}
+        jldopen(filename, ifelse(overwrite_existing, "w+", "w")) do file
+            file["filename"] = filename
+        end
+
         return new{FN, B}(filename, overwrite_existing)
     end
 end
@@ -13,15 +17,37 @@ function (save::SpeedyOutput)(simulation)
     t = time(simulation)
     iter = model.clock.iteration
 
-    file = jldopen(save.filename, ifelse(save.overwrite_existing, "w+", "w"))
+    jldopen(save.filename, "a+") do file
 
-    file["timeseries/t/$iter"] = t
+        file["timeseries/t/$iter"] = t
 
-    for (i, name) in enumerate(keys(model.fields))
-        file["timeseries/$name/$iter"] = model.field_values[i]
+        for (i, name) in enumerate(keys(model.fields))
+            file["timeseries/$name/$iter"] = model.field_values[i]
+        end
     end
 
-    close(file)
-
     return nothing
-end
\ No newline at end of file
+end
+
+function load_output(save::SpeedyOutput)
+    names = tuple()
+    
+    jldopen(save.filename) do f 
+        names = Symbol.(keys(f["timeseries"]))
+    end
+
+    sort_order = sortperm(load_output(save, "t"))
+
+    return NamedTuple{tuple(names...)}([load_output(save, name) for name in names])
+end
+    
+function load_output(save::SpeedyOutput, name)
+    sort_order = output_order(save)
+
+    return Float64.(values(load(save.filename, "timeseries/"*string(name))))[sort_order]
+end
+
+load_output(save::SpeedyOutput, names::Tuple) = 
+    NamedTuple{name}([load_output(save, name) for name in names])
+
+output_order(save) = sortperm(Float64.(values(load(save.filename, "timeseries/t"))))
\ No newline at end of file
diff --git a/src/Light/prescribed.jl b/src/Light/prescribed.jl
index 44c4d385c..b3949977e 100644
--- a/src/Light/prescribed.jl
+++ b/src/Light/prescribed.jl
@@ -46,8 +46,10 @@ function update_biogeochemical_state!(model, PAR::PrescribedPhotosyntheticallyAc
     return nothing
 end
 
-summary(::PrescribedPhotosyntheticallyActiveRadiation{FT}) where {FT} = string("Prescribed PAR")
-show(io::IO, model::PrescribedPhotosyntheticallyActiveRadiation{FT}) where {FT} = print(io, summary(model))
+summary(::PrescribedPhotosyntheticallyActiveRadiation) = string("Prescribed PAR")
+show(io::IO, model::PrescribedPhotosyntheticallyActiveRadiation{F}) where {F} = print(io, summary(model), "\n",
+                                                                                            "  Fields:", "\n",
+                                                                                            "    └── $(model.field_names)")
 
 biogeochemical_auxiliary_fields(par::PrescribedPhotosyntheticallyActiveRadiation) = NamedTuple{par.field_names}(par.fields)
 
diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl
index 1df55b2a1..fabea3538 100644
--- a/src/OceanBioME.jl
+++ b/src/OceanBioME.jl
@@ -11,7 +11,7 @@ export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, NPZD,
 export SLatissima
 
 # Box model
-export BoxModel, BoxModelGrid, SpeedyOutput, run!, set!
+export BoxModel, BoxModelGrid, SpeedyOutput, load_output
 
 # Particles
 export Particles
@@ -31,7 +31,7 @@ export ScaleNegativeTracers, ZeroNegativeTracers
 # Oceananigans extensions
 export ColumnField, isacolumn
 
-using Oceananigans.Architectures: architecture, device
+using Oceananigans.Architectures: architecture, device, CPU
 using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry
 using Oceananigans.Grids: RectilinearGrid, Flat
 
@@ -132,7 +132,7 @@ end
 
 update_biogeochemical_state!(model, modifiers::Tuple) = [update_biogeochemical_state!(model, modifier) for modifier in modifiers]
 
-const BoxModelGrid = RectilinearGrid(topology = (Flat, Flat, Flat), size = ())
+BoxModelGrid(; arch = CPU(), kwargs...) = RectilinearGrid(arch; topology = (Flat, Flat, Flat), size = (), kwargs...)
 
 @inline maximum_sinking_velocity(bgc) = 0.0
 
diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl
index 498130dcc..80c384f81 100644
--- a/src/Utils/negative_tracers.jl
+++ b/src/Utils/negative_tracers.jl
@@ -1,5 +1,5 @@
 using Oceananigans: fields, Simulation
-using KernelAbstractions
+using KernelAbstractions: @kernel, @index
 using Oceananigans.Utils: work_layout
 using Oceananigans.Architectures: device, architecture, on_architecture
 using Oceananigans.Biogeochemistry: AbstractBiogeochemistry
diff --git a/test/dependencies_for_runtests.jl b/test/dependencies_for_runtests.jl
index 5bcca955a..6b6250ad8 100644
--- a/test/dependencies_for_runtests.jl
+++ b/test/dependencies_for_runtests.jl
@@ -1,3 +1,3 @@
-using OceanBioME, Test, CUDA, Oceananigans, JLD2, Documenter
+using OceanBioME, Test, CUDA, Oceananigans, JLD2, Oceananigans.Units, Documenter
 
 architecture = CUDA.has_cuda() ? GPU() : CPU()
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index f18443d60..04cad701e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,6 +8,10 @@ include("test_NPZD.jl")
 include("test_gasexchange.jl")
 include("test_sediments.jl")
 
+if architecture == CPU() 
+    include("test_boxmodel.jl") # box models (probably) don't work on GPU, and it wouldn't be faster anyway
+end
+
 architecture == CPU() && @testset "Doctests" begin
     doctest(OceanBioME)
 end
diff --git a/test/test_boxmodel.jl b/test/test_boxmodel.jl
new file mode 100644
index 000000000..d55fc6721
--- /dev/null
+++ b/test/test_boxmodel.jl
@@ -0,0 +1,72 @@
+#include("dependencies_for_runtests.jl")
+
+using Oceananigans: AbstractGrid, AbstractModel, AbstractOutputWriter
+using Oceananigans.Fields: FunctionField
+
+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()
+    clock = Clock(time = zero(grid))
+
+    light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation(FunctionField{Center, Center, Center}(PAR, grid; clock))
+
+    biogeochemistry = LOBSTER(; grid, light_attenuation_model)
+
+    model = BoxModel(; grid, biogeochemistry, clock)
+
+    return model
+end
+
+@testset "Construct `BoxModel`" begin
+    model = simple_box_model()
+
+    @test grid isa AbstractGrid
+    @test model isa AbstractModel
+    @test set!(model; defaults...) isa Nothing
+end
+
+@testset "Timestep `BoxModel`" begin
+    model = simple_box_model()
+
+    set!(model; defaults...) 
+
+    for n in 1:10
+        time_step!(model, 20minutes)
+    end
+
+    # values are being updated
+    for (name, field) in enumerate(model.fields)
+        default = name in keys(defaults) ? defaults[name] : 0
+        @test field[1, 1, 1] != default
+    end
+end
+
+@testset "`BoxModel` simulation" begin
+    model = simple_box_model()
+
+    set!(model; defaults...) 
+
+    simulation = Simulation(model; Δt = 20minutes, stop_iteration = 1000, verbose = false)
+
+    fname = "box_model_test.jld2"
+
+    fast_output = SpeedyOutput(fname)
+
+    @test_broken fast_output isa AbstractOutputWriter
+
+    simulation.callbacks[:output] = Callback(fast_output, IterationInterval(20))
+
+    @test run!(simulation) isa Nothing
+
+    results = load_output(fast_output)
+
+    # outputting as expected
+    @test length(keys(results)) == length(model.fields) + 1
+    @test length(results.t) == 1000 / 20 + 1
+
+    # something is happening and we're not just writing the same value 
+    @test all([result[1] != result[end] for result in values(results)])
+end
\ No newline at end of file