diff --git a/Project.toml b/Project.toml
index 97ab802e38..071b494824 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,7 +20,6 @@ Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"
@@ -61,7 +60,6 @@ Colors = "0.12"
 Dates = "1"
 Dierckx = "0.5"
 DiffEqBase = "6"
-Distributions = "0.25"
 DocStringExtensions = "0.8, 0.9"
 FastGaussQuadrature = "0.4, 0.5, 1"
 Insolation = "0.9.2"
diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl
index 6de6fe1aff..17cfb2fd0b 100644
--- a/regression_tests/ref_counter.jl
+++ b/regression_tests/ref_counter.jl
@@ -1,4 +1,7 @@
-163
+164
+
+# 164:
+# - Changed approximation for erf in calculation of some boundary conditions for EDMF.
 
 # 163:
 # - Fixed bug introduced in 162
diff --git a/src/prognostic_equations/edmfx_boundary_condition.jl b/src/prognostic_equations/edmfx_boundary_condition.jl
index 8f177c1bbd..20810c9df9 100644
--- a/src/prognostic_equations/edmfx_boundary_condition.jl
+++ b/src/prognostic_equations/edmfx_boundary_condition.jl
@@ -2,8 +2,6 @@
 ##### EDMFX SGS boundary condition
 #####
 
-import Distributions
-
 function sgs_scalar_first_interior_bc(
     ᶜz_int::FT,
     ᶜρ_int::FT,
@@ -49,16 +47,27 @@ function get_first_interior_variance(
     end
 end
 
+function approximate_inverf(x::FT) where {FT <: Real}
+    # From Wikipedia
+    a = FT(0.147)
+    term1 = (2 / (π * a) + log(1 - x^2) / 2)
+    term2 = log(1 - x^2) / a
+    term3 = sqrt(term1^2 - term2)
+
+    return sign(x) * sqrt(term3 - term1)
+end
+
+function guass_quantile(p::FT) where {FT <: Real}
+    return sqrt(2) * approximate_inverf(2p - 1)
+end
+
 function percentile_bounds_mean_norm(
     low_percentile::FT,
     high_percentile::FT,
 ) where {FT <: Real}
     gauss_int(x) = -exp(-x * x / 2) / sqrt(2 * pi)
-    xp_low = Distributions.quantile(Distributions.Normal(), low_percentile)
-    xp_high = Distributions.quantile(Distributions.Normal(), high_percentile)
-    return (gauss_int(xp_high) - gauss_int(xp_low)) / max(
-        Distributions.cdf(Distributions.Normal(), xp_high) -
-        Distributions.cdf(Distributions.Normal(), xp_low),
-        eps(FT),
-    )
+    xp_high = guass_quantile(high_percentile)
+    xp_low = guass_quantile(low_percentile)
+    return (gauss_int(xp_high) - gauss_int(xp_low)) /
+           max(high_percentile - low_percentile, eps(FT))
 end
diff --git a/test/Project.toml b/test/Project.toml
index 8273983868..3a01493f43 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,7 +4,6 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a"
 AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
 ClimaAtmos = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
@@ -15,6 +14,7 @@ ClimaCorePlots = "cf7c7e5a-b407-4c48-9047-11a94a308626"
 ClimaCoreSpectra = "c2caaa1d-32ae-4754-ba0d-80e7561362e9"
 ClimaCoreTempestRemap = "d934ef94-cdd4-4710-83d6-720549644b70"
 ClimaCoreVTK = "c8b6d40d-e815-466f-95ae-c48aefa668fa"
+ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
 CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -23,6 +23,7 @@ Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/test/prognostic_equations.jl b/test/prognostic_equations.jl
new file mode 100644
index 0000000000..8d416af793
--- /dev/null
+++ b/test/prognostic_equations.jl
@@ -0,0 +1,33 @@
+using Test
+
+import Distributions
+
+import ClimaAtmos
+
+function percentile_bounds_mean_norm_distributions(
+    low_percentile::FT,
+    high_percentile::FT,
+) where {FT <: Real}
+    gauss_int(x) = -exp(-x * x / 2) / sqrt(2 * pi)
+    xp_low = Distributions.quantile(Distributions.Normal(), low_percentile)
+    xp_high = Distributions.quantile(Distributions.Normal(), high_percentile)
+    return (gauss_int(xp_high) - gauss_int(xp_low)) / max(
+        Distributions.cdf(Distributions.Normal(), xp_high) -
+        Distributions.cdf(Distributions.Normal(), xp_low),
+        eps(FT),
+    )
+end
+
+@testset "Gauss quantile" begin
+    for p in 0.0:0.1:1.0
+        @test ClimaAtmos.guass_quantile(p) ≈
+              Distributions.quantile(Distributions.Normal(), p) rtol = 1e-3
+    end
+
+    for p in 0.0:0.1:1.0
+        for q in 0.0:0.1:1.0
+            @test ClimaAtmos.percentile_bounds_mean_norm(p, q) ≈
+                  percentile_bounds_mean_norm_distributions(p, q) rtol = 1e-3
+        end
+    end
+end