diff --git a/test/Project.toml b/test/Project.toml
index f25faed39f..faabbafd01 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -12,6 +12,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f"
@@ -21,3 +22,4 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 [compat]
 CLIMAParameters = "0.7"
 NLsolve = "4.5"
+SafeTestsets = "0.1"
diff --git a/test/runtests.jl b/test/runtests.jl
index db3fd660fb..fd541b0356 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,27 +1,98 @@
-include("aqua.jl")
-include("integrated/lsms.jl")
-include("standalone/Snow/parameterizations.jl")
-include("shared_utilities/implicit_timestepping/richards_model.jl")
-include("standalone/Vegetation/test_bigleaf_parameterizations.jl")
-include("standalone/Vegetation/canopy_model.jl")
-include("standalone/Vegetation/test_two_stream.jl")
-include("standalone/Soil/Biogeochemistry/co2_parameterizations.jl")
-include("standalone/Soil/Biogeochemistry/biogeochemistry_module.jl")
-include("standalone/Bucket/soil_bucket_tests.jl")
-include("standalone/Bucket/snow_bucket_tests.jl")
-include("standalone/Bucket/albedo_types.jl")
-include("standalone/Vegetation/plant_hydraulics_test.jl")
-include("shared_utilities/domains.jl")
-include("shared_utilities/file_reader.jl")
-include("shared_utilities/regridder.jl")
-include("shared_utilities/utilities.jl")
-include("shared_utilities/variable_types.jl")
-include("standalone/SurfaceWater/pond_test.jl")
-include("integrated/pond_soil_lsm.jl")
-include("integrated/soil_energy_hydrology_biogeochemistry.jl")
-include("standalone/Soil/climate_drivers.jl")
-include("standalone/Soil/soil_test_3d.jl")
-include("standalone/Soil/soiltest.jl")
-include("standalone/Soil/soil_parameterizations.jl")
-include("standalone/Soil/soil_bc.jl")
-include("standalone/Soil/runoff.jl")
+using SafeTestsets
+
+# Performance and code quality tests
+@safetestset "Aqua tests" begin
+    include("aqua.jl")
+end
+
+# Integrated LSM tests
+@safetestset "Integrated LSM unit tests" begin
+    include("integrated/lsms.jl")
+end
+@safetestset "Integrated pond/soil LSM tests" begin
+    include("integrated/pond_soil_lsm.jl")
+end
+@safetestset "Integrated soil energy/hydrology/biogeochem LSM tests" begin
+    include("integrated/pond_soil_lsm.jl")
+end
+
+# Shared ClimaLSM utilities tests
+@safetestset "Richards model implicit timestepping tests" begin
+    include("shared_utilities/implicit_timestepping/richards_model.jl")
+end
+@safetestset "Domains module tests" begin
+    include("shared_utilities/domains.jl")
+end
+@safetestset "FileReader module tests" begin
+    include("shared_utilities/file_reader.jl")
+end
+@safetestset "Regridder module tests" begin
+    include("shared_utilities/regridder.jl")
+end
+@safetestset "General utilities tests" begin
+    include("shared_utilities/utilities.jl")
+end
+@safetestset "Variable types tests" begin
+    include("shared_utilities/variable_types.jl")
+end
+
+# Standalone Bucket model tests
+@safetestset "Bucket albedo types tests" begin
+    include("standalone/Bucket/albedo_types.jl")
+end
+@safetestset "Bucket snow tests" begin
+    include("standalone/Bucket/snow_bucket_tests.jl")
+end
+@safetestset "Bucket soil tests" begin
+    include("standalone/Bucket/soil_bucket_tests.jl")
+end
+
+# Standalone Snow model tests
+@safetestset "Snow parameterization tests" begin
+    include("standalone/Snow/parameterizations.jl")
+end
+
+# Standalone Soil model tests
+@safetestset "Soil Biogeochemistry module tests" begin
+    include("standalone/Soil/Biogeochemistry/biogeochemistry_module.jl")
+end
+@safetestset "Soil CO2 parameterization tests" begin
+    include("standalone/Soil/Biogeochemistry/co2_parameterizations.jl")
+end
+@safetestset "Soil climate drivers tests" begin
+    include("standalone/Soil/climate_drivers.jl")
+end
+@safetestset "Soil runoff tests" begin
+    include("standalone/Soil/runoff.jl")
+end
+@safetestset "Soil boundary condition tests" begin
+    include("standalone/Soil/soil_bc.jl")
+end
+@safetestset "Soil parameterization tests" begin
+    include("standalone/Soil/soil_parameterizations.jl")
+end
+@safetestset "Soil 3D domain tests" begin
+    include("standalone/Soil/soil_test_3d.jl")
+end
+@safetestset "Soil integration tests" begin
+    include("standalone/Soil/soiltest.jl")
+end
+
+# Standalone Surface Water model tests
+@safetestset "Pond module tests" begin
+    include("standalone/SurfaceWater/pond_test.jl")
+end
+
+# Standalone Vegetation model tests
+@safetestset "Canopy module tests" begin
+    include("standalone/Vegetation/canopy_model.jl")
+end
+@safetestset "Canopy PlantHydraulics tests" begin
+    include("standalone/Vegetation/plant_hydraulics_test.jl")
+end
+@safetestset "Big Leaf model tests" begin
+    include("standalone/Vegetation/test_bigleaf_parameterizations.jl")
+end
+@safetestset "Two Stream model tests" begin
+    include("standalone/Vegetation/test_two_stream.jl")
+end
diff --git a/test/standalone/Bucket/snow_bucket_tests.jl b/test/standalone/Bucket/snow_bucket_tests.jl
index 6179e837e9..5233b84788 100644
--- a/test/standalone/Bucket/snow_bucket_tests.jl
+++ b/test/standalone/Bucket/snow_bucket_tests.jl
@@ -1,6 +1,7 @@
 using Test
 
 using Statistics
