diff --git a/Project.toml b/Project.toml index 0212668476..8fc0834e77 100644 --- a/Project.toml +++ b/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" diff --git a/experiments/standalone/Snow/snowmip_simulation.jl b/experiments/standalone/Snow/snowmip_simulation.jl index 9fd7705a6b..82aea60c19 100644 --- a/experiments/standalone/Snow/snowmip_simulation.jl +++ b/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) +=# diff --git a/src/standalone/Snow/Snow.jl b/src/standalone/Snow/Snow.jl index 5d449a9378..30e589866c 100644 --- a/src/standalone/Snow/Snow.jl +++ b/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) diff --git a/src/standalone/Snow/boundary_fluxes.jl b/src/standalone/Snow/boundary_fluxes.jl index b97b608d9d..6e5cc0f7ad 100644 --- a/src/standalone/Snow/boundary_fluxes.jl +++ b/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 diff --git a/src/standalone/Snow/snow_parameterizations.jl b/src/standalone/Snow/snow_parameterizations.jl index bdbf0f5292..3f95f8dbfa 100644 --- a/src/standalone/Snow/snow_parameterizations.jl +++ b/src/standalone/Snow/snow_parameterizations.jl @@ -1,3 +1,8 @@ +using SurfaceFluxes +using StaticArrays +import SurfaceFluxes.Parameters as SFP +using RootSolvers + export snow_surface_temperature, snow_depth, specific_heat_capacity, @@ -62,10 +67,206 @@ end a helper function which returns the surface temperature for the snow model, which is stored in the aux state. """ -function ClimaLand.surface_temperature(model::SnowModel, Y, p, t) +function ClimaLand.surface_temperature(model::SnowModel{FT}, Y, p, t) where {FT} + # return max.(p.snow.T_sfc, FT(250)) return p.snow.T_sfc end +""" + partial_q_sat_partial_T_ice(P::FT, T::FT) where {FT} +Computes the quantity ∂q_sat∂T at temperature T and pressure P, +over ice. The temperature must be in Celsius. +Uses the polynomial approximation from Flatau et al. (1992). +""" +function partial_q_sat_partial_T_ice(P::FT, T::FT) where {FT} + T_celsius = T + esat = FT( + 6.11123516e2 + + 5.03109514e1 .* T_celsius + + 1.88369801 * T_celsius^2 + + 4.20547422e-2 * T_celsius^3 + + 6.14396778e-4 * T_celsius^4 + + 6.02780717e-6 * T_celsius^5 + + 3.87940929e-8 * T_celsius^6 + + 1.49436277e-10 * T_celsius^7 + + 2.62655803e-13 * T_celsius^8, + ) + desatdT = FT( + 5.03277922e1 + + 3.77289173 * T_celsius + + 1.26801703e-1 * T_celsius^2 + + 2.49468427e-3 * T_celsius^3 + + 3.13703411e-5 * T_celsius^4 + + 2.57180651e-7 * T_celsius^5 + + 1.33268878e-9 * T_celsius^6 + + 3.94116744e-12 * T_celsius^7 + + 4.98070196e-15 * T_celsius^8, + ) + + return FT(0.622) * P / (P - FT(0.378) * esat)^2 * desatdT +end + +function get_qsfc(thermo_params, T_s::FT, ρ_sfc::FT) where {FT} + if T_s <= 273.15 + q_sat = + Thermodynamics.q_vap_saturation_generic.( + Ref(thermo_params), + T_s, + ρ_sfc, + Ref(Thermodynamics.Ice()), + ) + else + q_sat = + Thermodynamics.q_vap_saturation_generic.( + Ref(thermo_params), + T_s, + ρ_sfc, + Ref(Thermodynamics.Liquid()), + ) + end + return q_sat +end + +function update_conditions( + T_s::FT, + thermo_params, + ρ_sfc::FT, + state_air, + z_0m::FT, + z_0b::FT, + surface_flux_params, +) where {FT} + q_sfc = get_qsfc(thermo_params, T_s, ρ_sfc) + thermal_state_sfc = + Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_s, q_sfc) + state_sfc = SurfaceFluxes.StateValues( + FT(0), + SVector{2, FT}(0, 0), + thermal_state_sfc, + ) + sc = + SurfaceFluxes.ValuesOnly(state_air, state_sfc, z_0m, z_0b, beta = FT(1)) + + conditions = SurfaceFluxes.surface_conditions( + surface_flux_params, + sc; + tol_neutral = SFP.cp_d(surface_flux_params) / 100000, + ) + + E0 = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch) + r_ae = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc)) + return E0, r_ae +end + +function solve_Ts(fns, T_initial::FT)::FT where {FT} + sol = RootSolvers.find_zero( + fns, + RootSolvers.NewtonsMethod(T_initial), + CompactSolution(), + nothing, + 10, + ) + T_s = sol.root + return T_s +end + +function calc_U(ρ_l::FT, d::FT, c::FT, T::FT, T0::FT)::FT where {FT} # assuming q_l = 1 + return ρ_l * d * c * (T - T0) +end + +function snow_surface_temperature( + u_air::FT, + T_air::FT, + P_air::FT, + z_0m::FT, + z_0b::FT, + q_air::FT, + h_air::FT, + SWE::FT, + T_bulk::FT, + ρ_sfc::FT, + energy_runoff::FT, + LW_d::FT, + SW_d::FT, + parameters, +) where {FT} + earth_param_set = parameters.earth_param_set + _LH_v0::FT = LP.LH_v0(earth_param_set) + _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set) + ρ_snow::FT = parameters.ρ_snow + h_sfc::FT = snow_depth(SWE, ρ_snow, _ρ_liq) + d = FT(0.2) + if h_sfc < d + return (T_s = T_bulk, ΔF = FT(0)) + end + thermo_params = LP.thermodynamic_parameters(earth_param_set) + thermal_state_air = + Thermodynamics.PhaseEquil_pTq(thermo_params, P_air, T_air, q_air) + cp_m::FT = Thermodynamics.cp_m(thermo_params, thermal_state_air) + ϵ_sfc::FT = parameters.ϵ_snow + α_sfc::FT = parameters.α_snow + d_sfc = FT(0) + + state_air = SurfaceFluxes.StateValues( + h_air - d_sfc - h_sfc, + SVector{2, FT}(u_air, 0), + thermal_state_air, + ) + surface_flux_params = LP.surface_fluxes_parameters(earth_param_set) + function conditions_func(T_s) + update_conditions( + T_s, + thermo_params, + ρ_sfc, + state_air, + z_0m, + z_0b, + surface_flux_params, + ) + end + E0(T_s) = conditions_func(T_s)[1] + r_ae(T_s) = conditions_func(T_s)[2] + LH(T_s) = _LH_v0 * E0(T_s) + + ρ_air = Thermodynamics.air_density(thermo_params, thermal_state_air) + SH(T_s) = cp_m * -ρ_air * (T_air - T_s) / r_ae(T_s) + + _σ = LP.Stefan(earth_param_set) + F_R(T_s) = -(1 - α_sfc) * SW_d - ϵ_sfc * (LW_d - _σ * T_s^4) + + F_h(T_s) = SH(T_s) + LH(T_s) + F_R(T_s) + energy_runoff + # κ = parameters.κ_ice + κ = snow_thermal_conductivity(ρ_snow, parameters) + f(T_s) = κ * (T_s - T_bulk) / (h_sfc / 2) + F_h(T_s) + + # derivatives for Newton's method + SH_prime(T_s) = cp_m * ρ_air / r_ae(T_s) + + ∂LHF∂qc(T_s) = ρ_air * _LH_v0 / r_ae(T_s) + ∂qc∂T(T_s) = partial_q_sat_partial_T_ice(P_air, T_s) + LH_prime(T_s) = ∂LHF∂qc(T_s) * ∂qc∂T(T_s) + + F_R_prime(T_s) = 4 * ϵ_sfc * _σ * T_s^3 + + F_h_prime(T_s) = SH_prime(T_s) + LH_prime(T_s) + F_R_prime(T_s) + f_prime(T_s) = (2 * κ / h_sfc) + F_h_prime(T_s) + + fns(T_s) = [f(T_s), f_prime(T_s)] + T_s::FT = solve_Ts(fns, T_bulk) + + # adjust T_s and internal energy accordingly if T_s > freezing temperature + if T_s > 273.15 + T_s_new = FT(273.15) + _T0::FT = FT(LP.T_0(earth_param_set)) + U_Ts::FT = calc_U(_ρ_liq, d, cp_m, T_s, _T0) + U_Tf::FT = calc_U(_ρ_liq, d, cp_m, T_s_new, _T0) + ΔU::FT = U_Ts - U_Tf + ΔF::FT = ΔU / parameters.Δt + return (T_s = T_s_new, ΔF = ΔF) + else + return (T_s = T_s, ΔF = FT(0)) + end +end """ ClimaLand.surface_specific_humidity(model::BucketModel, Y, p) @@ -109,7 +310,7 @@ end Returns the snow surface temperature assuming it is the same as the bulk temperature T. """ -snow_surface_temperature(T::FT) where {FT} = T +snow_surface_temperature_bulk(T::FT) where {FT} = T """