diff --git a/Project.toml b/Project.toml
index 97ab802e38c..071b4948243 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/src/prognostic_equations/edmfx_boundary_condition.jl b/src/prognostic_equations/edmfx_boundary_condition.jl
index 8f177c1bbdc..76a959692b8 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)
+    a = FT(0.147)
+    ln = log(1 - x^2)
+    term1 = (2 / (π * a) + ln / 2)
+    term2 = ln / a
+    term3 = sqrt(term1^2 - term2)
+
+    return sign(x) * sqrt(term3 - term1)
+end
+
+function guass_quantile(p::FT)
+    return sqrt(2) * approximate_inverf(2 * (1 - p))
+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