figure changes and freezing point depression (#963)

Author: kmdeck

docs/tutorials/standalone/Soil/freezing_front.jl

@@ -104,8 +104,8 @@ hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n);
 #c = FT(0.43)
 #hcm = BrooksCorey(;ψb = ψb, c = c);
 θ_r = FT(0.05)
-ν_ss_om = FT(0.3)
-ν_ss_quartz = FT(0.7)
+ν_ss_om = FT(0.4)
+ν_ss_quartz = FT(0.6)
 ν_ss_gravel = FT(0.0)
 params = Soil.EnergyHydrologyParameters(
     FT;

docs/tutorials/standalone/Soil/phase_change_analytic.jl

@@ -180,7 +180,7 @@ ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
         max_iters = 3,
-        update_j = CTS.UpdateEvery(CTS.NewTimeStep),
+        update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );
 
@@ -212,17 +212,12 @@ sol = SciMLBase.solve(
 sol_T = parent(sv.saveval[end].soil.T)[:]
 
 fig = Figure(size = (300, 300))
-ax1 = Axis(
-    fig[1, 1],
-    title = "Neumann Solution",
-    xlabel = L"T (K)",
-    ylabel = "Soil depth (m)",
-)
+ax1 = Axis(fig[1, 1], title = "", xlabel = L"T (K)", ylabel = "Soil depth (m)")
 limits!(ax1, 262, 276, -3, 0.0)
 
 z = parent(coords.subsurface.z)[:];
 
-lines!(ax1, sol_T, z, label = "Model", color = :green)
+lines!(ax1, sol_T, z, label = "Simulation", color = :green)
 
 # # Analytic Solution of Neumann
 # All details here are taken from [DallAmico2011](@cite) (see also [CarslawJaeger](@cite)), and the reader is referred to that
@@ -309,7 +304,7 @@ function implicit(ζ)
     return (term1 + term2 + term3)
 end
 ζ = find_zero(implicit, (0.25, 0.27), Bisection())
-depth = abs.(reverse(z))
+depth = Array(0:0.01:3)
 t = sol.t[end]
 zf = 2.0 * ζ * sqrt(d1 * t)
 analytic_unfrozen_profile(depth, zf) =
@@ -319,14 +314,14 @@ mask_unfrozen = depth .>= zf
 mask_frozen = depth .<= zf
 T_frozen = Ts .+ (0.0 - Ts) .* analytic_frozen_profile.(depth, zf) .+ 273.15
 T_unfrozen = Ti .- (Ti - 0.0) .* analytic_unfrozen_profile.(depth, zf) .+ 273.15
-scatter!(
+lines!(
     ax1,
     T_frozen[mask_frozen],
     -1 .* depth[mask_frozen],
-    label = "Analytic",
+    label = "Analytic Solution",
     color = "purple",
 )
-scatter!(
+lines!(
     ax1,
     T_unfrozen[mask_unfrozen],
     -1 .* depth[mask_unfrozen],

lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl

@@ -35,7 +35,7 @@ function ozark_default_args(;
     h_stem = FT(9),
     h_leaf = FT(9.5),
     t0 = Float64(120 * 3600 * 24),
-    dt = Float64(900),
+    dt = Float64(450),
     n = 15,
 )
     return (

src/standalone/Soil/soil_heat_parameterizations.jl

@@ -65,8 +65,9 @@ function phase_change_source(
     θtot = min(_ρ_i / _ρ_l * θ_i + θ_l, ν)
     # This is consistent with Equation (22) of Dall'Amico
     ψw0 = matric_potential(hydrology_cm, effective_saturation(ν, θtot, θ_r))
+    Tf_depressed = _T_freeze * exp(_grav * ψw0 / _LH_f0)
 
-    ψT = _LH_f0 / _T_freeze / _grav * (T - _T_freeze) * heaviside(_T_freeze - T)
+    ψT = _LH_f0 / _grav * log(T / Tf_depressed) * heaviside(Tf_depressed - T)
     # Equation (23) of Dall'Amico
     θstar = inverse_matric_potential(hydrology_cm, ψw0 + ψT) * (ν - θ_r) + θ_r
 

test/standalone/Soil/soil_parameterizations.jl

@@ -332,10 +332,13 @@ for FT in (Float32, Float64)
         @test τ == parameters.ρc_ds * Δz^2 / parameters.κ_dry
         θ_l = FT.([0.11, 0.15, ν])
         θ_i = FT(0.0)
-        T = FT(273)
+        T = FT(270)
         θtot = @.(_ρ_i / _ρ_l * θ_i + θ_l)
         ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r))
-        ψT = @.(_LH_f0 / _T_freeze / _grav * (T - _T_freeze))
+        Tf_depressed = _T_freeze * exp.(_grav * ψ0 / _LH_f0)
+        ψT = @. _LH_f0 / _grav *
+           log(T / Tf_depressed) *
+           ClimaLand.heaviside(Tf_depressed - T)
         θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r
         ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set)
         τ = thermal_time.(ρc_s, Δz, parameters.κ_dry)
@@ -361,14 +364,16 @@ for FT in (Float32, Float64)
                 hydrology_earth_params,
             ) .> 0.0,
         ) == 3
-        # try θ_l = 0.1
 
         θ_l = FT.([0.11, 0.15, ν])
         θ_i = FT(0.0)
         T = FT(274)
         θtot = @.(_ρ_i / _ρ_l * θ_i + θ_l)
         ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r))
-        ψT = FT(0.0)
+        Tf_depressed = _T_freeze * exp.(_grav * ψ0 / _LH_f0)
+        ψT = @. _LH_f0 / _grav *
+           log(T / Tf_depressed) *
+           ClimaLand.heaviside(Tf_depressed - T)
         θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r
         ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set)
         τ = thermal_time.(ρc_s, Δz, parameters.κ_dry)
@@ -386,12 +391,15 @@ for FT in (Float32, Float64)
         @test (θ_star ≈ θ_l)
 
 
-        θ_l = FT(0.01)
+        θ_l = FT(0.11)
         θ_i = FT.([0.05, 0.08])
         T = FT(274)
         θtot = @.(_ρ_i / _ρ_l * θ_i + θ_l)
         ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r))
-        ψT = FT(0.0)
+        Tf_depressed = _T_freeze * exp.(_grav * ψ0 / _LH_f0)
+        ψT = @. _LH_f0 / _grav *
+           log(T / Tf_depressed) *
+           ClimaLand.heaviside(Tf_depressed - T)
         θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r
         ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set)
         τ = thermal_time.(ρc_s, Δz, parameters.κ_dry)
@@ -424,7 +432,10 @@ for FT in (Float32, Float64)
         T = FT(260)
         θtot = @.(_ρ_i / _ρ_l * θ_i + θ_l)
         ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r))
-        ψT = @.(_LH_f0 / _T_freeze / _grav * (T - _T_freeze))
+        Tf_depressed = _T_freeze * exp.(_grav * ψ0 / _LH_f0)
+        ψT = @. _LH_f0 / _grav *
+           log(T / Tf_depressed) *
+           ClimaLand.heaviside(Tf_depressed - T)
         θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r
         ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set)
         τ = thermal_time.(ρc_s, Δz, parameters.κ_dry)