+using Insolation
 using ClimaCore
 using ClimaLSM.Bucket:
     BucketModel,
diff --git a/test/standalone/Bucket/soil_bucket_tests.jl b/test/standalone/Bucket/soil_bucket_tests.jl
index 88b566e5b9..289c92f34e 100644
--- a/test/standalone/Bucket/soil_bucket_tests.jl
+++ b/test/standalone/Bucket/soil_bucket_tests.jl
@@ -1,6 +1,7 @@
 using Test
 
 using Statistics
+using Insolation
 using ClimaCore
 import CLIMAParameters as CP
 using ClimaLSM.Bucket: BucketModel, BucketModelParameters, BulkAlbedoFunction
diff --git a/test/standalone/Soil/soil_bc.jl b/test/standalone/Soil/soil_bc.jl
index 2d2c898845..b19629894e 100644
--- a/test/standalone/Soil/soil_bc.jl
+++ b/test/standalone/Soil/soil_bc.jl
@@ -1,8 +1,11 @@
 using Test
+using Statistics
 using ClimaCore
 using ClimaLSM
+using ClimaLSM.Soil
+using ClimaLSM.Domains: HybridBox, SphericalShell, Column
 
-using ClimaLSM.Domains: HybridBox, SphericalShell
+include(joinpath(pkgdir(ClimaLSM), "parameters", "create_parameters.jl"))
 FT = Float64
 
 @testset "WVector usage in gradient" begin
@@ -89,3 +92,109 @@ FT = Float64
 
     #div_rpt = unique(parent(@. (divf2c(gradc2f(scalar)))))
 end
