use SF.jl sensible heat flux

Author: juliasloan25

src/shared_utilities/drivers.jl

@@ -304,6 +304,11 @@ Computes turbulent surface fluxes at a point on a surface given
 
 This returns an energy flux and a liquid water volume flux, stored in
 a tuple with self explanatory keys.
+
+Please note that this function, if r_sfc is set to zero, makes no alteration
+to the output of the `SurfaceFluxes.surface_conditions` call, aside
+from unit conversion of evaporation from a mass to a liquid water volume
+flux.
 """
 function turbulent_fluxes_at_a_point(
     T_sfc::FT,
@@ -341,7 +346,7 @@ function turbulent_fluxes_at_a_point(
     )
 
     # State containers
-    sc = SurfaceFluxes.ValuesOnly(
+    states = SurfaceFluxes.ValuesOnly(
         state_in,
         state_sfc,
         z_0m,
@@ -350,26 +355,33 @@ function turbulent_fluxes_at_a_point(
         gustiness = gustiness,
     )
     surface_flux_params = LP.surface_fluxes_parameters(earth_param_set)
-    conditions = SurfaceFluxes.surface_conditions(
-        surface_flux_params,
-        sc;
-        tol_neutral = SFP.cp_d(surface_flux_params) / 100000,
-    )
+    scheme = SurfaceFluxes.PointValueScheme()
+    conditions =
+        SurfaceFluxes.surface_conditions(surface_flux_params, states, scheme)
     _LH_v0::FT = LP.LH_v0(earth_param_set)
     _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set)
 
-    cp_m::FT = Thermodynamics.cp_m(thermo_params, ts_in)
-    T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in)
-    ΔT = T_in - T_sfc
-    r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc))
-    ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_in)
+    # 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
+    E = E0 * r_ae / (r_sfc + r_ae) # mass flux accounting for additional surface resistance
+    LH = _LH_v0 * E # Latent heat flux
 
-    E0::FT = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch)
-    E = E0 * r_ae / (r_sfc + r_ae)
+    # 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
     Ẽ = E / _ρ_liq
-    H = -ρ_air * cp_m * ΔT / r_ae
-    LH = _LH_v0 * E
-    return (lhf = LH, shf = H, vapor_flux = Ẽ, r_ae = r_ae)
+
+    return (lhf = LH, shf = SH, vapor_flux = Ẽ, r_ae = r_ae)
 end
 
 """

src/standalone/Soil/energy_hydrology.jl

@@ -859,6 +859,17 @@ Computes turbulent surface fluxes for soil at a point on a surface given
 
 This returns an energy flux and a liquid water volume flux, stored in
 a tuple with self explanatory keys.
