Snow surface temperature parameterization

.buildkite/Manifest.toml

@@ -393,7 +393,7 @@ uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 version = "0.2.8"
 
 [[deps.ClimaLand]]
-deps = ["ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"]
+deps = ["ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "RootSolvers", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 version = "0.15.2"

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

@@ -22,7 +22,7 @@ using DelimitedFiles
 # Site-specific quantities
 # Error if no site argument is provided
 if length(ARGS) < 1
-    @error("Please provide a site name as command line argument")
+    SITE_NAME = "cdp"#@error("Please provide a site name as command line argument")
 else
     SITE_NAME = ARGS[1]
 end
@@ -67,7 +67,7 @@ timestepper = CTS.ARS111()
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 3,
         update_j = CTS.UpdateEvery(CTS.NewTimeStep),
     ),
 );
@@ -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",
 )
+ylims!(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/shared_utilities/drivers.jl

@@ -262,7 +262,6 @@ function turbulent_fluxes(
     d_sfc = displacement_height(model, Y, p)
     u_air = p.drivers.u
     h_air = atmos.h
-
     return turbulent_fluxes_at_a_point.(
         T_sfc,
         q_sfc,

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) = (
@@ -252,6 +254,7 @@ auxiliary_domain_names(::SnowModel) = (
     :surface,
     :surface,
     :surface,
+    :surface,
 )
 
 
@@ -270,8 +273,6 @@ 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)
-
         @. p.snow.water_runoff =
             compute_water_runoff(Y.snow.S, p.snow.q_l, p.snow.T, parameters)
 

src/standalone/Snow/boundary_fluxes.jl

@@ -101,6 +101,43 @@ function snow_boundary_fluxes!(
     t,
 ) where {FT}
     bc = model.boundary_conditions
+    # original Tsfc code (Ts = Tbulk):
+    # @. p.snow.T_sfc = snow_surface_temperature_bulk(p.snow.T)
+    parameters = model.parameters
+    #  output =
+    #       snow_surface_temperature.(
+    #           p.drivers.u,
+    #           p.drivers.T,
+    #           p.drivers.P,
+    #           parameters.z_0m,
+    #          parameters.z_0b,
+    #         p.drivers.q,
+    #   bc.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,
+    #)
+
+    output =
+        kat_snow_surface_temperature.(
+            p.drivers.u,
+            p.drivers.thermal_state,
+            bc.atmos.h,
+            Y.snow.S,
+            p.snow.T,
+            p.drivers.LW_d,
+            p.drivers.SW_d,
+            p.snow.q_l,
+            t,
+            parameters,
+        )
+
+    @. p.snow.T_sfc = output.T_s
+    @. p.snow.ΔF = output.ΔF
     p.snow.turbulent_fluxes .= turbulent_fluxes(bc.atmos, model, Y, p, t)
     p.snow.R_n .= net_radiation(bc.radiation, model, Y, p, t)
     # How does rain affect the below?
@@ -117,13 +154,13 @@ function snow_boundary_fluxes!(
     _LH_f0 = FT(LP.LH_f0(model.parameters.earth_param_set))
     _ρ_liq = FT(LP.ρ_cloud_liq(model.parameters.earth_param_set))
     ρe_falling_snow = -_LH_f0 * _ρ_liq # per unit vol of liquid water
-
     # positive fluxes are TOWARDS atmos
     @. p.snow.total_energy_flux =
         P_snow * ρe_falling_snow +
         (
             p.snow.turbulent_fluxes.lhf +
             p.snow.turbulent_fluxes.shf +
-            p.snow.R_n - p.snow.energy_runoff
+            p.snow.R_n - p.snow.energy_runoff#) * p.snow.snow_cover_fraction
+            + p.snow.ΔF
         ) * p.snow.snow_cover_fraction
 end