+
+@testset "Test RRE state to flux BC calculations" begin
+    ν = FT(0.495)
+    K_sat = FT(0.0443 / 3600 / 100) # m/s
+    S_s = FT(1e-3) #inverse meters
+    vg_n = FT(2.0)
+    vg_α = FT(2.6) # inverse meters
+    hcm = vanGenuchten(; α = vg_α, n = vg_n)
+    θ_r = FT(0)
+    zmax = FT(0)
+    zmin = FT(-10)
+    nelems = 50
+    soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
+    z = ClimaCore.Fields.coordinate_field(soil_domain.space).z
+    Δz = abs(zmax - zmin) / nelems / 2.0
+
+    top_Δz, bottom_Δz = get_Δz(z)
+    @test (mean(abs.(parent(top_Δz) .- Δz)) < 1e-13) &&
+          (mean(abs.(parent(bottom_Δz) .- Δz)) < 1e-13)
+
+    ϑ_bc = ν / 2
+    ϑ_c = ν / 3
+
+    K_c = hydraulic_conductivity(
+        hcm,
+        K_sat,
+        Soil.effective_saturation(ν, ϑ_c, θ_r),
+    )
+
+    ψ_bc = pressure_head(hcm, θ_r, ϑ_bc, ν, S_s)
+    ψ_c = pressure_head(hcm, θ_r, ϑ_c, ν, S_s)
+
+    flux_int = diffusive_flux(K_c, ψ_bc + Δz, ψ_c, Δz)
+    flux_expected = -K_c * ((ψ_bc - ψ_c + Δz) / Δz)
+
+    @test abs(flux_expected - flux_int) < 1e-13
+end
+
+@testset "Test heat state to flux BC calculations" begin
+    earth_param_set = create_lsm_parameters(FT)
+    ν = FT(0.495)
+    K_sat = FT(0.0443 / 3600 / 100) # m/s
+    S_s = FT(1e-3) #inverse meters
+    vg_n = FT(2.0)
+    vg_α = FT(2.6) # inverse meters
+    hcm = vanGenuchten(; α = vg_α, n = vg_n)
+    θ_r = FT(0.1)
+    ν_ss_om = FT(0.0)
+    ν_ss_quartz = FT(1.0)
+    ν_ss_gravel = FT(0.0)
+    κ_minerals = FT(2.5)
+    κ_om = FT(0.25)
+    κ_quartz = FT(8.0)
+    κ_air = FT(0.025)
+    κ_ice = FT(2.21)
+    κ_liq = FT(0.57)
+    ρp = FT(2.66 / 1e3 * 1e6)
+    ρc_ds = FT(2e6 * (1.0 - ν))
+    κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals)
+    κ_dry = Soil.κ_dry(ρp, ν, κ_solid, κ_air)
+    κ_sat_frozen = Soil.κ_sat_frozen(κ_solid, ν, κ_ice)
+    κ_sat_unfrozen = Soil.κ_sat_unfrozen(κ_solid, ν, κ_liq)
+
+    hyd_on_en_on = Soil.EnergyHydrologyParameters{FT}(;
+        κ_dry = κ_dry,
+        κ_sat_frozen = κ_sat_frozen,
+        κ_sat_unfrozen = κ_sat_unfrozen,
+        ρc_ds = ρc_ds,
+        ν = ν,
+        ν_ss_om = ν_ss_om,
+        ν_ss_quartz = ν_ss_quartz,
+        ν_ss_gravel = ν_ss_gravel,
+        hydrology_cm = hcm,
+        K_sat = K_sat,
+        S_s = S_s,
+        θ_r = θ_r,
+        earth_param_set = earth_param_set,
+    )
+
+    zmax = FT(0)
+    zmin = FT(-1)
+    nelems = 200
+    Δz = abs(zmax - zmin) / nelems
+
+    ϑ_bc = ν / 2
+    ϑ_c = ν / 3
+    θ_i = FT(0)
+    T_bc = 298
+    T_c = 290
+
+    κ_c =
+        thermal_conductivity.(
+            κ_dry,
+            kersten_number.(
+                θ_i,
+                relative_saturation.(ϑ_c, θ_i, ν),
+                Ref(hyd_on_en_on),
+            ),
+            κ_sat.(ϑ_c, θ_i, κ_sat_unfrozen, κ_sat_frozen),
+        )
+
+    flux_int = diffusive_flux(κ_c, T_bc, T_c, Δz)
+    flux_expected = -κ_c * (T_bc - T_c) / Δz
+
+    @test abs(flux_expected - flux_int) < 1e-13
+end
diff --git a/test/standalone/Soil/soiltest.jl b/test/standalone/Soil/soiltest.jl
index c0418d94c3..d82e77221a 100644
--- a/test/standalone/Soil/soiltest.jl
+++ b/test/standalone/Soil/soiltest.jl
@@ -13,112 +13,6 @@ include(joinpath(pkgdir(ClimaLSM), "parameters", "create_parameters.jl"))
 
 FT = Float64
 
