Apply snow thermal properties by component

src/Physics/Snow/Snow.jl

@@ -24,17 +24,12 @@ abstract type SnowParameterization <: CryoGrid.Parameterization end
 
 Simple, bulk snow scheme where snowpack is represented as a single grid cell with homogenous state.
 """
-Base.@kwdef struct Bulk{Tdsn,Tthresh} <: SnowParameterization
+Base.@kwdef struct Bulk{Tdsn,Tρsn,Tthresh} <: SnowParameterization
     dsn::Tdsn = 0.0u"m" # snow depth
+    ρsn::Tρsn = 250.0u"kg/m^3" # bulk snow density
     thresh::Tthresh = 0.05u"m" # snow threshold
 end
 
-Base.@kwdef struct SnowThermalProperties{Tρ,Tk,Tc} <: HeatConduction.ThermalProperties
-    ρsn::Tρ = 350.0u"kg/m^3" # bulk snow density
-    ksn::Tk = 0.5u"W/m/K" # thermal conductivity of fresh snow (Westermann et al. 2011)
-    csn::Tc = 7.5e6u"J/K/m^3" # heat capacity of snow (Westermann et al. 2011)
-end
-
 struct SnowFree <: CryoGrid.Callback end
 CryoGrid.CallbackStyle(::Type{SnowFree}) = CryoGrid.Continuous()
 
@@ -45,7 +40,6 @@ Generic representation of a ground surface snow pack.
 """
 Base.@kwdef struct Snowpack{TPara<:SnowParameterization,TProp,TSp} <: CryoGrid.SubSurface
     para::TPara = Bulk()
-    prop::TProp = SnowThermalProperties()
     sp::TSp = nothing
 end
 

src/Physics/Snow/snow_bulk.jl

@@ -29,10 +29,7 @@ end
 
 snowdepth(snow::BulkSnowpack{<:Number}, state) = snow.para.dsn
 snowdepth(snow::BulkSnowpack{<:Forcing}, state) = snow.para.dsn(state.t)
-
-Physics.volumetricfractions(::BulkSnowpack, ::Heat, state, i) = (1.0,)
-HeatConduction.thermalconductivities(snow::Snowpack, ::Heat) = (snow.prop.ksn,)
-HeatConduction.heatcapacities(snow::Snowpack, ::Heat) = (snow.prop.csn,)
+snowdensity(snow::BulkSnowpack, state) = snow.para.ρsn
 
 # CryoGrid methods
 CryoGrid.thickness(::BulkSnowpack, state) = getscalar(state.dsn)
@@ -42,12 +39,14 @@ function CryoGrid.diagnosticstep!(snow::BulkSnowpack, state)
 end
 function CryoGrid.diagnosticstep!(snow::BulkSnowpack, heat::Heat{FreeWater,Enthalpy}, state)
     dsn = getscalar(state.dsn)
-    @. state.θwi = snow.prop.ρsn / heat.prop.ρw
+    ρsn = snowdensity(snow, state)
+    @. state.θwi = ρsn / heat.prop.ρw
     if dsn > snow.para.thresh
         @. state.θwi = snow.prop.ρsn / heat.prop.ρw
         HeatConduction.freezethaw!(snow, heat, state)
     else
         @. state.θwi = 0.0
+        @. state.θw = 0.0
         @. state.T = state.T_ub
         @. state.C = heat.prop.ca
     end