diff --git a/Project.toml b/Project.toml
index 734e3e1ac1..dc403c2e92 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.95.6"
+version = "0.95.7"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/Biogeochemistry.jl b/src/Biogeochemistry.jl
index 10ecc3f8e1..b0b27ef505 100644
--- a/src/Biogeochemistry.jl
+++ b/src/Biogeochemistry.jl
@@ -33,7 +33,7 @@ update_biogeochemical_state!(bgc, model) = nothing
 """
     AbstractBiogeochemistry
 
-Abstract type for biogeochemical models. To define a biogeochemcial relaionship
+Abstract type for biogeochemical models. To define a biogeochemcial relationship
 the following functions must have methods defined where `BiogeochemicalModel`
 is a subtype of `AbstractBioeochemistry`:
 
diff --git a/src/BoundaryConditions/continuous_boundary_function.jl b/src/BoundaryConditions/continuous_boundary_function.jl
index 922739641b..3055864d69 100644
--- a/src/BoundaryConditions/continuous_boundary_function.jl
+++ b/src/BoundaryConditions/continuous_boundary_function.jl
@@ -72,7 +72,7 @@ The regularization of `bc.condition::ContinuousBoundaryFunction` requries
    of the boundary.
 """
 function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction},
-                                       grid, loc, dim, Side, prognostic_field_names) where C
+                                       grid, loc, dim, Side, field_names) where C
 
     boundary_func = bc.condition
 
@@ -81,7 +81,7 @@ function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoun
 
     indices, interps = index_and_interp_dependencies(LX, LY, LZ,
                                                      boundary_func.field_dependencies,
-                                                     prognostic_field_names)
+                                                     field_names)
 
     regularized_boundary_func = ContinuousBoundaryFunction{LX, LY, LZ, Side}(boundary_func.func,
                                                                              boundary_func.parameters,
diff --git a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl
index 2c17fc47d5..dca5dc6a0e 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl
@@ -75,7 +75,7 @@ Return a flattened `NamedTuple` of the fields in `model.velocities`, `model.free
 `model.tracers`, and any auxiliary fields for a `HydrostaticFreeSurfaceModel` model.
 """
 @inline fields(model::HydrostaticFreeSurfaceModel) = 
-    merge(hydrostatic_fields(model.velocities, model.free_surface, model.tracers),
+    merge(hydrostatic_fields(model.velocities, model.free_surface, model.tracers), 
           model.auxiliary_fields,
           biogeochemical_auxiliary_fields(model.biogeochemistry))
 
@@ -107,20 +107,25 @@ Return a flattened `NamedTuple` of the prognostic fields associated with `Hydros
                                                                        v = velocities.v,
                                                                        w = velocities.w),
                                                                        tracers,
-                                                                       (; η = free_surface.η))
+                                                                       (; η = free_surface.η,
+                                                                          U = nothing,
+                                                                          V = nothing))
 
 @inline hydrostatic_fields(velocities, free_surface::SplitExplicitFreeSurface, tracers) = merge((u = velocities.u,
                                                                                                  v = velocities.v,
-                                                                                                 w = velocities.w,
-                                                                                                 η = free_surface.η,
-                                                                                                 U = free_surface.barotropic_velocities.U,
-                                                                                                 V = free_surface.barotropic_velocities.V),
-                                                                                                 tracers)
+                                                                                                 w = velocities.w),
+                                                                                                 tracers,
+                                                                                                 (; η = free_surface.η,
+                                                                                                    U = free_surface.barotropic_velocities.U,
+                                                                                                    V = free_surface.barotropic_velocities.V))
 
 @inline hydrostatic_fields(velocities, ::Nothing, tracers) = merge((u = velocities.u,
                                                                     v = velocities.v,
                                                                     w = velocities.w),
-                                                                    tracers)
+                                                                    tracers,
+                                                                    (; η = nothing, 
+                                                                       U = nothing,
+                                                                       V = nothing))
 
 displacement(free_surface) = free_surface.η
 displacement(::Nothing) = nothing
