diff --git a/src/shared_utilities/utils.jl b/src/shared_utilities/utils.jl
index 5629f3c660..7a3d4fb339 100644
--- a/src/shared_utilities/utils.jl
+++ b/src/shared_utilities/utils.jl
@@ -6,7 +6,7 @@ import SciMLBase
 Computes the heaviside function.
 """
 function heaviside(x::FT)::FT where {FT}
-    if x >= FT(0.0)
+    if x > eps(FT)
         return FT(1.0)
     else
         return FT(0.0)
diff --git a/test/standalone/Bucket/snow_bucket_tests.jl b/test/standalone/Bucket/snow_bucket_tests.jl
index e6c733d9e6..a94d92d478 100644
--- a/test/standalone/Bucket/snow_bucket_tests.jl
+++ b/test/standalone/Bucket/snow_bucket_tests.jl
@@ -59,9 +59,8 @@ bucket_domains = [
 ]
 init_temp(z::FT, value::FT) where {FT} = FT(value)
 orbital_data = Insolation.OrbitalData()
-for bucket_domain in bucket_domains
-
-    @testset "Conservation of water and energy" begin
+for i in 1:3
+    @testset "Conservation of water and energy I (snow present)" begin
         "Radiation"
         ref_time = DateTime(2005)
         SW_d = (t) -> eltype(t)(20.0)
@@ -69,15 +68,17 @@ for bucket_domain in bucket_domains
         bucket_rad =
             PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; orbital_data)
         "Atmos"
-        precip = (t) -> eltype(t)(0) # no precipitation
+        liquid_precip = (t) -> eltype(t)(1e-8) # precipitation
+        snow_precip = (t) -> eltype(t)(1e-7) # precipitation
+
         T_atmos = (t) -> eltype(t)(280.0)
         u_atmos = (t) -> eltype(t)(10.0)
         q_atmos = (t) -> eltype(t)(0.03)
         h_atmos = FT(3)
         P_atmos = (t) -> eltype(t)(101325) # Pa
         bucket_atmos = PrescribedAtmosphere(
-            precip,
-            precip,
+            liquid_precip,
+            snow_precip,
             T_atmos,
             u_atmos,
             q_atmos,
@@ -100,7 +101,7 @@ for bucket_domain in bucket_domains
         )
         model = BucketModel(
             parameters = bucket_parameters,
-            domain = bucket_domain,
+            domain = bucket_domains[i],
             atmosphere = bucket_atmos,
             radiation = bucket_rad,
         )
@@ -118,7 +119,114 @@ for bucket_domain in bucket_domains
         set_initial_aux_state! = make_set_initial_aux_state(model)
         set_initial_aux_state!(p, Y, t0)
         exp_tendency!(dY, Y, p, t0)
+        _LH_f0 = LSMP.LH_f0(model.parameters.earth_param_set)
+        _ρ_liq = LSMP.ρ_cloud_liq(model.parameters.earth_param_set)
+        _ρLH_f0 = _ρ_liq * _LH_f0 # Latent heat per unit volume
+        _T_freeze = LSMP.T_freeze(model.parameters.earth_param_set)
+        snow_cover_fraction(σS) = σS > eps(FT) ? FT(1.0) : FT(0.0)
+
+        partitioned_fluxes =
+            partition_surface_fluxes.(
+                Y.bucket.σS,
+                p.bucket.T_sfc,
+                model.parameters.τc,
+                snow_cover_fraction.(Y.bucket.σS),
+                p.bucket.evaporation,
+                p.bucket.turbulent_energy_flux .+ p.bucket.R_n,
+                _ρLH_f0,
+                _T_freeze,
+            )
+        F_melt = partitioned_fluxes.F_melt
+        F_into_snow =
+            partitioned_fluxes.F_into_snow .+ _ρLH_f0 .* snow_precip(t0)
+        G = partitioned_fluxes.G
+        F_sfc =
+            p.bucket.turbulent_energy_flux .+ p.bucket.R_n .+
+            _ρLH_f0 .* snow_precip(t0)
+        F_water_sfc =
+            liquid_precip(t0) + snow_precip(t0) .- p.bucket.evaporation
+
+
+        if i == 1
+            A_point = sum(ones(bucket_domains[i].surface.space))
+        else
+            A_point = 1
+        end
+
+        dIsnow = -_ρLH_f0 .* dY.bucket.σS
+        @test sum(dIsnow) / A_point ≈ sum(-1 .* F_into_snow) / A_point
+
+        de_soil = dY.bucket.T .* ρc_soil
+        @test sum(de_soil) ≈ sum(-1 .* G) / A_point
+
+        dWL = dY.bucket.W .+ dY.bucket.Ws .+ dY.bucket.σS
+        @test sum(dWL) / A_point ≈ sum(F_water_sfc) / A_point
+
+        dIL = sum(dIsnow) / A_point .+ sum(de_soil)
+        @test dIL ≈ sum(-1 .* F_sfc) / A_point
+    end
+end
+
+for i in 1:3
+    @testset "Conservation of water and energy II (no snow to start)" begin
+        "Radiation"
+        ref_time = DateTime(2005)
+        SW_d = (t) -> eltype(t)(20.0)
+        LW_d = (t) -> eltype(t)(20.0)
+        bucket_rad =
+            PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; orbital_data)
+        "Atmos"
+        liquid_precip = (t) -> eltype(t)(1e-8) # precipitation
+        snow_precip = (t) -> eltype(t)(1e-7) # precipitation
+
+        T_atmos = (t) -> eltype(t)(280.0)
+        u_atmos = (t) -> eltype(t)(10.0)
+        q_atmos = (t) -> eltype(t)(0.03)
+        h_atmos = FT(3)
+        P_atmos = (t) -> eltype(t)(101325) # Pa
+        bucket_atmos = PrescribedAtmosphere(
+            liquid_precip,
+            snow_precip,
+            T_atmos,
+            u_atmos,
+            q_atmos,
+            P_atmos,
+            ref_time,
+            h_atmos,
+        )
+        Δt = FT(1.0)
+        τc = FT(10.0)
+        bucket_parameters = BucketModelParameters(
+            κ_soil,
+            ρc_soil,
+            albedo,
+            σS_c,
+            W_f,
+            z_0m,
+            z_0b,
+            τc,
+            earth_param_set,
+        )
+        model = BucketModel(
+            parameters = bucket_parameters,
+            domain = bucket_domains[i],
+            atmosphere = bucket_atmos,
+            radiation = bucket_rad,
+        )
 
+        # run for some time
+        Y, p, coords = initialize(model)
+        Y.bucket.T .= init_temp.(coords.subsurface.z, 274.0)
+        Y.bucket.W .= 0.0 # no moisture
+        Y.bucket.Ws .= 0.0 # no runoff
+        Y.bucket.σS .= 0.0
+        t0 = FT(0.0)
+        dY = similar(Y)
+
+        exp_tendency! = make_exp_tendency(model)
+        set_initial_aux_state! = make_set_initial_aux_state(model)
+        set_initial_aux_state!(p, Y, t0)
+        exp_tendency!(dY, Y, p, t0)
         _LH_f0 = LSMP.LH_f0(model.parameters.earth_param_set)
         _ρ_liq = LSMP.ρ_cloud_liq(model.parameters.earth_param_set)
         _ρLH_f0 = _ρ_liq * _LH_f0 # Latent heat per unit volume
@@ -137,21 +245,134 @@ for bucket_domain in bucket_domains
                 _T_freeze,
             )
         F_melt = partitioned_fluxes.F_melt
-        F_into_snow = partitioned_fluxes.F_into_snow
+        F_into_snow =
+            partitioned_fluxes.F_into_snow .+ _ρLH_f0 .* snow_precip(t0)
         G = partitioned_fluxes.G
-        F_sfc = p.bucket.turbulent_energy_flux .+ p.bucket.R_n
-        F_water_sfc = precip(t0) .- p.bucket.evaporation
+        F_sfc =
+            p.bucket.turbulent_energy_flux .+ p.bucket.R_n .+
+            _ρLH_f0 .* snow_precip(t0)
+        F_water_sfc =
+            liquid_precip(t0) + snow_precip(t0) .- p.bucket.evaporation
+
+        if i == 1
+            A_point = sum(ones(bucket_domains[i].surface.space))
+        else
+            A_point = 1
+        end
 
         dIsnow = -_ρLH_f0 .* dY.bucket.σS
-        @test sum(dIsnow) ≈ sum(-1 .* F_into_snow)
+        @test sum(dIsnow) / A_point ≈ sum(-1 .* F_into_snow) / A_point
 
         de_soil = dY.bucket.T .* ρc_soil
-        @test sum(de_soil) ≈ sum(-1 .* G)
+        @test sum(de_soil) ≈ sum(-1 .* G) / A_point
 
         dWL = dY.bucket.W .+ dY.bucket.Ws .+ dY.bucket.σS
-        @test sum(dWL) ≈ sum(F_water_sfc)
+        @test sum(dWL) / A_point ≈ sum(F_water_sfc) / A_point
 
-        dIL = sum(dIsnow) .+ sum(de_soil)
-        @test dIL ≈ sum(-1 .* F_sfc)
+        dIL = sum(dIsnow) / A_point .+ sum(de_soil)
+        @test dIL ≈ sum(-1 .* F_sfc) / A_point
     end
 end
+
+
+
+@testset "Conservation of water and energy - nonuniform evaporation" begin
+    i = 3
+    "Radiation"
+    ref_time = DateTime(2005)
+    SW_d = (t) -> eltype(t)(20.0)
+    LW_d = (t) -> eltype(t)(20.0)
+    bucket_rad =
+        PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; orbital_data)
+    "Atmos"
+    liquid_precip = (t) -> eltype(t)(1e-8) # precipitation
+    snow_precip = (t) -> eltype(t)(1e-7) # precipitation
+
+    T_atmos = (t) -> eltype(t)(280.0)
+    u_atmos = (t) -> eltype(t)(10.0)
+    q_atmos = (t) -> eltype(t)(0.03)
+    h_atmos = FT(3)
+    P_atmos = (t) -> eltype(t)(101325) # Pa
+    bucket_atmos = PrescribedAtmosphere(
+        liquid_precip,
+        snow_precip,
+        T_atmos,
+        u_atmos,
+        q_atmos,
+        P_atmos,
+        ref_time,
+        h_atmos,
+    )
+    Δt = FT(1.0)
+    τc = FT(10.0)
+    bucket_parameters = BucketModelParameters(
+        κ_soil,
+        ρc_soil,
+        albedo,
+        σS_c,
+        W_f,
+        z_0m,
+        z_0b,
+        τc,
+        earth_param_set,
+    )
+    model = BucketModel(
+        parameters = bucket_parameters,
+        domain = bucket_domains[i],
+        atmosphere = bucket_atmos,
+        radiation = bucket_rad,
+    )
+
+    # run for some time
+    Y, p, coords = initialize(model)
+    Y.bucket.T .= init_temp.(coords.subsurface.z, 274.0)
+    Y.bucket.W .= 0.0 # no moisture
+    Y.bucket.Ws .= 0.0 # no runoff
+    Y.bucket.σS .= 0.0
+    t0 = FT(0.0)
+    dY = similar(Y)
+
+    compute_exp_tendency! = ClimaLSM.make_compute_exp_tendency(model)
+    set_initial_aux_state! = make_set_initial_aux_state(model)
+    set_initial_aux_state!(p, Y, t0)
+    random = zeros(bucket_domains[i].surface.space)
+    parent(random) .= rand(FT, size(parent(random)))
+    p.bucket.evaporation .= random .* 1e-7
+    compute_exp_tendency!(dY, Y, p, t0)
+    _LH_f0 = LSMP.LH_f0(model.parameters.earth_param_set)
+    _ρ_liq = LSMP.ρ_cloud_liq(model.parameters.earth_param_set)
+    _ρLH_f0 = _ρ_liq * _LH_f0 # Latent heat per unit volume
+    _T_freeze = LSMP.T_freeze(model.parameters.earth_param_set)
+    snow_cover_fraction(σS) = σS > eps(FT) ? FT(1.0) : FT(0.0)
+
+    partitioned_fluxes =
+        partition_surface_fluxes.(
+            Y.bucket.σS,
+            p.bucket.T_sfc,
+            model.parameters.τc,
+            snow_cover_fraction.(Y.bucket.σS),
+            p.bucket.evaporation,
+            p.bucket.turbulent_energy_flux .+ p.bucket.R_n,
+            _ρLH_f0,
+            _T_freeze,
+        )
+    F_melt = partitioned_fluxes.F_melt
+    F_into_snow = partitioned_fluxes.F_into_snow .+ _ρLH_f0 .* snow_precip(t0)
+    G = partitioned_fluxes.G
+    F_sfc =
+        p.bucket.turbulent_energy_flux .+ p.bucket.R_n .+
+        _ρLH_f0 .* snow_precip(t0)
+    F_water_sfc = liquid_precip(t0) + snow_precip(t0) .- p.bucket.evaporation
+
+    dIsnow = -_ρLH_f0 .* dY.bucket.σS
+    @test sum(dIsnow) ≈ sum(-1 .* F_into_snow)
+
+    de_soil = dY.bucket.T .* ρc_soil
+    @test sum(de_soil) ≈ sum(-1 .* G)
+
+    dWL = dY.bucket.W .+ dY.bucket.Ws .+ dY.bucket.σS
+    @test sum(dWL) ≈ sum(F_water_sfc)
+
+    dIL = sum(dIsnow) .+ sum(de_soil)
+    @test dIL ≈ sum(-1 .* F_sfc)
+end