From 9ae26d7ed02711761c746919c500c65d3f11e6a6 Mon Sep 17 00:00:00 2001 From: Sarah Date: Wed, 9 Oct 2024 13:59:37 -0700 Subject: [PATCH 1/2] add snow surface temp parameterization --- Project.toml | 1 + .../standalone/Snow/snowmip_simulation.jl | 53 ++++- src/standalone/Snow/Snow.jl | 26 ++- src/standalone/Snow/boundary_fluxes.jl | 2 +- src/standalone/Snow/snow_parameterizations.jl | 205 +++++++++++++++++- 5 files changed, 278 insertions(+), 9 deletions(-) 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 """ From 57b892057e3815e2ee9b21436ff1b9a6269b6347 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Tue, 15 Oct 2024 17:06:13 -0700 Subject: [PATCH 2/2] runs but looks poor --- .buildkite/Manifest.toml | 2 +- .../standalone/Snow/snowmip_simulation.jl | 6 +- src/shared_utilities/drivers.jl | 1 - src/standalone/Snow/Snow.jl | 23 +- src/standalone/Snow/boundary_fluxes.jl | 43 +++- src/standalone/Snow/snow_parameterizations.jl | 206 +++++++++++++++++- 6 files changed, 247 insertions(+), 34 deletions(-) diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml index 1b6dba9bab..fcafaa06e1 100644 --- a/.buildkite/Manifest.toml +++ b/.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" diff --git a/experiments/standalone/Snow/snowmip_simulation.jl b/experiments/standalone/Snow/snowmip_simulation.jl index 82aea60c19..ff731d97eb 100644 --- a/experiments/standalone/Snow/snowmip_simulation.jl +++ b/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), ), ); @@ -245,7 +245,7 @@ ax3 = CairoMakie.Axis( ylabel = "Average Snow Temperature", xlabel = "Day of Year", ) -xlims!(ax3, 245, 281) +ylims!(ax3, 245, 281) CairoMakie.lines!( ax3, mean_obs_df.doy, diff --git a/src/shared_utilities/drivers.jl b/src/shared_utilities/drivers.jl index 6aeb8cf6d8..45cf0abf76 100644 --- a/src/shared_utilities/drivers.jl +++ b/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, diff --git a/src/standalone/Snow/Snow.jl b/src/standalone/Snow/Snow.jl index 30e589866c..59f435e033 100644 --- a/src/standalone/Snow/Snow.jl +++ b/src/standalone/Snow/Snow.jl @@ -254,6 +254,7 @@ auxiliary_domain_names(::SnowModel) = ( :surface, :surface, :surface, + :surface, ) @@ -272,28 +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) - # 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 6e5cc0f7ad..0f88e25d84 100644 --- a/src/standalone/Snow/boundary_fluxes.jl +++ b/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.ΔF) * p.snow.snow_cover_fraction + 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 3f95f8dbfa..d6d50cbfec 100644 --- a/src/standalone/Snow/snow_parameterizations.jl +++ b/src/standalone/Snow/snow_parameterizations.jl @@ -68,7 +68,6 @@ 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{FT}, Y, p, t) where {FT} - # return max.(p.snow.T_sfc, FT(250)) return p.snow.T_sfc end @@ -174,6 +173,201 @@ function calc_U(ρ_l::FT, d::FT, c::FT, T::FT, T0::FT)::FT where {FT} # assuming return ρ_l * d * c * (T - T0) end +function snow_turbulent_fluxes_at_a_point( + T_sfc::FT, + h_sfc::FT, + d_sfc::FT, + thermal_state_air, + u::FT, + h::FT, + gustiness::FT, + z_0m::FT, + z_0b::FT, + earth_param_set::EP, +) where {FT <: AbstractFloat, EP} + thermo_params = LP.thermodynamic_parameters(earth_param_set) + surface_flux_params = LP.surface_fluxes_parameters(earth_param_set) + _T_freeze = LP.T_freeze(earth_param_set) + ρ_sfc = ClimaLand.compute_ρ_sfc(thermo_params, thermal_state_air, T_sfc) + qsat_over_ice = Thermodynamics.q_vap_saturation_generic( + thermo_params, + T_sfc, + ρ_sfc, + Thermodynamics.Ice(), + ) + qsat_over_liq = Thermodynamics.q_vap_saturation_generic( + thermo_params, + T_sfc, + ρ_sfc, + Thermodynamics.Liquid(), + ) + q_sfc = T_sfc >= _T_freeze ? qsat_over_liq : qsat_over_ice + thermal_state_sfc = + Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sfc) + ρ_air = Thermodynamics.air_density(thermo_params, thermal_state_air) + _LH_v0::FT = LP.LH_v0(earth_param_set) + _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set) + cp_m_sfc::FT = Thermodynamics.cp_m(thermo_params, thermal_state_sfc) + T_air = Thermodynamics.air_temperature(thermo_params, thermal_state_air) + P_sfc = Thermodynamics.air_pressure(thermo_params, thermal_state_sfc) + Rm_air = Thermodynamics.gas_constant_air(thermo_params, thermal_state_air) + + # SurfaceFluxes.jl expects a relative difference between where u = 0 + # and the atmosphere height. Here, we assume h and h_sfc are measured + # relative to a common reference. Then d_sfc + h_sfc + z_0m is the apparent + # source of momentum, and + # Δh ≈ h - d_sfc - h_sfc is the relative height difference between the + # apparent source of momentum and the atmosphere height. + + # In this we have neglected z_0m and z_0b (i.e. assumed they are small + # compared to Δh). + state_sfc = SurfaceFluxes.StateValues( + FT(0), + SVector{2, FT}(0, 0), + thermal_state_sfc, + ) + state_in = SurfaceFluxes.StateValues( + h - d_sfc - h_sfc, + SVector{2, FT}(u, 0), + thermal_state_air, + ) + + # State containers + states = SurfaceFluxes.ValuesOnly( + state_in, + state_sfc, + z_0m, + z_0b, + gustiness = gustiness, + ) + scheme = SurfaceFluxes.PointValueScheme() + conditions = + SurfaceFluxes.surface_conditions(surface_flux_params, states, scheme) + + # aerodynamic resistance + r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states)) + + # latent heat flux + E0::FT = + SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch) # mass flux at potential evaporation rate + LH = _LH_v0 * E0 # Latent heat flux + + # sensible heat flux + SH = SurfaceFluxes.sensible_heat_flux( + surface_flux_params, + conditions.Ch, + states, + scheme, + ) + + # vapor flux in volume of liquid water with density 1000kg/m^3 + Ẽ = E0 / _ρ_liq + + # Derivatives + # We ignore ∂r_ae/∂T_sfc, ∂u*/∂T_sfc + ∂ρsfc∂T = + ρ_air * + (Thermodynamics.cv_m(thermo_params, thermal_state_air) / Rm_air) * + ( + T_sfc / T_air + )^(Thermodynamics.cv_m(thermo_params, thermal_state_air) / Rm_air - 1) / + T_air + ∂cp_m_sfc∂T = 0 # Possibly can address at a later date + + ∂LHF∂q = ρ_sfc * _LH_v0 / r_ae + LH / ρ_sfc * ∂ρsfc∂T + + ∂SHF∂T = + ρ_sfc * cp_m_sfc / r_ae + + SH / ρ_sfc * ∂ρsfc∂T + + SH / cp_m_sfc * ∂cp_m_sfc∂T + + ∂q∂T = ClimaLand.Snow.partial_q_sat_partial_T_ice(P_sfc, T_sfc - _T_freeze) + + return ( + lhf = LH, + shf = SH, + vapor_flux = Ẽ, + ∂LHF∂T = ∂LHF∂q * ∂q∂T, + ∂SHF∂T = ∂SHF∂T, + ) +end + +function kat_snow_surface_temperature( + u_air::FT, + thermal_state_air::Thermodynamics.PhaseEquil, + h_air::FT, + SWE::FT, + T_bulk::FT, + LW_d::FT, + SW_d::FT, + q_l::FT, + t, + parameters, +) where {FT} + (; z_0m, z_0b, earth_param_set) = parameters + thermo_params = LP.thermodynamic_parameters(earth_param_set) + _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set) + _σ = LP.Stefan(earth_param_set) + _T_freeze = LP.T_freeze(earth_param_set) + _T_ref::FT = FT(LP.T_0(earth_param_set)) + _LH_f0::FT = LP.LH_f0(earth_param_set) + + T_air::FT = Thermodynamics.air_temperature(thermo_params, thermal_state_air) + ϵ_sfc::FT = parameters.ϵ_snow + α_sfc::FT = parameters.α_snow + d_sfc = FT(0) + + ρ_snow::FT = parameters.ρ_snow + h_sfc::FT = snow_depth(SWE, ρ_snow, _ρ_liq) + d = FT(0.1) + if h_sfc < d + return (T_s = T_bulk, ΔF = FT(0)) + end + gustiness = FT(1) + + function sum_fluxes(T_sfc) + turb_fluxes = snow_turbulent_fluxes_at_a_point( + T_sfc, + FT(0),#h_sfc, + d_sfc, + thermal_state_air, + u_air, + h_air, + gustiness, + z_0m, + z_0b, + earth_param_set, + ) + F_R = -(1 - α_sfc) * SW_d - ϵ_sfc * (LW_d - _σ * T_sfc^4) # positive up + F_h = turb_fluxes.shf + turb_fluxes.lhf + F_R # positive up + κ = ClimaLand.Snow.snow_thermal_conductivity( + parameters.ρ_snow, + parameters, + ) + f_total = κ * (T_sfc - T_bulk) / (h_sfc / 2) + F_h # Fh = - k(T_sfc - T_bulk)/(h/2) + + # Derivatives + F_R_prime = 4 * ϵ_sfc * _σ * T_sfc^3 + F_h_prime = turb_fluxes.∂SHF∂T + turb_fluxes.∂LHF∂T + F_R_prime + f_total_prime = (2 * κ / h_sfc) + F_h_prime + return (f_total, f_total_prime) + end + T_s = solve_Ts(sum_fluxes, T_air) + # adjust T_s and internal energy accordingly if T_s > freezing temperature + if T_s > _T_freeze + c_sfc = specific_heat_capacity(FT(1), parameters) + c_bulk = specific_heat_capacity(q_l, parameters) + U_Ts::FT = _ρ_liq * d * c_sfc * (T_s - _T_ref)#calc_U(_ρ_liq, d, c_sfc, T_s, _T0) + U_Tf::FT = + _ρ_liq * d * c_bulk * (_T_freeze - _T_ref) - (1 - q_l) * _LH_f0# calc_U(_ρ_liq, d, c_sfc, _T_freeze, _T0) + ΔU::FT = U_Tf - U_Ts + ΔF::FT = ΔU / parameters.Δt + @show(ΔF) + return (T_s = _T_freeze, ΔF = ΔF) + else + return (T_s = T_s, ΔF = FT(0)) + end +end function snow_surface_temperature( u_air::FT, T_air::FT, @@ -195,7 +389,7 @@ function snow_surface_temperature( _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set) ρ_snow::FT = parameters.ρ_snow h_sfc::FT = snow_depth(SWE, ρ_snow, _ρ_liq) - d = FT(0.2) + d = FT(0.1) if h_sfc < d return (T_s = T_bulk, ΔF = FT(0)) end @@ -263,11 +457,14 @@ function snow_surface_temperature( ΔU::FT = U_Ts - U_Tf ΔF::FT = ΔU / parameters.Δt return (T_s = T_s_new, ΔF = ΔF) + elseif T_s < 250.0 + return (T_s = FT(250), ΔF = FT(0)) else return (T_s = T_s, ΔF = FT(0)) end end + """ ClimaLand.surface_specific_humidity(model::BucketModel, Y, p) @@ -299,9 +496,10 @@ function ClimaLand.surface_specific_humidity( Ref(Thermodynamics.Liquid()), ) q_l = p.snow.q_l - return @. qsat_over_ice * (1 - q_l) + q_l * (qsat_over_liq) + helper.(T_sfc, qsat_over_liq, qsat_over_ice) + #return @. qsat_over_ice * (1 - q_l) + q_l * (qsat_over_liq) end - +helper(x, o1, o2) = x >= 273.15 ? o1 : o2 """