diff --git a/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/split_explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/split_explicit_free_surface.jl
index 9052c1eb7f..0cfe0e602b 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/split_explicit_free_surface.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/split_explicit_free_surface.jl
@@ -137,6 +137,10 @@ function split_explicit_substepping(cfl, ::Nothing, fixed_Δt, grid, averaging_k
     return substepping
 end
 
+# Disambiguation for default `SplitExplicitFreeSurface` constructors
+split_explicit_substepping(::Nothing, ::Nothing, ::Nothing, grid, averaging_kernel, gravitational_acceleration) = 
+    split_explicit_substepping(nothing, MINIMUM_SUBSTEPS, nothing, grid, averaging_kernel, gravitational_acceleration)
+
 # TODO: When open boundary conditions are online
 # We need to calculate the barotropic boundary conditions 
 # from the baroclinic boundary conditions by integrating the BC upwards
diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl
index 354b273017..588d322eb6 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl
@@ -128,8 +128,8 @@ function HydrostaticFreeSurfaceModel(; grid,
 
     # Validate biogeochemistry (add biogeochemical tracers automagically)
     tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example)
-    biogeochemical_fields = merge(auxiliary_fields, biogeochemical_auxiliary_fields(biogeochemistry))
-    tracers, auxiliary_fields = validate_biogeochemistry(tracers, biogeochemical_fields, biogeochemistry, grid, clock)
+    biogeochemical_fields = biogeochemical_auxiliary_fields(biogeochemistry)
+    tracers, biogeochemical_fields = validate_biogeochemistry(tracers, biogeochemical_fields, biogeochemistry, grid, clock)
 
     # Reduce the advection order in directions that do not have enough grid points
     @apply_regionally momentum_advection = validate_momentum_advection(momentum_advection, grid)
@@ -158,14 +158,17 @@ function HydrostaticFreeSurfaceModel(; grid,
                                          extract_boundary_conditions(diffusivity_fields))
 
     # Next, we form a list of default boundary conditions:
-    prognostic_field_names = (:u, :v, :w, tracernames(tracers)..., :η, keys(auxiliary_fields)...)
-    default_boundary_conditions = NamedTuple{prognostic_field_names}(Tuple(FieldBoundaryConditions()
-                                                                           for name in prognostic_field_names))
+    field_names = (:u, :v, :w, tracernames(tracers)..., :η, :U, :V, 
+                   keys(auxiliary_fields)..., 
+                   keys(biogeochemical_auxiliary_fields(biogeochemistry))...)
+                   
+    default_boundary_conditions = NamedTuple{field_names}(Tuple(FieldBoundaryConditions()
+                                                          for name in field_names))
 
     # Then we merge specified, embedded, and default boundary conditions. Specified boundary conditions
     # have precedence, followed by embedded, followed by default.
     boundary_conditions = merge(default_boundary_conditions, embedded_boundary_conditions, boundary_conditions)
-    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names)
+    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, field_names)
 
     # Finally, we ensure that closure-specific boundary conditions, such as
     # those required by CATKEVerticalDiffusivity, are enforced:
diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl
index 28cb289e58..af1771bc71 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl
@@ -120,7 +120,9 @@ where `c = C[tracer_index]`.
                                                           forcing) where tracer_index
 
     @inbounds c = tracers[tracer_index]
-    model_fields = merge(hydrostatic_fields(velocities, free_surface, tracers), auxiliary_fields)
+    model_fields = merge(hydrostatic_fields(velocities, free_surface, tracers), 
+                         auxiliary_fields,
+                         biogeochemical_auxiliary_fields(biogeochemistry))
 
     biogeochemical_velocities = biogeochemical_drift_velocity(biogeochemistry, val_tracer_name)
     closure_velocities = closure_turbulent_velocity(closure, diffusivities, val_tracer_name)
diff --git a/src/Models/HydrostaticFreeSurfaceModels/nothing_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/nothing_free_surface.jl
index 79a908a74e..66353ee04b 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/nothing_free_surface.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/nothing_free_surface.jl
@@ -2,5 +2,8 @@
 ##### There's no free surface with a rigid lid
 #####
 
+@inline materialize_free_surface(::Nothing, velocities, grid) = nothing
+
 @inline explicit_barotropic_pressure_x_gradient(i, j, k, grid, ::Nothing) = zero(grid)
 @inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, ::Nothing) = zero(grid)
+
diff --git a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl
index 1e44b9bea1..ac36afd25a 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl
@@ -100,6 +100,7 @@ extract_boundary_conditions(::PrescribedVelocityFields) = NamedTuple()
 free_surface_displacement_field(::PrescribedVelocityFields, ::Nothing, grid) = nothing
 HorizontalVelocityFields(::PrescribedVelocityFields, grid) = nothing, nothing
 
