Skip to content
Snippets Groups Projects
Commit e3208bb3 authored by Jan Nitzbon's avatar Jan Nitzbon Committed by Brian Groenke
Browse files

Add new surface energy balance parameterizations

parent dbd55040
No related branches found
No related tags found
1 merge request!34Add new surface energy balance parameterizations
Hide whitespace changes
Inline Side-by-side
examples/heat_vgfc_seb_samoylov_custom.jl
+ 12
− 19
View file @ e3208bb3
......@@ -12,17 +12,9 @@ soilprofile = SoilProfile(
0.4u"m" => soilparameters(xic=0.30,por=0.55,sat=1.0,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55),
3.0u"m" => soilparameters(xic=0.0,por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50),
10.0u"m" => soilparameters(xic=0.0,por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30),
),
);
# mid-winter temperature profile
tempprofile = TemperatureProfile(
0.01u"m" => -15.9u"°C",
0.05u"m" => -15.49u"°C",
0.10u"m" => -15.32u"°C",
0.20u"m" => -14.44u"°C",
0.30u"m" => -14.18u"°C",
0.40u"m" => -13.50u"°C",
1000.0u"m" => 10.2u"°C",
)
tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C");
Tair = TimeSeriesForcing(ustrip.(forcings.data.Tair), forcings.timestamps, :Tair);
# assume other forcings don't (yet) have units
......@@ -36,22 +28,23 @@ tspan = (DateTime(2010,1,1), DateTime(2011,1,1))
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
strat = Stratigraphy(
-2.0u"m" => top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z)),
Tuple(z => subsurface(:soil1, Soil(para=para), Heat(:H, freezecurve=SFCC(DallAmico()))) for (z,para) in soilprofile),
-z*u"m" => top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z,solscheme=Analytical(),stabfun=Businger())),
Tuple(knot.depth => subsurface(Symbol(:soil,i), Soil(para=knot.value), Heat(:H, freezecurve=SFCC(DallAmico()))) for (i,knot) in enumerate(soilprofile)),
1000.0u"m" => bottom(GeothermalHeatFlux(0.053u"J/s/m^2")),
);
grid = Grid(gridvals);
model = Tile(strat, grid, initT);
tile = Tile(strat, grid, initT);
# define time span
tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
p = parameters(model)
u0, du0 = initialcondition!(model, tspan, p)
p = parameters(tile)
u0, du0 = initialcondition!(tile, tspan, p)
# CryoGrid front-end for ODEProblem
prob = CryoGridProblem(model,u0,tspan,p,savevars=(:T,))
prob = CryoGridProblem(tile,u0,tspan,p,step_limiter=nothing,savevars=(:T,))
# solve with forward Euler (w/ CFL) and construct CryoGridOutput from solution
out = @time solve(prob, Euler(), dt=2*60.0, saveat=24*3600.0, progress=true) |> CryoGridOutput;
# Plot it!
zs = [1:10...,20:10:100...]
zs = [1,5,10,15,20:10:100...]
cg = Plots.cgrad(:copper,rev=true);
plot(out.H[Z(zs)], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150)
plot(out.T[Z(zs)], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
plot(ustrip.(out.H[Z(zs)]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150)
plot(ustrip.(out.T[Z(zs)]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
src/Physics/SEB/SEB.jl
+ 64
− 16
View file @ e3208bb3
module SEB
using ..HeatConduction: Heat
using ..Physics
using ..Boundaries
using CryoGrid.InputOutput: Forcing
using CryoGrid.Physics
using CryoGrid.Physics.Boundaries
using CryoGrid.Physics.HeatConduction
using CryoGrid.Physics.Soils
using CryoGrid.Numerics
using CryoGrid.Utils
import CryoGrid: BoundaryProcess, BoundaryStyle, Neumann, Top
import CryoGrid: initialcondition!, variables
import CryoGrid: initialcondition!, variables, boundaryvalue
using Unitful
Base.@kwdef struct SEBParams <: IterableStruct
export SurfaceEnergyBalance, SEBParams
export Businger, HøgstrømSHEBA, Iterative, Analytical, Numerical
abstract type SolutionScheme end
"""
Byun 1990
"""
struct Analytical <: SolutionScheme end
"""
equation by Westermann2016, but use Newton solver
- not implemented yet
"""
struct Numerical <: SolutionScheme end
"""
Westermann2016, use info from last time step
"""
struct Iterative <: SolutionScheme end
abstract type StabilityFunctions end
"""
Businger 1971
"""
struct Businger <: StabilityFunctions end
"""
Høgstrøm 1988 (unstable conditions)
SHEBA, Uttal et al., 2002, Grachev et al. 2007 (stable conditions)
"""
struct HøgstrømSHEBA <: StabilityFunctions end
Base.@kwdef struct SEBParams{TSolution,TStabFun}
# surface properties --> should be associated with the Stratigraphy and maybe made state variables
α::Float"1" = 0.2xu"1" # surface albedo [-]
ϵ::Float"1" = 0.97xu"1" # surface emissivity [-]
......@@ -30,28 +63,43 @@ Base.@kwdef struct SEBParams <: IterableStruct
# material properties (assumed to be constant)
ρₐ::Float"kg/m^3" = 1.293xu"kg/m^3" # density of air at standard pressure and 0°C [kg/m^3]
cₐ::Float"J/(m^3*K)" = 1005.7xu"J/(kg*K)" * ρₐ # volumetric heat capacity of dry air at standard pressure and 0°C [J/(m^3*K)]
Pr₀::Float64 = 0.74 # turbulent Prandtl number
βₘ::Float64 = 4.7
βₕ::Float64 = βₘ/Pr₀
γₘ::Float64 = 15.0
γₕ::Float64 = 9.0
# type-dependent parameters
solscheme::TSolution = Iterative()
stabfun::TStabFun = HøgstrømSHEBA()
end
struct SurfaceEnergyBalance{F} <: BoundaryProcess{Heat}
forcing::F
sebparams::SEBParams
"""
SurfaceEnergyBalance{TSolution,TStabFun,F} <: BoundaryProcess{Heat}
Surface energy balance upper boundary condition.
"""
struct SurfaceEnergyBalance{TSolution,TStabFun,F} <: BoundaryProcess{Heat}
forcings::F
sebparams::SEBParams{TSolution,TStabFun}
SurfaceEnergyBalance(forcings::NamedTuple, sebparams::SEBParams) = new{typeof(sebparams.solscheme),typeof(sebparams.stabfun),typeof(forcings)}(forcings, sebparams)
function SurfaceEnergyBalance(
Tair::Forcing,
Tair::Forcing{Float"°C"},
p::Forcing,
q::Forcing,
wind::Forcing,
Lin::Forcing,
Sin::Forcing,
z::Float"m",
params::SEBParams=SEBParams()
z::Float"m";
kwargs...
)
forcing = (Tair = Tair, p = p, q = q, wind = wind, Lin = Lin, Sin = Sin, z = z);
new{typeof(forcing)}(forcing, params)
forcings = (Tair = Tair, p = p, q = q, wind = wind, Lin = Lin, Sin = Sin, z = z);
sebparams = SEBParams(;kwargs...);
new{typeof(sebparams.solscheme),typeof(sebparams.stabfun),typeof(forcings)}(forcings, sebparams)
end
end
BoundaryStyle(::Type{<:SurfaceEnergyBalance}) = Neumann()
include("seb_heat.jl")
end
\ No newline at end of file
end
src/Physics/SEB/seb_heat.jl
+ 128
− 36
View file @ e3208bb3
......@@ -20,33 +20,35 @@ function initialcondition!(top::Top, seb::SurfaceEnergyBalance, state)
@setscalar state.ustar = 10.;
end
BoundaryStyle(::Type{<:SurfaceEnergyBalance}) = Neumann()
"""
Top interaction, ground heat flux from surface energy balance. (no snow, no water body, no infiltration)
"""
function (seb::SurfaceEnergyBalance)(top::Top, soil::Soil, heat::Heat, stop, ssoil)
function boundaryvalue(seb::SurfaceEnergyBalance, ::Top, ::Heat, ::Soil, stop, ssoil)
# TODO (optimize): pre-compute all forcings at time t here, then pass to functions
# 1. calculate radiation budget
# outgoing shortwave radiation as reflected
@setscalar stop.Sout = let α = seb.sebparams.α, Sin = seb.forcing.Sin(stop.t);
-α * Sin # Eq. (2) in Westermann et al. (2016)
@setscalar stop.Sout = let α = seb.sebparams.α, Sin = seb.forcings.Sin(stop.t);
-α * Sin # Eq. (2) in Westermann et al. (2016)
end
# outgoing longwave radiation composed of emitted and reflected radiation
@setscalar stop.Lout = let ϵ = seb.sebparams.ϵ, σ = seb.sebparams.σ, T₀ = ssoil.T[1], Lin = seb.forcing.Lin(stop.t);
-ϵ * σ * T₀^4 - (1 - ϵ) * Lin # Eq. (3) in Westermann et al. (2016)
@setscalar stop.Lout = let ϵ = seb.sebparams.ϵ, σ = seb.sebparams.σ, T₀ = ssoil.T[1], Lin = seb.forcings.Lin(stop.t);
-ϵ * σ * normalize_temperature(T₀)^4 - (1 - ϵ) * Lin # Eq. (3) in Westermann et al. (2016)
end
# net radiation budget
@setscalar stop.Qnet = let Sin = seb.forcing.Sin(stop.t), Lin = seb.forcing.Lin(stop.t), Sout = getscalar(stop.Sout), Lout = getscalar(stop.Lout);
@setscalar stop.Qnet = let Sin = seb.forcings.Sin(stop.t), Lin = seb.forcings.Lin(stop.t), Sout = getscalar(stop.Sout), Lout = getscalar(stop.Lout);
Sin + Sout + Lin + Lout
end
# 2. calcuate turbulent heat flux budget
# determine atmospheric stability conditions
@setscalar stop.ustar = ustar(seb, stop);
@setscalar stop.Lstar = Lstar(seb, stop, ssoil);
@setscalar stop.ustar = ustar(seb, stop);
# sensible heat flux
@setscalar stop.Qh = Q_H(seb, stop, ssoil);
......@@ -66,15 +68,12 @@ end
"""
Density of air at given tempeature and pressure
"""
density_air(seb::SurfaceEnergyBalance,T::Real"°C",p::Real"Pa") = p / (T * seb.sebparams.Rₐ);
density_air(seb::SurfaceEnergyBalance,T::Real"°C",p::Real"Pa") = p / (normalize_temperature(T) * seb.sebparams.Rₐ);
"""
Saturation pressure of water/ice according to the empirical August-Roche-Magnus formula
Note: T is passed [K] and converted to [°C]
"""
estar(T::Real"°C") = ( (T > 0) ? 611.2 * exp(17.62 * T / (243.12 + T))# Eq. (B3) in Westermann et al. (2016)
: 611.2 * exp(22.46 * T / (272.62 + T)) ;
)
estar(T::Real"°C") = (T > 0) ? 611.2 * exp(17.62 * T / (243.12 + T)) : 611.2 * exp(22.46 * T / (272.62 + T)); # Eq. (B3) in Westermann et al. (2016)
"""
Latent heat of evaporation/condensation of water in [J/kg]
......@@ -93,30 +92,88 @@ Friction velocity according to Monin-Obukhov theory
"""
function ustar(seb::SurfaceEnergyBalance, stop)
let κ = seb.sebparams.κ,
uz = seb.forcing.wind(stop.t), # wind speed at height z
z = seb.forcing.z, # height z of wind forcing
uz = seb.forcings.wind(stop.t), # wind speed at height z
z = seb.forcings.z, # height z of wind forcing
z₀ = seb.sebparams.z₀, # aerodynamic roughness length [m]
Lstar = stop.Lstar |> getscalar;
κ * uz ./ (log(z / z₀) - Ψ_M(z / Lstar, z₀ / Lstar)) # Eq. (7) in Westermann et al. (2016)
κ * uz ./ (log(z / z₀) - Ψ_M(seb, z / Lstar, z₀ / Lstar)) # Eq. (7) in Westermann et al. (2016)
end
end
"""
Obukhov length according to Monin-Obukhov theory
Obukhov length according to Monin-Obukhov theory, iterative determination as in CryoGrid3 / Westermann et al. 2016
- uses the turubulent fluxes Qe and Qh as well as the friction velocity of the previous time step
"""
function Lstar(seb::SurfaceEnergyBalance, stop, ssoil)
function Lstar(seb::SurfaceEnergyBalance{Iterative}, stop, ssoil)
res = let κ = seb.sebparams.κ,
g = seb.sebparams.g,
Rₐ = seb.sebparams.Rₐ,
cₚ = seb.sebparams.cₐ / seb.sebparams.ρₐ, # specific heat capacity of air at constant pressure
Tₕ = seb.forcing.Tair(stop.t), # air temperature at height z over surface
p = seb.forcing.p(stop.t), # atmospheric pressure at surface
Tair = seb.forcings.Tair(stop.t),
Tₕ = normalize_temperature(Tair), # air temperature at height z over surface
p = seb.forcings.p(stop.t), # atmospheric pressure at surface
ustar = stop.ustar |> getscalar,
Qₑ = stop.Qe |> getscalar,
Qₕ = stop.Qh |> getscalar,
Llg = L_lg(ssoil.T[1]),
ρₐ = density_air(seb, Tₕ, p); # density of air at surface air temperature and surface pressure [kg/m^3]
-ρₐ * cₚ * Tₕ / (κ * g) * ustar^3 / (Qₕ + 0.61 * cₚ / Llg * Tₕ * Qₑ) # Eq. (8) in Westermann et al. (2016)
ρₐ = density_air(seb, Tair, p); # density of air at surface air temperature and surface pressure [kg/m^3]
-ρₐ * cₚ * Tₕ / (κ * g) * ustar^3 / (Qₕ + 0.61 * cₚ / Llg * Tₕ * Qₑ) # Eq. (8) in Westermann et al. (2016)
end
# upper and lower limits for Lstar
res = (abs(res) < 1e-7) ? sign(res) * 1e-7 : res
res = (abs(res) > 1e+7) ? sign(res) * 1e+7 : res
end
"""
Obukhov length, analytical solution of Monin-Obukhov theory according to Byun 1990
- implicitly assumes the Businger 1971 stability functions
- only uses the current surface temperature as input
"""
function Lstar(seb::SurfaceEnergyBalance{Analytical}, stop, ssoil)
res = let g = seb.sebparams.g,
Rₐ = seb.sebparams.Rₐ,
cₚ = seb.sebparams.cₐ / seb.sebparams.ρₐ, # specific heat capacity of air at constant pressure
Tₕ = normalize_temperature(seb.forcings.Tair(stop.t)), # air temperature at height z over surface
T₀ = normalize_temperature(ssoil.T[1]), # surface temperature
p = seb.forcings.p(stop.t), # atmospheric pressure at surface (height z)
p₀ = seb.forcings.p(stop.t), # normal pressure (for now assumed to be equal to p)
uz = seb.forcings.wind(stop.t), # wind speed at height z
z = seb.forcings.z, # height z of wind forcing
z₀ = seb.sebparams.z₀, # aerodynamic roughness length [m]
Pr₀ = seb.sebparams.Pr₀, # turbulent Prandtl number
γₕ = seb.sebparams.γₕ,
γₘ = seb.sebparams.γₘ,
βₕ = seb.sebparams.βₕ,
βₘ = seb.sebparams.βₘ;
Θₕ = Tₕ * (p₀/p)^(Rₐ/cₚ); # potential temperature (for now identical to actual temperature)
Θ₀ = T₀ * (p₀/p)^(Rₐ/cₚ);
# calcuate bulk Richardson number
Ri_b = g / Θ₀ * (Θₕ - Θ₀) * (z - z₀) / uz^2; # bulk Richardson number, eq. (9) in Byun 1990
# calulate ζ
a = (z / (z-z₀)) * log(z/z₀);
if Ri_b>0 # stable conditions
b = 2 * βₕ * (βₘ * Ri_b - 1);
c = -(2 * βₕ * Ri_b - 1) - (1 + (4 * (βₕ - βₘ) * Ri_b) / Pr₀ )^0.5
ζ = a / b * c; # eq. (19) in Byun 1990
else # unstable conditions
s_b = Ri_b / Pr₀; # eq. (30) in Byun 1990
Q_b = 1/9 * ( 1/γₘ^2 + 3 * γₕ/γₘ * s_b^2 ); # eq. (31) in Byun 1990
P_b = 1/54 * ( -2/γₘ^3 + 9/γₘ * (3 - γₕ/γₘ) * s_b^2 ); # eq. (32) in Byun 1990
θ_b = acos( P_b / Q_b^(3/2) ); # eq. (33) in Byun 1990
T_b = ( (P_b^2 - Q_b^3)^0.5 + abs(P_b) )^(1/3); # eq. (34) in Byun 1990
if Q_b^3-P_b^2 >= 0
ζ = a * (-2 * Q_b^0.5 * cos(θ_b/3) + 1/(3*γₘ) ) # eq. (28) in Byun 1990
else
ζ = a * (-(T_b + Q_b/T_b ) + 1/(3*γₘ) ); # eq. (29) in Byun 1990
end
end
# calculate Lstar
z / ζ;
end
# upper and lower limits for Lstar
res = (abs(res) < 1e-7) ? sign(res) * 1e-7 : res
......@@ -129,17 +186,17 @@ Sensible heat flux, defined as positive if it is a flux towards the surface
function Q_H(seb::SurfaceEnergyBalance, stop, ssoil)
let κ = seb.sebparams.κ,
Rₐ = seb.sebparams.Rₐ,
Tₕ = seb.forcing.Tair(stop.t), # air temperature
Tₕ = seb.forcings.Tair(stop.t), # air temperature
T₀ = ssoil.T[1], # surface temperature
cₚ = seb.sebparams.cₐ / seb.sebparams.ρₐ, # specific heat capacity of air at constant pressure
z = seb.forcing.z, # height at which forcing data are provided
z = seb.forcings.z, # height at which forcing data are provided
Lstar = stop.Lstar |> getscalar,
ustar = stop.ustar |> getscalar,
p = seb.forcing.p(stop.t),
p = seb.forcings.p(stop.t),
z₀ = seb.sebparams.z₀,
ρₐ = density_air(seb, Tₕ, p); # density of air at surface air temperature and surface pressure [kg/m^3]
rₐᴴ = (κ * ustar)^-1 * (log(z / z₀) - Ψ_HW(z / Lstar, z₀ / Lstar)) # Eq. (6) in Westermann et al. (2016)
rₐᴴ = (κ * ustar)^-1 * (log(z / z₀) - Ψ_HW(seb, z / Lstar, z₀ / Lstar)) # Eq. (6) in Westermann et al. (2016)
# calculate Q_H
-ρₐ * cₚ * (Tₕ - T₀) / rₐᴴ # Eq. (4) in Westermann et al. (2016)
......@@ -155,28 +212,28 @@ function Q_E(seb::SurfaceEnergyBalance, stop, ssoil)
let κ = seb.sebparams.κ,
γ = seb.sebparams.γ,
Rₐ = seb.sebparams.Rₐ,
Tₕ = seb.forcing.Tair(stop.t), # air temperature at height z over surface
Tₕ = seb.forcings.Tair(stop.t), # air temperature at height z over surface
T₀ = ssoil.T[1], # surface temperature
p = seb.forcing.p(stop.t), # atmospheric pressure at surface
qₕ = seb.forcing.q(stop.t), # specific humidity at height h over surface
z = seb.forcing.z, # height at which forcing data are provided
p = seb.forcings.p(stop.t), # atmospheric pressure at surface
qₕ = seb.forcings.q(stop.t), # specific humidity at height h over surface
z = seb.forcings.z, # height at which forcing data are provided
rₛ = seb.sebparams.rₛ, # surface resistance against evapotranspiration / sublimation [1/m]
Lstar = stop.Lstar |> getscalar,
ustar = stop.ustar |> getscalar,
Llg = L_lg(ssoil.T[1]),
Lsg = L_sg(ssoil.T[1]),
z₀ = seb.sebparams.z₀, # aerodynamic roughness length [m]
ρₐ = density_air(seb, seb.forcing.Tair(stop.t), seb.forcing.p(stop.t)); # density of air at surface air temperature and surface pressure [kg/m^3]
ρₐ = density_air(seb, Tₕ, seb.forcings.p(stop.t)); # density of air at surface air temperature and surface pressure [kg/m^3]
q₀ = γ * estar(T₀) / p # saturation pressure of water/ice at the surface; Eq. (B1) in Westermann et al (2016)
rₐᵂ = (κ * ustar)^-1 * (log(z / z₀) - Ψ_HW(z / Lstar, z₀ / Lstar)) # aerodynamic resistance Eq. (6) in Westermann et al. (2016)
L = (T₀ <= 0.0) ? Lsg : Llg # latent heat of sublimation/resublimation or evaporation/condensation [J/kg]
rₐᵂ = (κ * ustar)^-1 * (log(z / z₀) - Ψ_HW(seb, z / Lstar, z₀ / Lstar)) # aerodynamic resistance Eq. (6) in Westermann et al. (2016)
L = (T₀ <= 0.0) ? Lsg : Llg # latent heat of sublimation/resublimation or evaporation/condensation [J/kg]
# calculate Q_E
res = (qₕ > q₀) ? -ρₐ * L * (qₕ - q₀) / (rₐᵂ) : # Eq. (5) in Westermann et al. (2016) # condensation / deposition (no aerodynamics resistance)
-ρₐ * L * (qₕ - q₀) / (rₐᵂ + rₛ) ; # evaporation / sublimation (account for surface resistance against evapotranspiration/sublimation)
res = (T₀ <= 0.0) ? 0 : res ; # for now: set sublimation and deposition to zero.
res = (T₀ <= 0.0) ? 0 : res ; # for now: set sublimation and deposition to zero.
end
end
......@@ -185,7 +242,7 @@ Integrated stability function for heat/water transport
Högström, 1988
SHEBA, Uttal et al., 2002, Grachev et al. 2007
"""
function Ψ_HW(ζ₁::Real, ζ₂::Real)
function Ψ_HW(seb::SurfaceEnergyBalance{T,HøgstrømSHEBA}, ζ₁::Float64, ζ₂::Float64) where T
if ζ₁ <= 0 # neutral and unstable conditions (according to Høgstrøm, 1988)
# computed using WolframAlpha command "Integrate[ (1-(0.95*(1-11.6x)^(-1/2)))/x ] assuming x<0"
res = real( log(Complex(ζ₁)) + 1.9 * atanh(Complex((1 - 11.6 * ζ₁)^0.5)) -
......@@ -205,8 +262,10 @@ end
"""
Integrated stability function for momentum transport
Högström, 1988
SHEBA, Uttal et al., 2002, Grachev et al. 2007
"""
function Ψ_M(ζ₁::Real, ζ₂::Real)
function Ψ_M(seb::SurfaceEnergyBalance{T,HøgstrømSHEBA}, ζ₁::Float64, ζ₂::Float64) where T
if ζ₁ <= 0 # neutral and unstable conditions (according to Høgstrøm, 1988)
# computed using WolframAlpha command "Integrate[ (1-(1-19.3x)^(-1/4))/x ] assuming x<0"
# log(ζ₁) - 2*atan((1-19.3*ζ₁)^(1/4)) + 2*atanh((1-19.3*ζ₁)^(1/4)) -
......@@ -230,4 +289,37 @@ function Ψ_M(ζ₁::Real, ζ₂::Real)
end
end
export SurfaceEnergyBalance
"""
Integrated stability function for heat/water transport
Businger 1971
"""
function Ψ_HW(seb::SurfaceEnergyBalance{T,Businger},ζ₁::Float64, ζ₂::Float64) where T
let γₕ = seb.sebparams.γₕ,
βₕ = seb.sebparams.βₕ;
if ζ₁ <= 0 # neutral and unstable conditions (according to Businger, 1971)
return 2*log( ( (1-γₕ*ζ₁)^0.5 + 1 ) / ( (1-γₕ*ζ₂)^0.5 + 1 ) ); # eq. (15) in Byun 1990
else # stable stratification (according to Busigner, 1971)
return -βₕ*(ζ₁-ζ₂); # eq. (13) in Byun 1990
end
end
end
"""
Integrated stability function for momentum transport
Busigner 1971
"""
function Ψ_M(seb::SurfaceEnergyBalance{T,Businger},ζ₁::Float64, ζ₂::Float64) where T
let γₘ = seb.sebparams.γₘ,
βₘ = seb.sebparams.βₘ;
if ζ₁ <= 0 # neutral and unstable conditions (according to Businger, 1971)
x=(1-γₘ*ζ₁)^0.25;
x₀=(1-γₘ*ζ₂)^0.25;
return 2*log( (1+x)/(1+x₀) ) + log( (1+x^2) / (1+x₀^2) ) -
2*atan(x) + 2*atan(x) ; # eq. (14) in Byun 1990
else # stable stratification (according to Busigner, 1971)
return -βₘ*(ζ₁-ζ₂); # eq. (12) in Byun 1990
end
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment