Physick & Garratt Roughness Sublayer Model

.dev/Project.toml

@@ -2,4 +2,4 @@
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 
 [compat]
-JuliaFormatter = "0.22"
+JuliaFormatter = "1"

Project.toml

@@ -4,10 +4,12 @@ authors = ["Climate Modeling Alliance"]
 version = "0.12.0"
 
 [deps]
+CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
-ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 
 [weakdeps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
@@ -18,6 +20,7 @@ CreateParametersExt = "ClimaParams"
 [compat]
 ClimaParams = "0.10"
 DocStringExtensions = "0.8, 0.9"
+QuadGK = "2.11.2"
 RootSolvers = "0.4"
 Thermodynamics = "0.12.8"
 julia = "1.9"

docs/bibliography.bib

@@ -2,6 +2,29 @@
 #    Last author name (titlecase), followed by
 #    (no characters in-between) the year.
 
+
+@article{Bonan2019,
+  title = {Climate Change and Terrestrial Ecosystem Modeling},
+  author = {Bonan, G},
+  journal = {Cambridge University Press},
+  doi = {10.1017/9781107339217},
+  year = {2019},
+  publisher = {Cambridge University Press}
+}
+
+
+@article{Physick1995,
+  title = {Incorporation of a high-roughness lower boundary into a mesoscale model for studies of dry deposition over complex terrain},
+  author = {Physick, W  and Garratt, R},
+  journal = {Boundary-Layer Meteorology},
+  volume = {74},
+  number = {1},
+  doi = {10.1007/BF00715710},
+  pages = {55-71},
+  year = {1995},
+  publisher = {Springer}
+}
+
 @article{Cheng2005,
   title = {Flux-profile Relationships for Wind Speed and Temperature in the Stable Atmospheric Boundary Layer},
   author = {Chenge, Y and Brutsaert, W},

docs/make.jl

@@ -11,7 +11,8 @@ pages = Any[
     "Home" => "index.md",
     "References" => "References.md",
     "Equations" => "SurfaceFluxes.md",
-    "Universal Functions" => "UniversalFunctions.md"
+    "Universal Functions" => "UniversalFunctions.md",
+    "Sample Profiles" => "SampleProfiles.md"
 ]
 
 mathengine = MathJax(Dict(

docs/src/SampleProfiles.md

@@ -0,0 +1,29 @@
+# Sample Vertical Profiles of Wind Speed and Temperature
+
+SurfaceFluxes.jl provides profile recovery functions for the roughness sublayer using the Physick and Garratt (1995) formulation. Here, we use these functions to reproduce vertical profiles of wind speed and temperature for different conditions of atmospheric stability, reproducing Bonan 2019 Figure 6.4 with and without canopy correction.
+
+```@example
+include("plot_profiles.jl")
+```
+
+# Fig 6.4 (a) 
+
+![](Fig6.4a_profile.svg)
+![](Fig6.4a_canopy_profile.svg)
+
+# Fig 6.4 (b) 
+
+![](Fig6.4b_profile.svg)
+![](Fig6.4b_canopy_profile.svg)
+
+# Fig 6.4 (c)
+
+![](Fig6.4c_profile.svg)
+![](Fig6.4c_canopy_profile.svg)
+
+# Fig 6.4 (d) 
+
+![](Fig6.4d_profile.svg)
+![](Fig6.4d_canopy_profile.svg)
+
+

docs/src/plot_profiles.jl

@@ -0,0 +1,292 @@
+using Plots
+
+using SurfaceFluxes
+const SF = SurfaceFluxes
+SurfaceFluxes.error_on_non_convergence() = true
+
+import SurfaceFluxes.UniversalFunctions as UF
+import Thermodynamics
+Thermodynamics.print_warning() = false
+
+include(joinpath(pkgdir(SurfaceFluxes), "parameters", "create_parameters.jl"))
+const SFP = SurfaceFluxes.Parameters
+const FT = Float64;
+
+### Generate parameter lists
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+param_set = create_parameters(toml_dict, UF.Gryanik())
+thermo_params = SFP.thermodynamics_params(param_set)
+uft = SFP.universal_func_type(param_set)
+### 
+
+""" 
+Here, we reproduce Figure 6.4 from Bonan (2019) Chapter 6.
+
+# References
+- [Bonan2019](@cite) (Chapter 6 Figure 6.4)
+
+# Original Research
+- [Physick2019](@cite)
+
+"""
+T₀ = FT(273.15)
+z_star = FT(49)
+d = FT(19)
+h_c = FT(22)
+z0m = FT(0.6) # Figure 6.4
+z0b = FT(0.135) * z0m
+
+u_star_stable = FT(0.13)
+u_star_unstable = FT(0.4)
+θ_star_stable = FT(0.06)
+θ_star_unstable = FT(-0.5)
+
+κ = SFP.von_karman_const(param_set)
+
+# For these measured values, see Chapter 6.6 page 89 of Bonan (2019)
+z_measured = FT[21, 29, 21];
+u_measured = FT[1.0, 2.1, 1.0];
+θ_measured = FT[29.0, 28.1, 29.5] .+ T₀
+
+p_sfc = FT(99340) # Pa
+q_sfc = FT(0.0107) # kg/kg
+
+θ_sfc = FT(289.7) # K
+ts_sfc_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 99340.0, 289.7, 0.0107)
+ts_int_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 95342.0, 298.0, 0.0085)
+state_in = SurfaceFluxes.InteriorValues(FT(100), (FT(1.0), FT(0)), ts_int_test)
+state_sfc = SurfaceFluxes.SurfaceValues(FT(0), (FT(0), FT(0)), ts_sfc_test)
+sc = SurfaceFluxes.ValuesOnly{FT}(; state_in, state_sfc, z0m, z0b)
+ts_sfc_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 100000.0, 289.7, 0.0)
+ts_int_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 99990.0, 298.0, 0.0)
+Z = collect(range(FT(d + 10^-10), stop = FT(100), length = 250))
+
+"""
+save_profile(param_set, sc, ca, L_MOs, Z, X_sfc, transport, uft, scheme, x_star, d)
+
+Saves profiles of variable X given values of Z coordinates. Follows Nishizawa equation (21,22)
+
+## Arguments
+  - param_set: Abstract Parameter Set containing physical, thermodynamic parameters.
+  - sc: Container for surface conditions based on known combination
+        of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment
+  - L_MOs: Monin-Obukhov length(s)
+  - Z: Z coordinate(s) (within surface layer) for which variable values are required
+  - X_sfc: For variable X, values at interior and surface nodes
+  - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients)
+  - uft: A Universal Function type, (returned by, e.g., Businger())
+  - scheme: Discretization scheme (currently supports FD and FV)
+  - rsl : Roughness Sublayer Formulation (e.g. NoRSL, PhysickRSL, DeRidderRSL)
+  - x_star: characteristic scale for variable x
+  - d: Displacement height (measure of the spatial lengthscale of the effect of the canopy roughness on near-wall turbulence)
+"""
+function save_profile(
+    param_set::SurfaceFluxes.APS,
+    sc::SurfaceFluxes.AbstractSurfaceConditions,
+    L_MOs::Array{FT, 1},
+    Z::Array{FT, 1},
+    X_sfc,
+    transport,
+    uft::UF.AUFT,
+    scheme::Union{SurfaceFluxes.FVScheme, SurfaceFluxes.FDScheme},
+    rsl::SurfaceFluxes.AbstractRoughnessSublayerType,
+    x_star,
+    d;
+    title = nothing,
+    xlims = nothing,
+    ylims = nothing,
+    xlabel = "u(z)",
+    ylabel = "z",
+    fig_prefix = "",
+    xaxis = :identity,
+    yaxis = :identity,
+)
+    Plots.plot()
+    for L_MO in L_MOs
+        x_i = map(Z) do z
+            Zi = typeof(rsl) == SurfaceFluxes.NoRSL ? FT(z - d) : FT(z)
+            state_in = SurfaceFluxes.InteriorValues(FT(Zi), (FT(1.0), FT(0)), ts_int_test)
+            state_sfc = SurfaceFluxes.SurfaceValues(FT(0), (FT(0), FT(0)), ts_sfc_test)
+            sc = SurfaceFluxes.ValuesOnly{FT}(; state_in, state_sfc, z0m, z0b)
+            rsc = SurfaceFluxes.surface_conditions(param_set, sc, SurfaceFluxes.FDScheme())
+            dx = SurfaceFluxes.recover_profile(param_set, sc, L_MO, Zi, X_sfc, x_star, transport, uft, scheme, rsl)
+        end
+
+        uf = UF.universal_func(uft, L_MO, SFP.uf_params(param_set))
+        _π_group = FT(UF.π_group(uf, transport))
+
+        Δx = @. (x_i - X_sfc)
+
+        Plots.plot!(Δx, Z, label = "L_MO = $L_MO")
+        Plots.plot!(; title, xlabel, ylabel, ylims, xlims, grid = :off, legend = :outerright, titlefontalign = :center)
+
+        Plots.savefig("$(fig_prefix)_profile.png")
+    end
+end
+
+
+
+# Save profile for unstable and stable L_MO set
+Stable_L_MOs = FT[30, 50, 1000]
+Unstable_L_MOs = FT[-10, -50, -1000]
+
+testcanopy = SurfaceFluxes.SparseCanopy{FT}(d, z_star)
+PhysickRSL = SurfaceFluxes.PhysickRSL(testcanopy)
+NoRSL = SurfaceFluxes.NoRSL()
+
+
+# Bonan2019 Fig. 6.4a (With PG95 RSL)
+save_profile(
+    param_set,
+    sc,
+    Unstable_L_MOs,
+    Z,
+    FT(0),
+    UF.MomentumTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    PhysickRSL,
+    u_star_unstable,
+    d;
+    xlims = (0, 4),
+    ylims = (15, 50),
+    fig_prefix = "Fig6.4a_canopy",
+    title = "Vertical Profile of Wind Velocity (PG95 RSL)",
+)
+
+# Bonan2019 Fig. 6.4a (No RSL Model)
+save_profile(
+    param_set,
+    sc,
+    Unstable_L_MOs,
+    Z,
+    FT(0),
+    UF.MomentumTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    NoRSL,
+    u_star_unstable,
+    d;
+    xlims = (0, 4),
+    ylims = (15, 50),
+    fig_prefix = "Fig6.4a",
+    title = "Vertical Profile of Wind Velocity (No RSL Model)",
+)
+
+# Bonan2019 Fig. 6.4b (With PG95 RSL)
+save_profile(
+    param_set,
+    sc,
+    Unstable_L_MOs,
+    Z,
+    θ_sfc,
+    UF.HeatTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    PhysickRSL,
+    θ_star_unstable,
+    d;
+    xlims = (-8, 0),
+    ylims = (15, 50),
+    xlabel = "θ - θ_sfc",
+    fig_prefix = "Fig6.4b_canopy",
+    title = "Vertical Profile of Temperature (PG95 RSL)",
+)
+
+# Bonan2019 Fig. 6.4b (No RSL Model)
+save_profile(
+    param_set,
+    sc,
+    Unstable_L_MOs,
+    Z,
+    θ_sfc,
+    UF.HeatTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    NoRSL,
+    θ_star_unstable,
+    d;
+    xlims = (-8, 0),
+    ylims = (15, 50),
+    xlabel = "θ- θ_sfc",
+    fig_prefix = "Fig6.4b",
+    title = "Vertical Profile of Temperature (No RSL Model)",
+)
+
+# Bonan2019 Fig. 6.4c (With PG95 RSL)
+save_profile(
+    param_set,
+    sc,
+    Stable_L_MOs,
+    Z,
+    FT(0),
+    UF.MomentumTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    PhysickRSL,
+    u_star_stable,
+    d;
+    xlims = (0, 4),
+    ylims = (15, 50),
+    fig_prefix = "Fig6.4c_canopy",
+    title = "Vertical Profile of Wind Velocity (PG95 RSL)",
+)
+
+# Bonan2019 Fig. 6.4c (No RSL Model)
+save_profile(
+    param_set,
+    sc,
+    Stable_L_MOs,
+    Z,
+    FT(0),
+    UF.MomentumTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    NoRSL,
+    u_star_stable,
+    d;
+    xlims = (0, 4),
+    ylims = (15, 50),
+    fig_prefix = "Fig6.4c",
+    title = "Vertical Profile of Wind Velocity (No RSL Model)",
+)
+
+# Bonan2019 Fig. 6.4d (With PG95 RSL)
+save_profile(
+    param_set,
+    sc,
+    Stable_L_MOs,
+    Z,
+    θ_sfc,
+    UF.HeatTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    PhysickRSL,
+    θ_star_stable,
+    d;
+    xlims = (0, 2),
+    ylims = (15, 50),
+    xlabel = "θ-θ_sfc",
+    fig_prefix = "Fig6.4d_canopy",
+    title = "Vertical Profile of Temperature (PG95 RSL)",
+)
+
+# Bonan2019 Fig. 6.4d (No RSL Model)
+save_profile(
+    param_set,
+    sc,
+    Stable_L_MOs,
+    Z,
+    θ_sfc,
+    UF.HeatTransport(),
+    uft,
+    SurfaceFluxes.FDScheme(),
+    NoRSL,
+    θ_star_stable,
+    d;
+    xlims = (0, 2),
+    ylims = (15, 50),
+    xlabel = "θ-θ_sfc",
+    fig_prefix = "Fig6.4d",
+    title = "Vertical Profile of Temperature (No RSL Model)",
+)