Use SoilProfile for heterogeneous soil layer types

examples/heat_freeW_seb_snow_bucketW_samoylov.jl

@@ -31,11 +31,7 @@ water = WaterBalance(BucketScheme(), DampedET())
 strat = @Stratigraphy(
     -z => Top(upperbc), 
     -z => Snowpack(heat=HeatBalance(), water=water),
-    soilprofile[1].depth => Ground(soilprofile[1].value; heat, water),
-    soilprofile[2].depth => Ground(soilprofile[2].value; heat, water),
-    soilprofile[3].depth => Ground(soilprofile[3].value; heat, water),
-    soilprofile[4].depth => Ground(soilprofile[4].value; heat, water),
-    soilprofile[5].depth => Ground(soilprofile[5].value; heat, water),
+    0.0u"m" => Ground(soilprofile; heat, water),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")),
 );
 ## create Tile
@@ -44,7 +40,7 @@ tile = Tile(strat, modelgrid, initT);
 # Set up the problem and solve it!
 tspan = (DateTime(2010,10,30), DateTime(2011,10,30))
 ## generate initial condition and set up CryoGridProblem
-@time u0, du0 = initialcondition!(tile, tspan)
+u0, du0 = @time initialcondition!(tile, tspan)
 
 prob = CryoGridProblem(
     tile,
@@ -55,7 +51,7 @@ prob = CryoGridProblem(
 )
 integrator = init(prob, Euler(), dt=60.0)
 ## step forwards 24 hours and check for NaN/Inf values
-@time step!(integrator, 24*3600)
+@time step!(integrator); integrator.dt
 @assert all(isfinite.(integrator.u))
 ## iterate over remaining timespan at fixed points using `TimeChoiceIterator`
 @time for (u,t) in TimeChoiceIterator(integrator, convert_t.(tspan[1]:Day(1):tspan[end]))

src/Numerics/profile.jl

@@ -1,14 +1,14 @@
-struct ProfileKnot{D,T}
+struct ProfileKnot{D,V}
     depth::D
-    value::T
+    value::V
 end
 Base.iterate(knot::ProfileKnot) = iterate((knot.depth, knot.value))
 Base.iterate(knot::ProfileKnot, state) = iterate((knot.depth, knot.value), state)
 Base.show(io::IO, knot::ProfileKnot) = print(io, "$(knot.depth): $(knot.value)")
-struct Profile{N,TKnots}
-    knots::TKnots
-    Profile(::Tuple{}) = new{0,Tuple{}}(())
-    Profile(knots::Tuple{Vararg{ProfileKnot,N}}) where N = new{N,typeof(knots)}(knots)
+struct Profile{N,V,D}
+    knots::NTuple{N,ProfileKnot{D,V}}
+    Profile(::Tuple{}) = new{0,Nothing,Nothing}(())
+    Profile(knots::NTuple{N,ProfileKnot{D,V}}) where {N,D,V} = new{N,V,D}(knots)
     Profile(pairs::Tuple{Vararg{Pair}}) = Profile(map(Base.splat(ProfileKnot), pairs))
     Profile(pairs::Pair...) = Profile(pairs)
 end
@@ -21,6 +21,7 @@ end
 Base.length(::Profile{N}) where N = N
 Base.iterate(profile::Profile) = iterate(profile.knots)
 Base.iterate(profile::Profile, state) = iterate(profile.knots, state)
+Base.map(f, profile::Profile) = Profile(map(knot -> ProfileKnot(knot.depth, f(knot.value)), profile.knots))
 Base.getindex(profile::Profile, itrv::Interval) = Profile(Tuple(knot for knot in profile.knots if knot.depth ∈ itrv))
 Base.getindex(profile::Profile, i::Int) = profile.knots[i]
 Base.getindex(profile::Profile, i) = Profile(profile.knots[i])

src/Physics/Heat/heat_methods.jl

@@ -6,8 +6,20 @@ Returns the thermal properties for the given subsurface layer.
 """
 thermalproperties(::SubSurface) = error("not implemented")
 thermalproperties(sub::SubSurface, state) = thermalproperties(sub)
-thermalproperties(sub::SubSurface, state, i) = thermalproperties(sub, state, 1)
+thermalproperties(sub::SubSurface, state, i) = thermalproperties(sub, state)
 
+"""
+    thermalconductivity(::SubSurface, state, i)
+
+Computes the thermal conductivity for the given `SubSurface` layer at grid cell `i`.
+"""
+thermalconductivity(::SubSurface, state, i) = error("not implemented")
+
+"""
+    heatcapacity(::SubSurface, state, i)
+
+Computes the heat capacity for the given `SubSurface` layer at grid cell `i`.
+"""
 heatcapacity(::SubSurface, state, i) = error("not implemented")
 
 """

src/Physics/Heat/heat_water.jl

@@ -19,7 +19,7 @@ Adds advective energy fluxes for all internal grid cell faces.
 """
 function water_energy_advection!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
     water, heat = ps
-    cw = heatcapacity_water(sub)
+    cw = heatcapacitywater(sub, state)
     @inbounds for i in 2:length(state.jw)-1
         let jw = state.jw[i],
             T₁ = state.T[i-1],
@@ -42,7 +42,7 @@ function CryoGrid.interact!(top::Top, bc::WaterBC, sub::SubSurface, heat::HeatBa
     if heat.advection
         T_ub = getscalar(stop.T_ub)
         Ts = ssub.T[1]
-        cw = heatcapacity_water(sub)
+        cw = heatcapacitywater(sub, ssub)
         L = heat.prop.L
         jw = ssub.jw[1]
         jH_w = advectiveflux(jw, T_ub, Ts, cw, L)
@@ -54,7 +54,7 @@ function CryoGrid.interact!(sub::SubSurface, heat::HeatBalance, bot::Bottom, bc:
     if heat.advection
         Ts = ssub.T[end]
         T_lb = getscalar(sbot.T_ub)
-        cw = heatcapacity_water(sub)
+        cw = heatcapacitywater(sub, ssub)
         L = heat.prop.L
         jw = ssub.jw[end]
         jH_w = advectiveflux(jw, Ts, T_lb, cw, L)
@@ -80,8 +80,8 @@ function CryoGrid.interact!(
         T₂ = state2.T[1]
         jw = state1.jw[end]
         L = p1[2].prop.L
-        cw1 = heatcapacity_water(sub1)
-        cw2 = heatcapacity_water(sub2)
+        cw1 = heatcapacitywater(sub1, state1)
+        cw2 = heatcapacitywater(sub2, state2)
         cw = (jw >= zero(jw))*cw1 + (jw < zero(jw))*cw2
         jH_w = advectiveflux(jw, T₁, T₂, cw, L)
         state1.jH[end] = state2.jH[1] += jH_w
@@ -98,7 +98,15 @@ function CryoGrid.computefluxes!(sub::SubSurface, ps::Coupled(WaterBalance, Heat
     CryoGrid.computefluxes!(sub, heat, state)
 end
 
-function heatcapacity_water(sub::SubSurface)
-    @unpack ch_w = thermalproperties(sub)
-    return ch_w
+function heatcapacitywater(sub::SubSurface, state)
+    if hasproperty(state, :ch_w)
+        # TODO: this is a temporary solution that works so long as the
+        # heat capacity of water is constant throughout the layer.
+        # This should generally be the case, but it's not a very elegant
+        # solution and it would be better to treat such constants more properly.
+        return state.ch_w[1]
+    else
+        @unpack ch_w = thermalproperties(sub)
+        return ch_w
+    end
 end

src/Physics/Hydrology/water_balance.jl

@@ -1,3 +1,40 @@
+"""
+    watercontent!(::SubSurface, ::WaterBalance, state)
+
+Computes the volumetric water content from current saturation or pressure state.
+"""
+function watercontent!(sub::SubSurface, water::WaterBalance, state)
+    @inbounds for i in eachindex(state.sat)
+        state.θsat[i] = maxwater(sub, water, state, i)
+        state.θwi[i] = state.sat[i]*state.θsat[i]
+        # initially set unfrozen water = water+ice;
+        # when heat conduction is included, it should update this based on temperature.
+        state.θw[i] = state.θwi[i]
+    end
+end
+
+"""
+    hydraulicconductivity!(sub::SubSurface, water::WaterBalance, state)
+
+Computes hydraulic conductivities for the given subsurface layer and water balance scheme.
+"""
+function hydraulicconductivity!(sub::SubSurface, water::WaterBalance, state)
+    @inbounds for i in eachindex(state.kwc)
+        let θsat = Hydrology.maxwater(sub, water, state, i),
+            θw = state.θw[i],
+            θwi = state.θwi[i],
+            kw_sat = kwsat(sub, state, i);
+            state.kwc[i] = hydraulicconductivity(water, kw_sat, θw, θwi, θsat)
+            if i > 1
+                state.kw[i] = min(state.kwc[i-1], state.kwc[i])
+            end
+        end
+    end
+    # set hydraulic conductivity at boundaries equal to cell conductivities
+    state.kw[1] = state.kwc[1]
+    state.kw[end] = state.kwc[end]
+end
+
 """
     limit_upper_flux(water::WaterBalance, jw, θw, θwi, θsat, sat, Δz)
 

src/Physics/Hydrology/water_methods.jl

@@ -1,10 +1,3 @@
-"""
-    hydraulicproperties(::SubSurface)
-
-Retrieves the hydraulic properties from the given subsurface layer.
-"""
-hydraulicproperties(::SubSurface) = error("not implemented")
-
 """
     waterdensity(sub::SubSurface)
 
@@ -14,11 +7,22 @@ Retrieves the density of water `ρw` from the given `SubSurface` layer. Default
 waterdensity(sub::SubSurface) = sub.water.prop.ρw
 
 """
-    kwsat(::SubSurface, ::WaterBalance)
+    hydraulicproperties(::SubSurface)
+
+Retrieves the hydraulic properties from the given subsurface layer.
+"""
+hydraulicproperties(::SubSurface) = error("not implemented")
+hydraulicproperties(sub::SubSurface, state) = hydraulicproperties(sub)
+hydraulicproperties(sub::SubSurface, state, i) = hydraulicproperties(sub, state)
+
+"""
+    kwsat(::SubSurface, state, i)
 
 Hydraulic conductivity at saturation.
 """
-kwsat(sub::SubSurface, ::WaterBalance) = hydraulicproperties(sub).kw_sat
+kwsat(sub::SubSurface) = hydraulicproperties(sub).kw_sat
+kwsat(sub::SubSurface, state) = hydraulicproperties(sub, state).kw_sat
+kwsat(sub::SubSurface, state, i) = hydraulicproperties(sub, state, i).kw_sat
 
 """
     interact_ET!(::Top, ::WaterBC, ::SubSurface, ::WaterBalance, state1, state2)
@@ -49,7 +53,7 @@ maxwater(sub::SubSurface, water::WaterBalance, state, i) = Utils.getscalar(maxwa
 
 Returns the minimum volumetric water content (typically field capacity for simplified schemes) for grid cell `i`. Defaults to zero.
 """
-minwater(sub::SubSurface, ::WaterBalance) = hydraulicproperties(sub).fieldcapacity
+minwater(sub::SubSurface, ::WaterBalance) = sqrt(eps())
 minwater(sub::SubSurface, water::WaterBalance, state) = minwater(sub, water)
 minwater(sub::SubSurface, water::WaterBalance, state, i) = Utils.getscalar(minwater(sub, water, state), i)
 
@@ -63,45 +67,9 @@ watercontent(sub::SubSurface, state) = state.θwi
 watercontent(sub::SubSurface, state, i) = Utils.getscalar(watercontent(sub, state), i)
 
 """
-    watercontent!(::SubSurface, ::WaterBalance, state)
-
-Computes the volumetric water content from current saturation or pressure state.
-"""
-function watercontent!(sub::SubSurface, water::WaterBalance, state)
-    @inbounds for i in eachindex(state.sat)
-        state.θsat[i] = maxwater(sub, water, state, i)
-        state.θwi[i] = state.sat[i]*state.θsat[i]
-        # initially set unfrozen water = water+ice;
-        # when heat conduction is included, it should update this based on temperature.
-        state.θw[i] = state.θwi[i]
-    end
-end
-
-"""
-    hydraulicconductivity(sub::SubSurface, water::WaterBalance, θw, θwi, θsat)
+    hydraulicconductivity(water::WaterBalance, kw_sat, θw, θwi, θsat)
 
-Computes the hydraulic conductivity for the given layer and water balance configuration, current unfrozen
+Computes the hydraulic conductivity for the given water balance configuration, current unfrozen
 water content `θw`, total water/ice content `θwi`, and saturated (maximum) water content `θsat`.
 """
-hydraulicconductivity(sub::SubSurface, water::WaterBalance, θw, θwi, θsat) = kwsat(sub, water)*θw / θsat
-
-"""
-    hydraulicconductivity!(sub::SubSurface, water::WaterBalance, state)
-
-Computes hydraulic conductivities for the given subsurface layer and water balance scheme.
-"""
-function hydraulicconductivity!(sub::SubSurface, water::WaterBalance, state)
-    @inbounds for i in eachindex(state.kwc)
-        let θsat = Hydrology.maxwater(sub, water, state, i),
-            θw = state.θw[i],
-            θwi = state.θwi[i];
-            state.kwc[i] = hydraulicconductivity(sub, water, θw, θwi, θsat)
-            if i > 1
-                state.kw[i] = min(state.kwc[i-1], state.kwc[i])
-            end
-        end
-    end
-    # set hydraulic conductivity at boundaries equal to cell conductivities
-    state.kw[1] = state.kwc[1]
-    state.kw[end] = state.kwc[end]
-end
+hydraulicconductivity(::WaterBalance, kw_sat, θw, θwi, θsat) = kw_sat*θw / θsat

src/Physics/Hydrology/water_types.jl

@@ -7,8 +7,6 @@ Base.@kwdef struct WaterBalanceProperties{Tρw,TLsg,Trfs}
     ρw::Tρw = CryoGrid.Constants.ρw
     Lsg::TLsg = CryoGrid.Constants.Lsg
     rf_smoothness::Trfs = 0.3
-    # r_β::Trb = 1e3 # reduction factor scale parameter
-    # r_c::Trc = 0.96325 # reduction factor shift parameter
 end
 
 # do not parameterize water balance constants
@@ -21,7 +19,6 @@ Default material hydraulic properties.
 """
 Utils.@properties HydraulicProperties(
     kw_sat = 1e-5u"m/s",
-    fieldcapacity = 0.05,
 )
 
 function CryoGrid.parameterize(prop::HydraulicProperties)

src/Physics/Salt/salt_diffusion.jl

@@ -7,10 +7,9 @@ function Heat.freezethaw!(
 ) where {THeat<:HeatBalance{<:DallAmicoSalt,<:TemperatureBased}}
     salt, heat = ps
     sfcc = heat.freezecurve
-    thermalprops = Heat.thermalproperties(soil)
-    @unpack ch_w, ch_i = thermalprops
     let L = heat.prop.L;
         @inbounds @fastmath for i in eachindex(state.T)
+            @unpack ch_w, ch_i = Heat.thermalproperties(soil, state, i)
             T = state.T[i]
             c = state.c[i]
             θsat = Soils.porosity(soil, state, i)

src/Physics/Snow/snow_bulk.jl

@@ -21,26 +21,6 @@ const EnthalpyBased = Heat.EnthalpyBased
 
 threshold(snow::BulkSnowpack) = snow.para.thresh
 
-function snowdensity!(
-    snow::BulkSnowpack{<:ConstantDensity},
-    mass::DynamicSnowMassBalance,
-    state
-)
-    ρsn = snow.para.density.ρsn
-    ρw = waterdensity(snow)
-    state.ρsn .= ρsn
-    state.por .= 1 - ρsn / ρw
-    return nothing
-end
-
-function snowdepth!(
-    ::BulkSnowpack,
-    ::DynamicSnowMassBalance,
-    state
-)
-    @setscalar state.dsn = getscalar(state.Δz)
-end
-
 # implement ablation! for DegreeDayMelt
 function ablation!(
     ::Top,
@@ -148,7 +128,7 @@ function CryoGrid.trigger!(
     # using current upper boundary temperature.
     heatcapacity!(snow, snow.heat, state)
     state.T .= state.T_ub
-    state.H .= state.T.*C
+    state.H .= state.T.*state.C
     state.por .= 1 - getscalar(state.ρsn) / waterdensity(snow)
     state.sat .= zero(eltype(state.sat))
     return nothing

src/Physics/Snow/snow_density.jl

@@ -4,3 +4,15 @@ abstract type SnowDensityScheme end
 Base.@kwdef struct ConstantDensity{Tρsn} <: SnowDensityScheme
     ρsn::Tρsn = 250.0u"kg/m^3" # constant snow density
 end
+
+function snowdensity!(
+    snow::Snowpack{<:ConstantDensity},
+    mass::DynamicSnowMassBalance,
+    state
+)
+    ρsn = snow.para.density.ρsn
+    ρw = waterdensity(snow)
+    state.ρsn .= ρsn
+    state.por .= 1 - ρsn / ρw
+    return nothing
+end

src/Physics/Snow/snow_methods.jl

@@ -70,6 +70,14 @@ snowvariables(::Snowpack) = (
     Diagnostic(:T_ub, Scalar, u"°C"),
 )
 
+function snowdepth!(
+    ::Snowpack,
+    ::DynamicSnowMassBalance,
+    state
+)
+    @setscalar state.dsn = getscalar(state.Δz)
+end
+
 ### Default implmentations of CryoGrid methods for Snowpack ###
 
 # implement CryoGrid.Volume for prescribed vs. dynamic snow mass balance

src/Physics/Soils/Soils.jl

@@ -38,11 +38,11 @@ include("soil_methods.jl")
 export SoilTexture
 include("soil_texture.jl")
 
-export MineralOrganic, soilcomponent
-include("para/mineral_organic.jl")
+export Heterogeneous, MineralOrganic, soilcomponent
+include("soil_para.jl")
 
 export SURFEX
-include("para/surfex.jl")
+include("soil_surfex.jl")
 
 export RichardsEq
 include("soil_water.jl")