+
+Notes:
+In this function, we call SF twice - once for ice, and once for water. 
+We then combine the fluxes from them to get the total SH, LH. 
+The vapor fluxes are applied to ice and water differently - as either a source term
+or a boundary condition-  so they are kept separate.
+
+For ice, we supply a beta factor and do not need an additional surface resistance, 
+and so we use the output of surface_conditions directly. 
+For water, we do, and it's a little complicated how we handle it. 
+But once we correct E, we compute SH and LH the same as SurfaceFluxes.jl does.
 """
 function soil_turbulent_fluxes_at_a_point(
     T_sfc::FT,
@@ -881,24 +892,31 @@ function soil_turbulent_fluxes_at_a_point(
     γT_ref::FT,
     earth_param_set::P,
 ) where {FT <: AbstractFloat, P}
+    # Parameters
+    surface_flux_params = LP.surface_fluxes_parameters(earth_param_set)
     thermo_params = LP.thermodynamic_parameters(earth_param_set)
+    _LH_v0::FT = LP.LH_v0(earth_param_set)
+    _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set)
+
+    # Atmos air state
+    state_air = SurfaceFluxes.StateValues(
+        h_air - d_sfc - h_sfc,
+        SVector{2, FT}(u_air, 0),
+        thermal_state_air,
+    )
+    q_air::FT =
+        Thermodynamics.total_specific_humidity(thermo_params, thermal_state_air)
+
     # Estimate surface air density
     ρ_sfc::FT = ClimaLand.compute_ρ_sfc(thermo_params, thermal_state_air, T_sfc)
-    # Compute saturated specific humidities at surface over ice and liquid water
-    q_sat_ice::FT = Thermodynamics.q_vap_saturation_generic(
-        thermo_params,
-        T_sfc,
-        ρ_sfc,
-        Thermodynamics.Ice(),
-    )
+
+    # Now compute surface fluxes from liquid water component:
     q_sat_liq::FT = Thermodynamics.q_vap_saturation_generic(
         thermo_params,
         T_sfc,
         ρ_sfc,
         Thermodynamics.Liquid(),
     )
-
-    # Evaporation:
     thermal_state_sfc::Thermodynamics.PhaseEquil{FT} =
         Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sat_liq)# use to get potential evaporation E0, r_ae (weak dependence on q).
 
@@ -916,12 +934,6 @@ function soil_turbulent_fluxes_at_a_point(
         SVector{2, FT}(0, 0),
         thermal_state_sfc,
     )
-    state_air = SurfaceFluxes.StateValues(
-        h_air - d_sfc - h_sfc,
-        SVector{2, FT}(u_air, 0),
-        thermal_state_air,
-    )
-
     # State containers
     states = SurfaceFluxes.ValuesOnly(
         state_air,
@@ -930,21 +942,9 @@ function soil_turbulent_fluxes_at_a_point(
         z_0b,
         gustiness = gustiness,
     )
-    surface_flux_params = LP.surface_fluxes_parameters(earth_param_set)
-    conditions = SurfaceFluxes.surface_conditions(
-        surface_flux_params,
-        states;
-        tol_neutral = SFP.cp_d(surface_flux_params) / 100000,
-    )
-    _LH_v0::FT = LP.LH_v0(earth_param_set)
-    _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set)
-    cp_m::FT = Thermodynamics.cp_m(thermo_params, thermal_state_air)
-    T_air::FT = Thermodynamics.air_temperature(thermo_params, thermal_state_air)
-    ρ_air::FT = Thermodynamics.air_density(thermo_params, thermal_state_air)
-    q_air::FT =
-        Thermodynamics.total_specific_humidity(thermo_params, thermal_state_air)
-
-    ΔT::FT = T_air - T_sfc
+    scheme = SurfaceFluxes.PointValueScheme()
+    conditions =
+        SurfaceFluxes.surface_conditions(surface_flux_params, states, scheme)
     r_ae_liq::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states))
 
     E0::FT =
@@ -970,9 +970,22 @@ function soil_turbulent_fluxes_at_a_point(
     else
         Ẽ *= 0 # condensation, set to zero
     end
-    H_liq::FT = -ρ_air * cp_m * ΔT / r_ae_liq
 
-    # Sublimation:
+    # sensible heat flux
+    SH_liq = SurfaceFluxes.sensible_heat_flux(
+        surface_flux_params,
+        conditions.Ch,
+        states,
+        scheme,
+    )
+
+    # Repeat for ice component:
+    q_sat_ice::FT = Thermodynamics.q_vap_saturation_generic(
+        thermo_params,
+        T_sfc,
+        ρ_sfc,
+        Thermodynamics.Ice(),
+    )
     thermal_state_sfc =
         Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sat_ice)
     β_i::FT = FT(1)
@@ -996,25 +1009,28 @@ function soil_turbulent_fluxes_at_a_point(
         beta = β_i,
         gustiness = gustiness,
     )
-    conditions = SurfaceFluxes.surface_conditions(
-        surface_flux_params,
-        states;
-        tol_neutral = SFP.cp_d(surface_flux_params) / 100000,
-    )
+    conditions =
+        SurfaceFluxes.surface_conditions(surface_flux_params, states, scheme)
     S::FT =
         SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch)# sublimation rate, mass flux
     S̃::FT = S / _ρ_liq
 
     r_ae_ice::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states))
-    H_ice::FT = -ρ_air * cp_m * ΔT / r_ae_ice
+    # sensible heat flux from ice
+    SH_ice = SurfaceFluxes.sensible_heat_flux(
+        surface_flux_params,
+        conditions.Ch,
+        states,
+        scheme,
+    )
 
-    # Heat fluxes
+    # Heat fluxes for soil
     LH::FT = _LH_v0 * (Ẽ + S̃) * _ρ_liq
-    H::FT = (H_ice + H_liq) / 2
+    SH::FT = (SH_ice + SH_liq) / 2
     r_ae::FT = 2 * r_ae_liq * r_ae_ice / (r_ae_liq + r_ae_ice) # implied by definition of H
     return (
         lhf = LH,
-        shf = H,
+        shf = SH,
         vapor_flux_liq = Ẽ,
         r_ae = r_ae,
         vapor_flux_ice = S̃,

src/standalone/Vegetation/canopy_boundary_fluxes.jl

@@ -286,6 +286,10 @@ end
 
 Computes the turbulent surface fluxes for the canopy at a point
 and returns the fluxes in a named tuple.
+
+Note that an additiontal resistance is used in computing both
+evaporation and sensible heat flux, and this modifies the output
+of `SurfaceFluxes.surface_conditions`.
 """
 function canopy_turbulent_fluxes_at_a_point(
     T_sfc::FT,
@@ -315,7 +319,7 @@ function canopy_turbulent_fluxes_at_a_point(
     )
 
     # State containers
-    sc = SurfaceFluxes.ValuesOnly(
+    states = SurfaceFluxes.ValuesOnly(
         state_in,
         state_sfc,
         z_0m,
@@ -324,37 +328,62 @@ function canopy_turbulent_fluxes_at_a_point(
         gustiness = gustiness,
     )
     surface_flux_params = LP.surface_fluxes_parameters(earth_param_set)
-    conditions = SurfaceFluxes.surface_conditions(
-        surface_flux_params,
-        sc;
-        tol_neutral = SFP.cp_d(surface_flux_params) / 100000,
-    )
+    scheme = SurfaceFluxes.PointValueScheme()
+    conditions =
+        SurfaceFluxes.surface_conditions(surface_flux_params, states, scheme)
     _LH_v0::FT = LP.LH_v0(earth_param_set)
     _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set)
-
-    cp_m::FT = Thermodynamics.cp_m(thermo_params, ts_in)
-    T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in)
-    ΔT = T_in - T_sfc
-    r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc))
-    ρ_sfc::FT = Thermodynamics.air_density(thermo_params, ts_sfc)
+    cp_m_sfc::FT = Thermodynamics.cp_m(thermo_params, ts_sfc)
+    r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states))
     ustar::FT = conditions.ustar
