Initial implementation of bulk snow scheme

examples/heat_freeW_samoylov_snow.jl

+using CryoGrid
+using Plots
+
+forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA5_fitted_daily_1979_2020, :Tair => u"°C", :Dsn => u"m"; spec=JsonSpec{2});
+# use air temperature as upper boundary forcing;
+tair = TimeSeriesForcing(ustrip.(forcings.data.Tair), forcings.timestamps, :Tair);
+snowdepth = TimeSeriesForcing(ustrip.(forcings.data.Dsn), forcings.timestamps, :Dsn);
+# use default profiles for samoylov
+soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
+# "simple" heat conduction model w/ 5 cm grid spacing (defaults to free water freezing scheme)
+modelgrid = Grid(vcat([-1.0u"m"], CryoGrid.Presets.DefaultGrid_5cm))
+initT = initializer(:T, tempprofile)
+z_top = -2.0u"m"
+z_sub = map(knot -> knot.depth, soilprofile)
+z_bot = modelgrid[end]
+strat = @Stratigraphy(
+    z_top => top(TemperatureGradient(tair)),
+    z_top => subsurface(:snowpack, Snowpack(para=Snow.Bulk(dsn=snowdepth)), Heat(:H)),
+    z_sub[1] => subsurface(:topsoil1, Soil(para=soilprofile[1].value), Heat(:H)),
+    z_sub[2] => subsurface(:topsoil2, Soil(para=soilprofile[2].value), Heat(:H)),
+    z_sub[3] => subsurface(:sediment1, Soil(para=soilprofile[3].value), Heat(:H)),
+    z_sub[4] => subsurface(:sediment2, Soil(para=soilprofile[4].value), Heat(:H)),
+    z_sub[5] => subsurface(:sediment3, Soil(para=soilprofile[5].value), Heat(:H)),
+    z_bot => bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
+);
+tile = Tile(strat, modelgrid, initT)
+# define time span
+tspan = (DateTime(2010,10,30),DateTime(2012,12,30))
+p = parameters(tile)
+u0, du0 = initialcondition!(tile, tspan, p)
+# CryoGrid front-end for ODEProblem
+prob = CryoGridProblem(tile,u0,tspan,p,savevars=(:T,:C,:kc))
+# solve with forward Euler, 10-minute time steps
+out = @time solve(prob, Euler(), dt=600.0, saveat=24*3600.0, progress=true) |> CryoGridOutput;
+# Plot it!
+zs = [1,10,20,30,50,100,200,500,1000]u"cm"
+cg = Plots.cgrad(:copper,rev=true);
+plot(ustrip(out.T[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, dpi=150)
+plot(ustrip(out.H[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150)
+
+integrator = init(prob, Euler(), dt=300.0)
+step!(integrator)
+snow_state = getstate(:snowpack, integrator)
+snow_state.kc

src/Physics/HeatConduction/HeatConduction.jl

 module HeatConduction
 
 import CryoGrid: SubSurfaceProcess, BoundaryStyle, Dirichlet, Neumann, BoundaryProcess, Layer, Top, Bottom, SubSurface, Callback
-import CryoGrid: diagnosticstep!, prognosticstep!, interact!, initialcondition!, boundaryflux, boundaryvalue, variables, callbacks, criterion, affect!, thickness, midpoints
+import CryoGrid: diagnosticstep!, prognosticstep!, interact!, initialcondition!, boundaryflux, boundaryvalue, variables, callbacks, criterion, affect!, thickness, midpoints, volumetricfractions
 
 using CryoGrid.InputOutput: Forcing
 using CryoGrid.Physics
@@ -74,7 +74,7 @@ freezecurve(heat::Heat) = heat.freezecurve
 # Default implementation of `variables` for freeze curve
 variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
 
-export ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor, Damping
+export HeatBC, ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor
 include("heat_bc.jl")
 
 export heatconduction!, enthalpy, totalwater, liquidwater, freezethaw!, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!

src/Physics/HeatConduction/heat.jl

@@ -37,7 +37,7 @@ Get heat capacities for generic `SubSurface` layer.
 """
 @inline heatcapacities(::SubSurface, heat::Heat) = (heat.prop.cw, heat.prop.ci, heat.prop.ca)
 """
-    volumetricfractions(::SubSurface, heat::Heat)
+    volumetricfractions(::SubSurface, heat::Heat, state, i)
 
 Get constituent volumetric fractions for generic `SubSurface` layer. Default implementation assumes
 the only constitutents are liquid water, ice, and air: `(θl,θi,θa)`.
@@ -270,18 +270,15 @@ function interact!(sub1::SubSurface, ::Heat, sub2::SubSurface, ::Heat, s1, s2)
     Δk₁ = thickness(sub1, s1)
     Δk₂ = thickness(sub2, s2)
     # thermal conductivity between cells
-    @inbounds let k₁ = s1.kc[end],
-        k₂ = s2.kc[1],
-        Δ₁ = Δk₁[end],
-        Δ₂ = Δk₂[1];
-        k = harmonicmean(k₁, k₂, Δ₁, Δ₂);
-        s1.k[end] = s2.k[1] = k
-    end
+    k = s1.k[end] = s2.k[1] =
+        @inbounds let k₁ = s1.kc[end],
+            k₂ = s2.kc[1],
+            Δ₁ = Δk₁[end],
+            Δ₂ = Δk₂[1];
+            harmonicmean(k₁, k₂, Δ₁, Δ₂)
+        end
     # calculate heat flux between cells
-    Qᵢ = @inbounds let k₁ = s1.k[end],
-        k₂ = s2.k[1],
-        k = k₁ = k₂,
-        z₁ = midpoints(sub1, s1),
+    Qᵢ = @inbounds let z₁ = midpoints(sub1, s1),
         z₂ = midpoints(sub2, s2),
         δ = z₂[1] - z₁[end];
         k*(s2.T[1] - s1.T[end]) / δ

src/Physics/Physics.jl

@@ -25,6 +25,7 @@ include("coupled.jl")
 include("Boundaries/Boundaries.jl")
 include("HeatConduction/HeatConduction.jl")
 include("WaterBalance/WaterBalance.jl")
+include("Snow/Snow.jl")
 include("Soils/Soils.jl")
 include("SEB/SEB.jl")
 include("Sources/Sources.jl")
@@ -32,6 +33,7 @@ include("Sources/Sources.jl")
 @reexport using .Boundaries
 @reexport using .HeatConduction
 @reexport using .WaterBalance
+@reexport using .Snow
 @reexport using .Soils
 @reexport using .SEB
 @reexport using .Sources

src/Physics/Snow/Snow.jl

+module Snow
+
+using CryoGrid: Top, SubSurfaceProcess, BoundaryProcess, BoundaryStyle, Dirichlet
+using CryoGrid.InputOutput: Forcing
+using CryoGrid.Physics.HeatConduction
+using CryoGrid.Numerics
+using CryoGrid.Utils
+
+import CryoGrid
+import CryoGrid.Physics.HeatConduction
+
+using Unitful
+
+export Snowpack, SnowThermalProperties
+
+"""
+Base type for snow paramterization schemes.
+"""
+abstract type SnowParameterization <: CryoGrid.Parameterization end
+
+"""
+    Bulk{Tdsn} <: SnowParameterization
+
+Simple, bulk snow scheme where snowpack is represented as a single grid cell with homogenous state.
+"""
+Base.@kwdef struct Bulk{Tdsn,Tthresh} <: SnowParameterization
+    dsn::Tdsn = 0.0u"m" # snow depth
+    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()
+
+"""
+    Snowpack{TPara<:SnowParameterization,TProp,TSp} <: CryoGrid.SubSurface
+
+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
+
+include("snow_bulk.jl")
+
+end
\ No newline at end of file

src/Physics/Snow/snow_bulk.jl

+"""
+    BulkSnowpack{Tdsn} = Snowpack{<:Bulk{Tdsn}} where {Tdsn}
+
+Type alias for Snowpack with `Bulk` parameterization.
+"""
+const BulkSnowpack{Tdsn} = Snowpack{<:Bulk{Tdsn}} where {Tdsn}
+# Local alias for HeatConduction Enthalpy type
+const Enthalpy = HeatConduction.Enthalpy
+
+CryoGrid.variables(snow::BulkSnowpack) = (
+    Diagnostic(:θw, OnGrid(Cells)),
+    Diagnostic(:T_ub, Scalar, u"°C"),
+    Diagnostic(:dsn, Scalar, u"m"),
+)
+
+CryoGrid.callbacks(snow::BulkSnowpack, heat::Heat) = (
+    SnowFree(),
+)
+function CryoGrid.criterion(::SnowFree, snow::BulkSnowpack, heat::Heat, state)
+    dsn = snowdepth(snow, state)
+    # snow depth shifted down by 1 cm (TODO: should be configurable?);
+    # callback will fire when snow depth drops below or rises above this threshold.
+    return dsn - one(dsn)*snow.para.thresh
+end
+function CryoGrid.affect!(::SnowFree, snow::BulkSnowpack, heat::Heat, state)
+    @. state.θw = 0.0
+    @. state.H = 0.0
+end
+
+snowdepth(snow::BulkSnowpack{<:Number}, state) = snow.para.dsn
+snowdepth(snow::BulkSnowpack{<:Forcing}, state) = snow.para.dsn(state.t)
+
+# HeatConduction methods
+function HeatConduction.freezethaw!(snow::Snowpack, heat::Heat{FreeWater,Enthalpy}, state)
+    @inbounds for i in 1:length(state.H)
+        # liquid water content = (total water content) * (liquid fraction)
+        state.θl[i] = HeatConduction.liquidwater(snow, heat, state, i)
+        # update heat capacity
+        state.C[i] = C = snow.prop.csn
+        # enthalpy inverse function
+        state.T[i] = HeatConduction.enthalpyinv(snow, heat, state, i)
+        # set dHdT (a.k.a dHdT)
+        state.dHdT[i] = state.T[i] ≈ 0.0 ? 1e8 : 1/C
+    end
+    return nothing
+end
+@inline HeatConduction.thermalconductivities(snow::Snowpack, ::Heat) = (snow.prop.ksn,)
+@inline HeatConduction.heatcapacities(snow::Snowpack, ::Heat) = (snow.prop.csn,)
+
+# CryoGrid methods
+CryoGrid.volumetricfractions(::BulkSnowpack, ::Heat, state, i) = (1.0,)
+CryoGrid.thickness(::BulkSnowpack, state) = getscalar(state.dsn)
+CryoGrid.midpoints(::BulkSnowpack, state) = -getscalar(state.dsn) / 2
+function CryoGrid.diagnosticstep!(snow::BulkSnowpack, state)
+    @setscalar state.dsn = snowdepth(snow, state)
+end
+function CryoGrid.diagnosticstep!(snow::BulkSnowpack, heat::Heat{FreeWater,Enthalpy}, state)
+    dsn = getscalar(state.dsn)
+    @. state.θw = snow.prop.ρsn / heat.prop.ρw
+    if dsn > snow.para.thresh
+        @. state.θw = snow.prop.ρsn / heat.prop.ρw
+        HeatConduction.freezethaw!(snow, heat, state)
+    else
+        @. state.θw = 0.0
+        @. state.T = state.T_ub
+        @. state.C = heat.prop.ca
+    end
+    HeatConduction.thermalconductivity!(snow, heat, state)
+    @. state.k = state.kc
+    return nothing
+end
+function CryoGrid.prognosticstep!(snow::BulkSnowpack, heat::Heat{FreeWater,Enthalpy}, state)
+    dsn = getscalar(state.dsn)
+    if dsn < snow.para.thresh
+        # set energy flux to zero if there is no snow
+        @. state.dH = zero(eltype(state.H))
+    end
+end
+# Special implementation of interact! only for Dirichlet boundary temperature;
+# We do this only so that T_ub can be set to the value of the boundary condition
+CryoGrid.interact!(top::Top, bc::HeatBC, snow::BulkSnowpack, heat::Heat, stop, ssnow) = CryoGrid.interact!(BoundaryStyle(bc), top, bc, snow, heat, stop, ssnow)
+function CryoGrid.interact!(::Dirichlet, top::Top, bc::HeatBC, sub::BulkSnowpack, heat::Heat, stop, ssnow)
+    Δk = CryoGrid.thickness(sub, ssnow) # using `thickness` allows for generic layer implementations
+    @setscalar ssnow.T_ub = CryoGrid.boundaryvalue(bc, top, heat, sub, stop, ssnow)
+    # boundary flux
+    @setscalar ssnow.dH_upper = CryoGrid.boundaryflux(bc, top, heat, sub, stop, ssnow)
+    @inbounds ssnow.dH[1] += getscalar(ssnow.dH_upper) / Δk[1]
+    return nothing # ensure no allocation
+end

src/Physics/Soils/Soils.jl

 module Soils
 
 import CryoGrid: SubSurface, Parameterization
-import CryoGrid: initialcondition!, variables
+import CryoGrid: initialcondition!, variables, volumetricfractions
 import CryoGrid.Physics: totalwater
-import CryoGrid.Physics.HeatConduction: Enthalpy, Temperature, liquidwater, thermalconductivity, heatcapacity, freezethaw!, enthalpyinv
+import CryoGrid.Physics.HeatConduction: Enthalpy, Temperature, liquidwater, thermalconductivity, heatcapacity, freezethaw!, enthalpyinv, heatcapacities, thermalconductivities
 
 using CryoGrid.Numerics
 using CryoGrid.Numerics: heaviside

src/methods.jl

@@ -142,3 +142,12 @@ midpoints(::Layer, state) = cells(state.grid)
 Get thickness (m) of layer or all grid cells within the layer.
 """
 thickness(::Layer, state) = Δ(state.grid)
+# Composition
+"""
+    volumetricfractions(::Layer, ::Process, state)
+    volumetricfractions(::Layer, ::Process, state, i)
+
+Get the volumetric fractions of each constituent in the volume (at grid cell `i`, if specificed).
+"""
+volumetricfractions(::Layer, ::Process, state) = ()
+volumetricfractions(::Layer, ::Process, state, i) = ()