add snow surface temp parameterization

Project.toml

@@ -16,6 +16,7 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f"

experiments/standalone/Snow/snowmip_simulation.jl

@@ -107,6 +107,7 @@ sol = SciMLBase.solve(
 # Plotting
 q_l = [parent(sv.saveval[k].snow.q_l)[1] for k in 1:length(sol.t)];
 T = [parent(sv.saveval[k].snow.T)[1] for k in 1:length(sol.t)];
+T_sfc = [parent(sv.saveval[k].snow.T_sfc)[1] for k in 1:length(sol.t)];
 evaporation = [
     parent(sv.saveval[k].snow.turbulent_fluxes.vapor_flux)[1] for
     k in 1:length(sol.t)
@@ -126,7 +127,6 @@ t = sol.t;
 start_day = 1
 days = start_day .+ floor.(t ./ 3600 ./ 24)
 doys = days .% 365 # doesn't account for leap year
-
 obs_swes = Vector{Union{Float64, Missing}}(missing, length(doys))
 obs_swes[snow_data_avail] .= mass[snow_data_avail] ./ 1000
 
@@ -137,13 +137,30 @@ obs_df = DataFrame(
     doy = doys,
     model_swe = S,
     obs_swe = obs_swes,
-    model_tsnow = T,
+    model_tsnow = T_sfc,
     obs_tsnow = obs_tsnows,
 )
 function missingmean(x)
     return mean(skipmissing(x))
 end
 
+function rmse(model_vals, obs_vals)
+    return sqrt(mean((model_vals .- obs_vals) .^ 2))
+end
+
+filtered_df = obs_df[snow_data_avail .== 1, :]
+println(
+    SITE_NAME,
+    " SWE RMSE: ",
+    rmse(filtered_df[!, :model_swe], filtered_df[!, :obs_swe]),
+)
+snowtemp_obs = filtered_df[.!ismissing.(filtered_df[!, :obs_tsnow]), :]
+println(
+    SITE_NAME,
+    " Tsfc RMSE: ",
+    rmse(snowtemp_obs[!, :model_tsnow], snowtemp_obs[!, :obs_tsnow]),
+)
+
 mean_obs_df = combine(
     groupby(obs_df, :doy),
     [:model_swe, :obs_swe, :model_tsnow, :obs_tsnow] .=> missingmean,
@@ -165,6 +182,7 @@ CairoMakie.scatter!(
     label = "Data",
     color = :red,
 )
+
 CairoMakie.axislegend(ax1, position = :rt)
 
 ax2 = CairoMakie.Axis(
@@ -198,12 +216,12 @@ ax1 = CairoMakie.Axis(fig[1, 1], ylabel = "Temperature")
 
 xlims!(ax1, 0, ndays)
 CairoMakie.hidexdecorations!(ax1, ticks = false)
-CairoMakie.lines!(ax1, daily, T, label = "Model")
+CairoMakie.lines!(ax1, daily, T_sfc, label = "Snow Surface, Model")
 CairoMakie.scatter!(
     ax1,
     seconds[snow_data_avail] ./ 24 ./ 3600,
     T_snow .+ 273.15,
-    label = "Snow T",
+    label = "Snow Bulk, Model",
     color = :red,
 )
 CairoMakie.scatter!(
@@ -227,6 +245,7 @@ ax3 = CairoMakie.Axis(
     ylabel = "Average Snow Temperature",
     xlabel = "Day of Year",
 )
+xlims!(ax3, 245, 281)
 CairoMakie.lines!(
     ax3,
     mean_obs_df.doy,
@@ -274,6 +293,7 @@ CairoMakie.lines!(
     label = "Runoff/Melt",
     color = :blue,
 )
+
 CairoMakie.axislegend(ax1, position = :lb)
 
 CairoMakie.save(joinpath(savedir, "water_fluxes_$(SITE_NAME).png"), fig)
@@ -336,3 +356,28 @@ CairoMakie.lines!(
 )
 
 CairoMakie.save(joinpath(savedir, "conservation_$(SITE_NAME).png"), fig)
+
+#=
+plot4 = plot(
+    xlabel = "Time (days)",
+    ylabel = "Snow Temperature",
+    xlim = (420, 434),
+    ylim = (260, 280),
+)
+plot!(
+    plot4,
+    t ./ 3600 ./ 24,
+    T_sfc,
+    legend = :bottomright,
+    label = "Model",
+    ylabel = "Snow Temperature",
+)
+scatter!(
+    plot4,
+    seconds[snow_data_avail] ./ 3600 ./ 24,
+    T_snow .+ 273.15,
+    label = "Data",
+)
+savefig(joinpath(savedir, "tsfc_1week_$(SITE_NAME).png"))
+>>>>>>> 7d39ec5de (add snow surface temp parameterization)
+=#

src/standalone/Snow/Snow.jl

@@ -12,7 +12,7 @@ using ClimaLand:
     net_radiation,
     AbstractModel,
     heaviside
-
+using SurfaceFluxes
 import ClimaLand:
     AbstractBC,
     make_update_aux,
@@ -218,6 +218,7 @@ auxiliary_vars(::SnowModel) = (
     :applied_energy_flux,
     :applied_water_flux,
     :snow_cover_fraction,
+    :ΔF,
 )
 
 auxiliary_types(::SnowModel{FT}) where {FT} = (
@@ -235,6 +236,7 @@ auxiliary_types(::SnowModel{FT}) where {FT} = (
     FT,
     FT,
     FT,
+    FT,
 )
 
 auxiliary_domain_names(::SnowModel) = (
@@ -270,7 +272,27 @@ function ClimaLand.make_update_aux(model::SnowModel{FT}) where {FT}
         @. p.snow.T =
             snow_bulk_temperature(Y.snow.U, Y.snow.S, p.snow.q_l, parameters)
 
-        @. p.snow.T_sfc = snow_surface_temperature(p.snow.T)
+        # original Tsfc code (Ts = Tbulk):
+        # @. p.snow.T_sfc = snow_surface_temperature_bulk(p.snow.T)
+        output =
+            snow_surface_temperature.(
+                p.drivers.u,
+                p.drivers.T,
+                p.drivers.P,
+                parameters.z_0m,
+                parameters.z_0b,
+                p.drivers.q,
+                model.atmos.h,
+                Y.snow.S,
+                p.snow.T,
+                p.drivers.thermal_state.ρ,
+                p.snow.energy_runoff,
+                p.drivers.LW_d,
+                p.drivers.SW_d,
+                parameters,
+            )
+        @. p.snow.T_sfc = output.T_s
+        @. p.snow.ΔF = output.ΔF
 
         @. p.snow.water_runoff =
             compute_water_runoff(Y.snow.S, p.snow.q_l, p.snow.T, parameters)

src/standalone/Snow/boundary_fluxes.jl

@@ -125,5 +125,5 @@ function snow_boundary_fluxes!(
             p.snow.turbulent_fluxes.lhf +
             p.snow.turbulent_fluxes.shf +
             p.snow.R_n - p.snow.energy_runoff
-        ) * p.snow.snow_cover_fraction
+         + p.snow.ΔF) * p.snow.snow_cover_fraction
 end