+    ρ_sfc = Thermodynamics.air_density(thermo_params, ts_sfc)
+    T_int = Thermodynamics.air_temperature(thermo_params, ts_in)
+    Rm_int = Thermodynamics.gas_constant_air(thermo_params, ts_in)
+    ρ_air = Thermodynamics.air_density(thermo_params, ts_in)
+
     r_b::FT = FT(1 / 0.01 * (ustar / 0.04)^(-1 / 2)) # CLM 5, tech note Equation 5.122
     leaf_r_b = r_b / LAI
     canopy_r_b = r_b / (LAI + SAI)
-    E0::FT = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch)
+
+    E0::FT =
+        SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch)
     E = E0 * r_ae / (leaf_r_b + leaf_r_stomata + r_ae) # CLM 5, tech note Equation 5.101, and fig 5.2b, assuming all sunlit, f_wet = 0
     Ẽ = E / _ρ_liq
-    H = -ρ_sfc * cp_m * ΔT / (canopy_r_b + r_ae) # CLM 5, tech note Equation 5.88, setting H_v = H and solving to remove T_s
+
+    SH =
+        SurfaceFluxes.sensible_heat_flux(
+            surface_flux_params,
+            conditions.Ch,
+            states,
+            scheme,
+        ) * r_ae / (canopy_r_b + r_ae)
+    # The above follows from CLM 5, tech note Equation 5.88, setting H_v = SH and solving to remove T_s, ignoring difference between cp in atmos and above canopy. 
     LH = _LH_v0 * E
-    # We ignore ∂r_ae/∂T_sfc, ∂u*/∂T_sfc
-    ∂LHF∂qc = ρ_sfc * _LH_v0 / (leaf_r_b + leaf_r_stomata + r_ae) # Note that our estimate of ρ_sfc depends on T_sfc, but we ignore the derivative here.
-    ∂SHF∂Tc = ρ_sfc * cp_m / (canopy_r_b + r_ae) # Same approximation re: ρ_sfc(T_sfc)
+
+    # Derivatives
+    # We ignore ∂r_ae/∂T_sfc, ∂u*/∂T_sfc, ∂r_stomata∂Tc
+    ∂ρsfc∂Tc =
+        ρ_air *
+        (Thermodynamics.cv_m(thermo_params, ts_in) / Rm_int) *
+        (T_sfc / T_int)^(Thermodynamics.cv_m(thermo_params, ts_in) / Rm_int - 1) /
+        T_int
+    ∂cp_m_sfc∂Tc = 0 # Possibly can address at a later date
+
+    ∂LHF∂qc =
+        ρ_sfc * _LH_v0 / (leaf_r_b + leaf_r_stomata + r_ae) +
+        LH / ρ_sfc * ∂ρsfc∂Tc
+
+    ∂SHF∂Tc =
+        ρ_sfc * cp_m_sfc / (canopy_r_b + r_ae) +
+        SH / ρ_sfc * ∂ρsfc∂Tc +
+        SH / cp_m_sfc * ∂cp_m_sfc∂Tc
     return (
         lhf = LH,
-        shf = H,
+        shf = SH,
         vapor_flux = Ẽ,
         r_ae = r_ae,
         ∂LHF∂qc = ∂LHF∂qc,
         ∂SHF∂Tc = ∂SHF∂Tc,
     )
+
 end

test/standalone/Vegetation/canopy_model.jl

@@ -1070,12 +1070,12 @@ end
         finitediff_SHF = (p_2.canopy.energy.shf .- p.canopy.energy.shf) ./ ΔT
         estimated_SHF = p.canopy.energy.∂SHF∂Tc
         @test parent(abs.(finitediff_SHF .- estimated_SHF) ./ finitediff_SHF)[1] <
-              0.05
+              0.15
 
         finitediff_LHF = (p_2.canopy.energy.lhf .- p.canopy.energy.lhf) ./ ΔT
         estimated_LHF = p.canopy.energy.∂LHF∂qc .* p.canopy.energy.∂qc∂Tc
         @test parent(abs.(finitediff_LHF .- estimated_LHF) ./ finitediff_LHF)[1] <
-              0.5
+              0.3
 
         # Error in `q` derivative is large
         finitediff_q = (q_sfc2 .- q_sfc) ./ ΔT