+materialize_free_surface(::Nothing,                      ::PrescribedVelocityFields, grid) = nothing
 materialize_free_surface(::ExplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, grid) = nothing
 materialize_free_surface(::ImplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, grid) = nothing
 materialize_free_surface(::SplitExplicitFreeSurface,     ::PrescribedVelocityFields, grid) = nothing
diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl
index ae986f0ec2..0eb1f422e9 100644
--- a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl
+++ b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl
@@ -194,14 +194,14 @@ function NonhydrostaticModel(; grid,
                                          extract_boundary_conditions(diffusivity_fields))
 
     # Next, we form a list of default boundary conditions:
-    prognostic_field_names = (:u, :v, :w, tracernames(tracers)..., keys(auxiliary_fields)...)
-    default_boundary_conditions = NamedTuple{prognostic_field_names}(FieldBoundaryConditions()
-                                                                     for name in prognostic_field_names)
+    field_names = (:u, :v, :w, tracernames(tracers)..., keys(auxiliary_fields)...)
+    default_boundary_conditions = NamedTuple{field_names}(FieldBoundaryConditions()
+                                                          for name in field_names)
 
     # Finally, we merge specified, embedded, and default boundary conditions. Specified boundary conditions
     # have precedence, followed by embedded, followed by default.
     boundary_conditions = merge(default_boundary_conditions, embedded_boundary_conditions, boundary_conditions)
-    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names)
+    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, field_names)
 
     # Ensure `closure` describes all tracers
     closure = with_tracers(tracernames(tracers), closure)
