diff --git a/docs/Project.toml b/docs/Project.toml
index dc84d5869..6e582add3 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -9,4 +9,4 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
-CloudMicrophysics = "0.5"
\ No newline at end of file
+CloudMicrophysics = "0.5"
diff --git a/docs/make.jl b/docs/make.jl
index 8cd3eed9a..c6c786eee 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -5,6 +5,8 @@ using EnsembleKalmanProcesses,
     Documenter,
     Plots,  # so that Literate.jl does not capture precompilation output
     Literate
+using Downloads
+using DelimitedFiles
 
 # Gotta set this environment variable when using the GR run-time on CI machines.
 # This happens as examples will use Plots.jl to make plots and movies.
@@ -19,6 +21,7 @@ examples_for_literation = [
     "LossMinimization/loss_minimization.jl",
     "SparseLossMinimization/loss_minimization_sparse_eki.jl",
     "AerosolActivation/aerosol_activation.jl",
+    "SurfaceFluxExample/kappa_calibration.jl",
 ]
 
 if isempty(get(ENV, "CI", ""))
@@ -42,6 +45,31 @@ for example in examples_for_literation
     )
 end
 
+# The purpose of the code below is to remove the occurrences of the Documenter @example block
+# created by Literate in order to prevent Documenter from evaluating the code blocks from 
+# kappa_calibration.jl. The reason the code blocks do not work is due to a package compatibility
+# mismatch, i.e. aerosol_activation is only compatible with CloudMicrophysics v0.5, and CloudMicrophysics v0.5
+# is only compatible with SurfaceFluxes v0.3. However, kappa_calibration.jl requires a new version of SurfaceFluxes, 
+# version 0.6, making it impossible to evaluate the code blocks in the markdown file kappa_calibration.md without
+# running into errors. Thus, we filter out the @example blocks and merely display the code on the docs.
+
+# Another reason we cannot evaluate the code blocks in kappa_calibration is that kappa_calibration depends
+# on locally downloaded data. Because we cannot download data to the remote repository, it is never plausible
+# to run kappa_calibration remotely.
+
+# read file and copy over modified
+kappa_md_file = open("docs/src/literated/kappa_calibration.md", "r")
+all_lines = string("")
+while (!eof(kappa_md_file))
+    line = readline(kappa_md_file) * "\n"
+    line = replace(line, "@example kappa_calibration" => "")
+    global all_lines *= line
+end
+
+# write to file
+kappa_md_file = open("docs/src/literated/kappa_calibration.md", "w")
+write(kappa_md_file, all_lines)
+close(kappa_md_file)
 #----------
 
 api = [
@@ -66,6 +94,7 @@ examples = [
     "Aerosol activation" => "literated/aerosol_activation.md",
     "TOML interface" => "examples/sinusoid_example_toml.md",
     "HPC interfacing example: ClimateMachine" => "examples/ClimateMachine_example.md",
+    "Surface Fluxes" => "literated/kappa_calibration.md",
     "Template" => "examples/template_example.md",
 ]
 
diff --git a/docs/src/assets/kappa_calibration_plot1.png b/docs/src/assets/kappa_calibration_plot1.png
new file mode 100644
index 000000000..4d0914193
Binary files /dev/null and b/docs/src/assets/kappa_calibration_plot1.png differ
diff --git a/docs/src/assets/kappa_calibration_plot2.png b/docs/src/assets/kappa_calibration_plot2.png
new file mode 100644
index 000000000..0fe31af33
Binary files /dev/null and b/docs/src/assets/kappa_calibration_plot2.png differ
diff --git a/examples/SurfaceFluxExample/Project.toml b/examples/SurfaceFluxExample/Project.toml
new file mode 100644
index 000000000..2c2319fe4
--- /dev/null
+++ b/examples/SurfaceFluxExample/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
diff --git a/examples/SurfaceFluxExample/alt_kappa_calibration.jl b/examples/SurfaceFluxExample/alt_kappa_calibration.jl
new file mode 100644
index 000000000..0dfac3dda
--- /dev/null
+++ b/examples/SurfaceFluxExample/alt_kappa_calibration.jl
@@ -0,0 +1,131 @@
+# # Alternative Kappa Calibration Example
+# ## Overview
+#=
+In this example, just like in kappa_calibration.jl, we use the inverse problem to calibrate the von-karman constant, κ in
+the equation: u(z) = u^* / κ log (z / z0), 
+which represents the wind profile in Monin-Obukhov
+Similarity Theory (MOST) formulations. We use the same dataset: https://turbulence.pha.jhu.edu/Channel_Flow.aspx
+
+Instead of using u^* as an observable, we use the dataset's u, and each ensemble member will estimate u
+through the profile equation u(z) = u^* / κ log (z / z0).
+=#
+
+# ## Prerequisites
+#=
+[EnsembleKalmanProcess.jl](https://github.com/CliMA/EnsembleKalmanProcesses.jl),
+=#
+
+# ## Example
+
+# First, we import relevant modules.
+using LinearAlgebra, Random
+using Distributions, Plots
+using EnsembleKalmanProcesses
+using EnsembleKalmanProcesses.ParameterDistributions
+const EKP = EnsembleKalmanProcesses
+
+using Downloads
+using DelimitedFiles
+
+FT = Float64
+
+mkpath(joinpath(@__DIR__, "data")) # create data folder if not exists
+web_datafile_path = "https://turbulence.oden.utexas.edu/channel2015/data/LM_Channel_5200_mean_prof.dat"
+localfile = "data/profiles.dat"
+Downloads.download(web_datafile_path, localfile)
+data_mean_velocity = readdlm("data/profiles.dat", skipstart = 112) ## We skip 72 lines (header) and 40(laminar layer)
+
+web_datafile_path = "https://turbulence.oden.utexas.edu/channel2015/data/LM_Channel_5200_mean_stdev.dat"
+localfile = "data/vel_stdev.dat"
+Downloads.download(web_datafile_path, localfile)
+# We skip 72 lines (header) and 40(laminar layer)
+data_stdev_velocity = readdlm("data/vel_stdev.dat", skipstart = 112)
+
+# We extract the required info for this problem
+u_star_obs = 4.14872e-02 # add noise later
+z0 = FT(0.0001)
+κ = 0.4
+
+# turn u into distributions
+u = data_mean_velocity[:, 3] * u_star_obs
+z = data_mean_velocity[:, 1]
+u = u[1:(length(u) - 1)] # filter out last element because σᵤ is only of length 727, not 728
+z = z[1:(length(z) - 1)]
+
+σᵤ = data_stdev_velocity[:, 4] * u_star_obs
+dist_u = Array{Normal{Float64}}(undef, length(u))
+for i in 1:length(u)
+    dist_u[i] = Normal(u[i], σᵤ[i])
+end
+
+# u(z) = u^* / κ log (z / z0)
+function physical_model(parameters, inputs)
+    κ = parameters[1] # this is being updated by the EKP iterator
+    (; u_star_obs, z, z0) = inputs
+    u_profile = u_star_obs ./ κ .* log.(z ./ z0)
+    return u_profile
+end
+
+function G(parameters, inputs, u_profile = nothing)
+    if (isnothing(u_profile))
+        u_profile = physical_model(parameters, inputs)
+    end
+    return [maximum(u_profile) - minimum(u_profile), mean(u_profile)]
+end
+
+Γ = 0.0001 * I
+η_dist = MvNormal(zeros(2), Γ)
+noisy_u_profile = [rand(dist_u[i]) for i in 1:length(u)]
+y = G(nothing, nothing, noisy_u_profile)
+
+parameters = (; κ)
+inputs = (; u_star_obs, z, z0)
+# y = G(parameters, inputs) .+ rand(η_dist)
+
+# Assume that users have prior knowledge of approximate truth.
+# (e.g. via physical models / subset of obs / physical laws.)
+prior_u1 = constrained_gaussian("κ", 0.35, 0.25, 0, Inf);
+prior = combine_distributions([prior_u1])
+
+# Set up the initial ensembles
+N_ensemble = 5;
+N_iterations = 10;
+
+rng_seed = 41
+rng = Random.MersenneTwister(rng_seed)
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble);
+
+# Define EKP and run iterative solver for defined number of iterations
+ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, y, Γ, Inversion(); rng = rng)
+
+for n in 1:N_iterations
+    params_i = get_ϕ_final(prior, ensemble_kalman_process)
+    G_ens = hcat([G(params_i[:, m], inputs) for m in 1:N_ensemble]...)
+    EKP.update_ensemble!(ensemble_kalman_process, G_ens)
+end
+
+# Mean values in final ensemble for the two parameters of interest reflect the "truth" within some degree of 
+# uncertainty that we can quantify from the elements of `final_ensemble`.
+final_ensemble = get_ϕ_final(prior, ensemble_kalman_process)
+mean(final_ensemble[1, :]) # [param, ens_no]
+
+
+ENV["GKSwstype"] = "nul"
+zrange = z
+plot(zrange, noisy_u_profile, c = :black, label = "Truth", linewidth = 2, legend = :bottomright)
+plot!(zrange, physical_model(parameters, inputs), c = :green, label = "Model truth", linewidth = 2)# reshape to convert from vector to matrix)
+plot!(
+    zrange,
+    [physical_model(get_ϕ(prior, ensemble_kalman_process, 1)[:, i]..., inputs) for i in 1:N_ensemble],
+    c = :red,
+    label = reshape(vcat(["Initial ensemble"], ["" for i in 1:(N_ensemble - 1)]), 1, N_ensemble), # reshape to convert from vector to matrix
+)
+plot!(
+    zrange,
+    [physical_model(final_ensemble[:, i]..., inputs) for i in 1:N_ensemble],
+    c = :blue,
+    label = reshape(vcat(["Final ensemble"], ["" for i in 1:(N_ensemble - 1)]), 1, N_ensemble),
+)
+xlabel!("Z")
+ylabel!("U")
+png("profile_plot")
diff --git a/examples/SurfaceFluxExample/kappa_calibration.jl b/examples/SurfaceFluxExample/kappa_calibration.jl
new file mode 100644
index 000000000..a00f5ece0
--- /dev/null
+++ b/examples/SurfaceFluxExample/kappa_calibration.jl
@@ -0,0 +1,265 @@
+# # Kappa Calibration Example
+# ## Overview
+#=
+In this example, we use the inverse problem to calibrate the von-karman constant, κ in
+the equation: u(z) = u^* / κ log (z / z0), 
+which represents the wind profile in Monin-Obukhov
+Similarity Theory (MOST) formulations. In order to recover the empirically determined κ = 0.4, 
+we use data from the John Hopkins Tubulence Channel Flow, which offers DNS simulations of a channel flow with 
+smooth wall boundary conditions, i.e. z0m ≈ 0 m. The dataset can be found here: https://turbulence.pha.jhu.edu/Channel\_Flow.aspx.
+We use the dataset's u^* as an observable, and each ensemble member estimates u^* through the
+SurfaceFluxes.jl function surface\_conditions, see https://github.com/CliMA/SurfaceFluxes.jl
+In order to calculate u^*, the function surface\_conditions is provided a set of thermodynamic params,
+a functional form for stability functions (Businger, Gryanick, Grachev), and the constants corresponding
+to that functional form. In this example, we elect the Businger functions. 
+=#
+
+# ## Prerequisites
+#=
+This example depends on standard Julia packages as well as CliMA packages:
+[EnsembleKalmanProcess.jl](https://github.com/CliMA/EnsembleKalmanProcesses.jl),
+[CLIMAParameters.jl](https://github.com/CliMA/CLIMAParameters.jl),
+[SurfaceFluxes v0.6](https://github.com/CliMA/SurfaceFluxes.jl),
+[Thermodynamics v0.10](https://github.com/CliMA/Thermodynamics.jl)
+Note that this example is only compatible with certain versions of these packages.
+=#
+
+# ## Example
+
+# First, we import relevant modules.
+using LinearAlgebra, Random
+using Distributions, Plots
+using EnsembleKalmanProcesses
+using EnsembleKalmanProcesses.ParameterDistributions
+const EKP = EnsembleKalmanProcesses
+
+using Downloads
+using DelimitedFiles
+
+using CLIMAParameters
+const CP = CLIMAParameters
+FT = Float64
+
+import SurfaceFluxes as SF
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TP
+import SurfaceFluxes.UniversalFunctions as UF
+import SurfaceFluxes.Parameters as SFP
+using StaticArrays: SVector
+include("setup_parameter_set.jl")
+
+#=
+Next, we download and read data from the John Hopkins Tubulence Channel Flow dataset,
+specifically those concerning mean velocity and its variance over various heights.
+The parameters defining the dataset are given by:
+- u\_star = 4.14872e-02
+- δ = 1.000
+- ν = 8.00000e-06
+- Re\_tau = 5185.897
+=#
+mkpath(joinpath(@__DIR__, "data")) # create data folder if not exists
+web_datafile_path = "https://turbulence.oden.utexas.edu/channel2015/data/LM_Channel_5200_mean_prof.dat"
+localfile = "data/profiles.dat"
+Downloads.download(web_datafile_path, localfile)
+data_mean_velocity = readdlm("data/profiles.dat", skipstart = 112) ## We skip 72 lines (header) and 40(laminar layer)
+
+# We extract the required info for this problem
+u_star_obs = 4.14872e-02
+z = data_mean_velocity[:, 1]
+u = data_mean_velocity[:, 3] * u_star_obs
+
+# Next, we define our physical model, where we first define thermodynamic parameters
+# and MOST parameters to pass into the surface\_conditions function from SurfaceFluxes.jl.
+# We define the MOST stability functions to be of the Businger type.
+"""
+    physical_model(inputs, parameters)
+
+Takes in kappa and mean states, producing a ustar (or u(z) profile) ensemble for each horizontal point.
+
+Inputs:
+    inputs{Array}(len)          in this case containing u and z
+    parameters{NamedTuple}      in this case κ
+
+Example:
+
+`inputs` are the measured profiles, which would be inputs to `SurfaceFluxes.jl` -> See JH DNS Database.
+EKP requires `parameters` in `Vector` form - `Tuples` are not permitted.
+"""
+function physical_model(parameters, inputs)
+    κ = parameters[1] # this is being updated by the EKP iterator
+    (; u, z) = inputs
+
+    # First, we set up thermodynamic parameters
+    ## This line initializes a toml dict, where we will extract parameters from
+    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+    param_set = create_parameters(toml_dict, UF.BusingerType())
+    thermo_params = SFP.thermodynamics_params(param_set)
+
+    ## in this idealized case, we assume dry isothermal conditions
+    ts_sfc = TD.PhaseEquil_ρθq(thermo_params, FT(1), FT(300), FT(0))
+    ts_in = TD.PhaseEquil_ρθq(thermo_params, FT(1), FT(300), FT(0))
+
+    # Next, we set up SF parameters
+    ## An alias for each constant we need
+    aliases = ["Pr_0_Businger", "a_m_Businger", "a_h_Businger", "ζ_a_Businger", "γ_Businger"]
+    sf_pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions")
+    sf_pairs = (; sf_pairs...) # convert parameter pairs to NamedTuple
+    ## change the keys from their alias to more concise keys
+    sf_pairs = (;
+        Pr_0 = sf_pairs.Pr_0_Businger,
+        a_m = sf_pairs.a_m_Businger,
+        a_h = sf_pairs.a_h_Businger,
+        ζ_a = sf_pairs.ζ_a_Businger,
+        γ = sf_pairs.γ_Businger,
+    )
+    ufp = UF.BusingerParams{FT}(; sf_pairs...) # initialize Businger params
+
+    κ_nt = (; von_karman_const = κ)
+
+    # Now, we initialize the variable surf\_flux\_params, which we will eventually pass into 
+    # surface_conditions along with mean wind data
+    UFP = typeof(ufp)
+    TPtype = typeof(thermo_params)
+    surf_flux_params = SF.Parameters.SurfaceFluxesParameters{FT, UFP, TPtype}(; κ_nt..., ufp, thermo_params)
+
+    # Now, we loop over all the observations and call SF.surface_conditions to estimate u^*
+    u_star = zeros(length(u))
+    for i in 1:lastindex(u)
+        u_in = u[i]
+        v_in = FT(0)
+        z_in = z[i]
+        u_in = SVector{2, FT}(u_in, v_in)
+        u_sfc = SVector{2, FT}(FT(0), FT(0))
+
+        state_sfc = SF.SurfaceValues(FT(0), u_sfc, ts_sfc)
+        state_in = SF.InteriorValues(z_in, u_in, ts_in)
+
+        # We provide a few additional parameters for SF.surface_conditions
+        z0m = z0b = FT(0.0001)
+        gustiness = FT(0)
+        kwargs = (; state_in, state_sfc, z0m, z0b, gustiness)
+        sc = SF.ValuesOnly{FT}(; kwargs...)
+
+        # Now, we call surface_conditions and store the calculated ustar:
+        sf = SF.surface_conditions(surf_flux_params, sc)
+        u_star[i] = sf.ustar  # TODO: also try for u_profiles[i, :] = sf.u_profile(z)
+    end
+
+    return u_star
+end
+
+# Here, we define G, which returns observable values given the parameters and inputs
+# from the dataset. The observable we elect is the mean of the calculated ustar across
+# all z, which is eventually compared to the actual observed ustar.
+function G(parameters, inputs)
+    u_star = physical_model(parameters, inputs)
+    u_star_mean = mean(u_star) # H map
+    return [u_star_mean]
+end
+
+# Here, we define the true value of the parameters we wish to recover
+κ = 0.4
+theta_true = [κ]
+u_star_obs = 4.14872e-02
+
+# Define the arguments to be passed into G:
+parameters = (; κ)
+u = u[1:(end - 1)] # we remove the last line because we want different surface state conditions
+z = z[1:(end - 1)]
+inputs = (; u, z)
+
+# Next, we define y, which is the noisy observation. Because we already have the truth observation,
+# we add noise to the observed u^* and store it in y. Refer to Cleary et al 2021 for more details.
+# We choose a noise scaling constant of 0.0005, which is small because the order of magnitude of 
+# u_* is 10^-2, and we don't want the noise to be greater than our observation.
+Γ = 0.0005 * I
+η_dist = MvNormal(zeros(1), Γ)
+y = [u_star_obs] .+ rand(η_dist) # (H ⊙ Ψ ⊙ T^{-1})(θ) + η from Cleary et al 2021
+## we can try these definitions of y later: y = G(inputs, parameters) .+ rand(η_dist)
+## y = u - u_star/κ (log(z/z0))
+
+# Assume that users have prior knowledge of approximate truth.
+# (e.g. via physical models / subset of obs / physical laws.)
+prior_u1 = constrained_gaussian("κ", 0.35, 0.25, 0, Inf);
+prior = combine_distributions([prior_u1])
+
+# Set up the initial ensembles
+N_ensemble = 10;
+N_iterations = 50;
+
+rng_seed = 41
+rng = Random.MersenneTwister(rng_seed)
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble);
+
+# Define EKP and run iterative solver for defined number of iterations
+ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, y, Γ, Inversion(); rng = rng)
+
+for n in 1:N_iterations
+    ## get_ϕ_final returns the most recently updated constrained parameters, which it used to make the
+    ## next model forward and thus the next update 
+    params_i = get_ϕ_final(prior, ensemble_kalman_process)
+    ## calculate the forwarded model values
+    G_ens = hcat([G(params_i[:, m], inputs) for m in 1:N_ensemble]...)
+    EKP.update_ensemble!(ensemble_kalman_process, G_ens)
+end
+
+# Mean values in final ensemble for the two parameters of interest reflect the "truth" within some degree of 
+# uncertainty that we can quantify from the elements of `final_ensemble`.
+final_ensemble = get_ϕ_final(prior, ensemble_kalman_process)
+mean(final_ensemble[1, :]) # [param, ens_no]
+
+# To visualize the success of the inversion, we plot model with 3 different forms of the truth: 
+# - The absolute truth of u^* given by the dataset
+# - y, the noisy observation we used to calibrate our model parameter κ
+# - The output of the physical model given the true κ = 0.4
+
+# We then compare them to the initial ensemble and the final ensemble.
+zrange = z
+initial_ensemble = get_ϕ(prior, ensemble_kalman_process, 1)
+ENV["GKSwstype"] = "nul"
+plot(
+    zrange,
+    physical_model(theta_true..., inputs),
+    c = :black,
+    label = "Model Truth",
+    legend = :bottomright,
+    linewidth = 2,
+    linestyle = :dash,
+)
+plot!(
+    zrange,
+    ones(length(zrange)) .* y,
+    c = :black,
+    label = "y",
+    legend = :bottomright,
+    linewidth = 2,
+    linestyle = :dot,
+)
+plot!(zrange, ones(length(zrange)) .* u_star_obs, c = :black, label = "Truth u*", legend = :bottomright, linewidth = 2)
+plot!(
+    zrange,
+    [physical_model(get_ϕ(prior, ensemble_kalman_process, 1)[:, i]..., inputs) for i in 1:N_ensemble],
+    c = :red,
+    label = reshape(vcat(["Initial ensemble"], ["" for i in 1:(N_ensemble - 1)]), 1, N_ensemble), # reshape to convert from vector to matrix
+)
+plot!(
+    zrange,
+    [physical_model(final_ensemble[:, i]..., inputs) for i in 1:N_ensemble],
+    c = :blue,
+    label = reshape(vcat(["Final ensemble"], ["" for i in 1:(N_ensemble - 1)]), 1, N_ensemble),
+)
+xlabel!("Z")
+ylabel!("U^*")
+png("our_plot")
+# ![see plot: ](../assets/kappa_calibration_plot1.png)
+
+# We also plot the constrained κ values across all ensembles before and after the EKI process in a histogram.
+histogram(initial_ensemble[1, :], label = "initial")
+histogram!(final_ensemble[1, :], label = "final")
+xlabel!("κ")
+ylabel!("# of Ensembles")
+png("final_and_initial_ensemble")
+# ![see plot: ](../assets/kappa_calibration_plot2.png)
+
+# Evidently, EKI was highly successful at recovering the von karman constant κ = 0.4. This process will be extended to recover
+# stability function parameters such as a\_m, a\_h, b\_m, b\_h, and Pr\_0. 
diff --git a/examples/SurfaceFluxExample/setup_parameter_set.jl b/examples/SurfaceFluxExample/setup_parameter_set.jl
new file mode 100644
index 000000000..d377e2a80
--- /dev/null
+++ b/examples/SurfaceFluxExample/setup_parameter_set.jl
@@ -0,0 +1,130 @@
+# sets up parameters
+# const sf_params = SF.Parameters.SurfaceFluxesParameters{
+#     FT,
+#     SF.UniversalFunctions.BusingerParams{FT},
+#     TP.ThermodynamicsParameters{FT},
+# }(
+#     0.4f0,
+#     SF.UniversalFunctions.BusingerParams{FT}(0.74f0, 4.7f0, 4.7f0, 2.5f0, 4.45f0),
+#     TP.ThermodynamicsParameters{FT}(
+#         273.16f0,
+#         100000.0f0,
+#         1859.0f0,
+#         4181.0f0,
+#         2100.0f0,
+#         2.5008f6,
+#         2.8344f6,
+#         611.657f0,
+#         273.16f0,
+#         273.15f0,
+#         150.0f0,
+#         1000.0f0,
+#         298.15f0,
+#         6864.8f0,
+#         10513.6f0,
+#         0.2857143f0,
+#         8.31446f0,
+#         0.02897f0,
+#         0.01801528f0,
+#         290.0f0,
+#         220.0f0,
+#         9.80616f0,
+#         233.0f0,
+#         1.0f0,
+#     ),
+# )
+
+
+import CLIMAParameters as CP
+import SurfaceFluxes as SF
+import SurfaceFluxes.UniversalFunctions as UF
+import Thermodynamics as TD
+
+function create_uf_parameters(toml_dict, ::UF.GryanikType)
+    FT = CP.float_type(toml_dict)
+
+    aliases = ["Pr_0_Gryanik", "a_m_Gryanik", "a_h_Gryanik", "b_m_Gryanik", "b_h_Gryanik", "ζ_a_Gryanik", "γ_Gryanik"]
+
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions")
+    pairs = (; pairs...) # convert to NamedTuple
+
+    pairs = (;
+        Pr_0 = pairs.Pr_0_Gryanik,
+        a_m = pairs.a_m_Gryanik,
+        a_h = pairs.a_h_Gryanik,
+        b_m = pairs.b_m_Gryanik,
+        b_h = pairs.b_h_Gryanik,
+        ζ_a = pairs.ζ_a_Gryanik,
+        γ = pairs.γ_Gryanik,
+    )
+    return UF.GryanikParams{FT}(; pairs...)
+end
+
+function create_uf_parameters(toml_dict, ::UF.BusingerType)
+    FT = CP.float_type(toml_dict)
+    aliases = ["Pr_0_Businger", "a_m_Businger", "a_h_Businger", "ζ_a_Businger", "γ_Businger"]
+
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions")
+    pairs = (; pairs...) # convert to NamedTuple
+
+    pairs = (;
+        Pr_0 = pairs.Pr_0_Businger,
+        a_m = pairs.a_m_Businger,
+        a_h = pairs.a_h_Businger,
+        ζ_a = pairs.ζ_a_Businger,
+        γ = pairs.γ_Businger,
+    )
+    return UF.BusingerParams{FT}(; pairs...)
+end
+
+function create_uf_parameters(toml_dict, ::UF.GrachevType)
+    FT = CP.float_type(toml_dict)
+    aliases = [
+        "Pr_0_Grachev",
+        "a_m_Grachev",
+        "a_h_Grachev",
+        "b_m_Grachev",
+        "b_h_Grachev",
+        "c_h_Grachev",
+        "ζ_a_Grachev",
+        "γ_Grachev",
+    ]
+
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions")
+    pairs = (; pairs...) # convert to NamedTuple
+
+    pairs = (;
+        Pr_0 = pairs.Pr_0_Grachev,
+        a_m = pairs.a_m_Grachev,
+        a_h = pairs.a_h_Grachev,
+        b_m = pairs.b_m_Grachev,
+        b_h = pairs.b_h_Grachev,
+        c_h = pairs.c_h_Grachev,
+        ζ_a = pairs.ζ_a_Grachev,
+        γ = pairs.γ_Grachev,
+    )
+    return UF.GrachevParams{FT}(; pairs...)
+end
+
+function create_parameters(toml_dict, ufpt)
+    FT = CP.float_type(toml_dict)
+
+    ufp = create_uf_parameters(toml_dict, ufpt)
+    AUFP = typeof(ufp)
+
+    aliases = string.(fieldnames(TD.Parameters.ThermodynamicsParameters))
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+    thermo_params = TD.Parameters.ThermodynamicsParameters{FT}(; pairs...)
+    TP = typeof(thermo_params)
+
+    aliases = ["von_karman_const"]
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "SurfaceFluxesParameters")
+    override = return SFP.SurfaceFluxesParameters{FT, AUFP, TP}(; pairs..., ufp, thermo_params)
+end
+
+function override_climaatmos_defaults(defaults::NamedTuple, overrides::NamedTuple)
+    intersect_keys = intersect(keys(defaults), keys(overrides))
+    intersect_vals = getproperty.(Ref(overrides), intersect_keys)
+    intersect_overrides = (; zip(intersect_keys, intersect_vals)...)
+    return merge(defaults, intersect_overrides)
+end