-@testset "Test RRE state to flux BC calculations" begin
-    ν = FT(0.495)
-    K_sat = FT(0.0443 / 3600 / 100) # m/s
-    S_s = FT(1e-3) #inverse meters
-    vg_n = FT(2.0)
-    vg_α = FT(2.6) # inverse meters
-    hcm = vanGenuchten(; α = vg_α, n = vg_n)
-    θ_r = FT(0)
-    zmax = FT(0)
-    zmin = FT(-10)
-    nelems = 50
-    soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems)
-    z = ClimaCore.Fields.coordinate_field(soil_domain.space).z
-    Δz = abs(zmax - zmin) / nelems / 2.0
-
-    top_Δz, bottom_Δz = get_Δz(z)
-    @test (mean(abs.(parent(top_Δz) .- Δz)) < 1e-13) &&
-          (mean(abs.(parent(bottom_Δz) .- Δz)) < 1e-13)
-
-    ϑ_bc = ν / 2
-    ϑ_c = ν / 3
-
-    K_c = hydraulic_conductivity(
-        hcm,
-        K_sat,
-        Soil.effective_saturation(ν, ϑ_c, θ_r),
-    )
-
-    ψ_bc = pressure_head(hcm, θ_r, ϑ_bc, ν, S_s)
-    ψ_c = pressure_head(hcm, θ_r, ϑ_c, ν, S_s)
-
-    flux_int = diffusive_flux(K_c, ψ_bc + Δz, ψ_c, Δz)
-    flux_expected = -K_c * ((ψ_bc - ψ_c + Δz) / Δz)
-
-    @test abs(flux_expected - flux_int) < 1e-13
-end
-
-@testset "Test heat state to flux BC calculations" begin
-    earth_param_set = create_lsm_parameters(FT)
-    ν = FT(0.495)
-    K_sat = FT(0.0443 / 3600 / 100) # m/s
-    S_s = FT(1e-3) #inverse meters
-    vg_n = FT(2.0)
-    vg_α = FT(2.6) # inverse meters
-    hcm = vanGenuchten(; α = vg_α, n = vg_n)
-    θ_r = FT(0.1)
-    ν_ss_om = FT(0.0)
-    ν_ss_quartz = FT(1.0)
-    ν_ss_gravel = FT(0.0)
-    κ_minerals = FT(2.5)
-    κ_om = FT(0.25)
-    κ_quartz = FT(8.0)
-    κ_air = FT(0.025)
-    κ_ice = FT(2.21)
-    κ_liq = FT(0.57)
-    ρp = FT(2.66 / 1e3 * 1e6)
-    ρc_ds = FT(2e6 * (1.0 - ν))
-    κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals)
-    κ_dry = Soil.κ_dry(ρp, ν, κ_solid, κ_air)
-    κ_sat_frozen = Soil.κ_sat_frozen(κ_solid, ν, κ_ice)
-    κ_sat_unfrozen = Soil.κ_sat_unfrozen(κ_solid, ν, κ_liq)
-
-    hyd_on_en_on = Soil.EnergyHydrologyParameters{FT}(;
-        κ_dry = κ_dry,
-        κ_sat_frozen = κ_sat_frozen,
-        κ_sat_unfrozen = κ_sat_unfrozen,
-        ρc_ds = ρc_ds,
-        ν = ν,
-        ν_ss_om = ν_ss_om,
-        ν_ss_quartz = ν_ss_quartz,
-        ν_ss_gravel = ν_ss_gravel,
-        hydrology_cm = hcm,
-        K_sat = K_sat,
-        S_s = S_s,
-        θ_r = θ_r,
-        earth_param_set = earth_param_set,
-    )
-
-    zmax = FT(0)
-    zmin = FT(-1)
-    nelems = 200
-    Δz = abs(zmax - zmin) / nelems
-
-    ϑ_bc = ν / 2
-    ϑ_c = ν / 3
-    θ_i = FT(0)
-    T_bc = 298
-    T_c = 290
-
-    κ_c =
-        thermal_conductivity.(
-            κ_dry,
-            kersten_number.(
-                θ_i,
-                relative_saturation.(ϑ_c, θ_i, ν),
-                Ref(hyd_on_en_on),
-            ),
-            κ_sat.(ϑ_c, θ_i, κ_sat_unfrozen, κ_sat_frozen),
-        )
-
-    flux_int = diffusive_flux(κ_c, T_bc, T_c, Δz)
-    flux_expected = -κ_c * (T_bc - T_c) / Δz
-
-    @test abs(flux_expected - flux_int) < 1e-13
-end
-
 @testset "Richards equation with flux BCs" begin
     ν = FT(0.495)
     K_sat = FT(0.0443 / 3600 / 100) # m/s