diff --git a/test/test_biogeochemistry.jl b/test/test_biogeochemistry.jl
index 46b92b8edc..18018161fa 100644
--- a/test/test_biogeochemistry.jl
+++ b/test/test_biogeochemistry.jl
@@ -69,7 +69,7 @@ const MB = Union{MinimalDiscreteBiogeochemistry, MinimalContinuousBiogeochemistr
 
 @inline          required_biogeochemical_tracers(::MB) = tuple(:P)
 @inline required_biogeochemical_auxiliary_fields(::MB) = tuple(:Iᴾᴬᴿ)
-@inline      biogeochemical_auxiliary_fields(bgc::MB) = (; Iᴾᴬᴿ = bgc.photosynthetic_active_radiation)
+@inline       biogeochemical_auxiliary_fields(bgc::MB) = (; Iᴾᴬᴿ = bgc.photosynthetic_active_radiation)
 @inline   biogeochemical_drift_velocity(bgc::MB, ::Val{:P}) = bgc.sinking_velocity
 
 # Update state test (won't actually change between calls but here to check it gets called)
diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl
index de169e0a2d..63e639b7d4 100644
--- a/test/test_boundary_conditions_integration.jl
+++ b/test/test_boundary_conditions_integration.jl
@@ -3,14 +3,14 @@ include("dependencies_for_runtests.jl")
 using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction, FlatExtrapolationOpenBoundaryCondition, fill_halo_regions!
 using Oceananigans: prognostic_fields
 
-function test_boundary_condition(arch, FT, topo, side, field_name, boundary_condition)
+function test_boundary_condition(arch, FT, Model, topo, side, field_name, boundary_condition)
     grid = RectilinearGrid(arch, FT, size=(1, 1, 1), extent=(1, π, 42), topology=topo)
 
     boundary_condition_kwarg = (; side => boundary_condition)
     field_boundary_conditions = FieldBoundaryConditions(; boundary_condition_kwarg...)
     bcs = (; field_name => field_boundary_conditions)
-    model = NonhydrostaticModel(; grid, boundary_conditions=bcs,
-                                buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
+    model = Model(; grid, boundary_conditions=bcs,
+                    buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
 
     success = try
         time_step!(model, 1e-16)
@@ -257,15 +257,15 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType),
             topo = (Bounded, Bounded, Bounded)
 
             for C in (Gradient, Flux, Value), boundary_condition in test_boundary_conditions(C, FT, array_type(arch))
-                @test test_boundary_condition(arch, FT, topo, :east, :T, boundary_condition)
-                @test test_boundary_condition(arch, FT, topo, :south, :T, boundary_condition)
-                @test test_boundary_condition(arch, FT, topo, :top, :T, boundary_condition)
+                @test test_boundary_condition(arch, FT, NonhydrostaticModel, topo, :east, :T, boundary_condition)
+                @test test_boundary_condition(arch, FT, NonhydrostaticModel, topo, :south, :T, boundary_condition)
+                @test test_boundary_condition(arch, FT, NonhydrostaticModel, topo, :top, :T, boundary_condition)
             end
 
             for boundary_condition in test_boundary_conditions(Open, FT, array_type(arch))
-                @test test_boundary_condition(arch, FT, topo, :east, :u, boundary_condition)
-                @test test_boundary_condition(arch, FT, topo, :south, :v, boundary_condition)
-                @test test_boundary_condition(arch, FT, topo, :top, :w, boundary_condition)
+                @test test_boundary_condition(arch, FT, NonhydrostaticModel, topo, :east, :u, boundary_condition)
+                @test test_boundary_condition(arch, FT, NonhydrostaticModel, topo, :south, :v, boundary_condition)
+                @test test_boundary_condition(arch, FT, NonhydrostaticModel, topo, :top, :w, boundary_condition)
             end
         end
     end
diff --git a/test/test_hydrostatic_free_surface_models.jl b/test/test_hydrostatic_free_surface_models.jl
index bffc166703..dfd703fbe0 100644
--- a/test/test_hydrostatic_free_surface_models.jl
+++ b/test/test_hydrostatic_free_surface_models.jl
@@ -102,7 +102,6 @@ topos_3d = ((Periodic, Periodic, Bounded),
             # SingleColumnGrid tests
             @test grid isa SingleColumnGrid
             @test isnothing(model.free_surface)
-            @test !(:η ∈ keys(fields(model))) # doesn't include free surface
         end
     end
 
@@ -113,7 +112,6 @@ topos_3d = ((Periodic, Periodic, Bounded),
                 grid = RectilinearGrid(arch, FT, topology=topo, size=(1, 1), extent=(1, 2))
                 model = HydrostaticFreeSurfaceModel(; grid)
                 @test model isa HydrostaticFreeSurfaceModel
-                @test :η ∈ keys(fields(model)) # contrary to the SingleColumnGrid case
             end
         end
     end
@@ -129,6 +127,17 @@ topos_3d = ((Periodic, Periodic, Bounded),
         end
     end
 
+    for FreeSurface in (ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface, Nothing)
+        @testset "$FreeSurface model construction" begin
+            @info "  Testing $FreeSurface model construction..."
+            for arch in archs, FT in float_types
+                grid = RectilinearGrid(arch, FT, size=(1, 1, 1), extent=(1, 2, 3))
+                model = HydrostaticFreeSurfaceModel(; grid, free_surface=FreeSurface())
+                @test model isa HydrostaticFreeSurfaceModel
+            end
+        end
+    end
+
     @testset "Halo size check in model constructor" begin
         for topo in topos_3d
             grid = RectilinearGrid(topology=topo, size=(1, 1, 1), extent=(1, 2, 3), halo=(1, 1, 1))
diff --git a/test/utils_for_runtests.jl b/test/utils_for_runtests.jl
index 11dc6ecd5e..7edcaeecd1 100644
--- a/test/utils_for_runtests.jl
+++ b/test/utils_for_runtests.jl
@@ -205,12 +205,11 @@ end
 ##### Boundary condition utils
 #####
 
-discrete_func(i, j, grid, clock, model_fields) = - model_fields.u[i, j, grid.Nz]
-parameterized_discrete_func(i, j, grid, clock, model_fields, p) = - p.μ * model_fields.u[i, j, grid.Nz]
-
-parameterized_fun(ξ, η, t, p) = p.μ * cos(p.ω * t)
-field_dependent_fun(ξ, η, t, u, v, w) = - w * sqrt(u^2 + v^2)
-exploding_fun(ξ, η, t, T, S, p) = - p.μ * cosh(S - p.S0) * exp((T - p.T0) / p.λ)
+@inline parameterized_discrete_func(i, j, grid, clock, model_fields, p) = - p.μ * model_fields.u[i, j, grid.Nz]
+@inline discrete_func(i, j, grid, clock, model_fields) = - model_fields.u[i, j, grid.Nz]
+@inline parameterized_fun(ξ, η, t, p) = p.μ * cos(p.ω * t)
+@inline field_dependent_fun(ξ, η, t, u, v, w) = - w * sqrt(u^2 + v^2)
+@inline exploding_fun(ξ, η, t, T, S, p) = - p.μ * cosh(S - p.S0) * exp((T - p.T0) / p.λ)
 
 # Many bc. Very many
                  integer_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, 1)