diff --git a/docs/src/supportmatrix.md b/docs/src/supportmatrix.md
index 82fcd5f3..4265c06a 100644
--- a/docs/src/supportmatrix.md
+++ b/docs/src/supportmatrix.md
@@ -1,58 +1,61 @@
 # Support Matrix
 
-While this library aims to support all possible integration rules and [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl)
-geometry types, some combinations are ill-suited and some others are simply not yet
-implemented. The following Support Matrix aims to capture the current development state of
-all geometry/rule combinations. Entries with a green check mark are fully supported
-and have passing unit tests that provide some confidence they produce accurate results.
+This library aims to enable users to calculate the value of integrals over all [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl)
+geometry types using an array of numerical integration rules and techniques. However, some
+combinations of geometry types and integration rules are ill-suited, and some others are simply
+not yet implemented. The following Support Matrix captures the current state of support for
+all geometry/rule combinations. Entries with a green check mark are fully supported and pass
+unit tests designed to check for accuracy.
 
 In general, Gauss-Kronrod integration rules are recommended (and the default) for geometries
-with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. Gauss-Kronrod
-rules can also be applied to some geometries with more dimensions by nesting multiple
-integration solves, but this is inefficient. These Gauss-Kronrod rules are supported (but
-not recommended) for surface-like geometries, but not for volume-like geometries. For
-geometries with more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive
-Cubature integration rules are recommended (and the default).
+with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. For geometries with
+more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive Cubature rules are
+recommended (and the default).
+
+While it is possible to apply nested Gauss-Kronrod rules to numerically integrate geometries
+with more than one parametric dimension, this produces results that are strictly inferior to
+using an equivalent H-Adaptive Cubature rule, so support for this usage is not recommended.
 
 | Symbol | Support Level |
 |--------|---------|
-| ✅ | Supported, passes unit tests |
+| ✅ | Supported |
 | 🎗️ | Planned to support in the future |
+| ⚠️ | Deprecated |
 | 🛑 | Not supported |
 
 | `Meshes.Geometry` | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
 |----------|----------------|---------------|---------------------|
-| `Ball` in `𝔼{2}` | ✅ | ✅ | ✅ |
+| `Ball` in `𝔼{2}` | ✅ | ⚠️ | ✅ |
 | `Ball` in `𝔼{3}` | ✅ | 🛑 | ✅ |
 | `BezierCurve` | ✅ | ✅ | ✅ |
 | `Box` in `𝔼{1}` | ✅ | ✅ | ✅ |
-| `Box` in `𝔼{2}` | ✅ | ✅ | ✅ |
+| `Box` in `𝔼{2}` | ✅ | ⚠️ | ✅ |
 | `Box` in `𝔼{≥3}` | ✅ | 🛑 | ✅ |
 | `Circle` | ✅ | ✅ | ✅ |
-| `Cone` | ✅ | ✅ | ✅ |
-| `ConeSurface` | ✅ | ✅ | ✅ |
+| `Cone` | ✅ | 🛑 | ✅ |
+| `ConeSurface` | ✅ | ⚠️ | ✅ |
 | `Cylinder` | ✅ | 🛑 | ✅ |
-| `CylinderSurface` | ✅ | ✅ | ✅ |
-| `Disk` | ✅ | ✅ | ✅ |
+| `CylinderSurface` | ✅ | ⚠️ | ✅ |
+| `Disk` | ✅ | ⚠️ | ✅ |
 | `Ellipsoid` | ✅ | ✅ | ✅ |
 | `Frustum` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
-| `FrustumSurface` | ✅ | ✅ | ✅ |
+| `FrustumSurface` | ✅ | ⚠️ | ✅ |
 | `Hexahedron` | ✅ | ✅ | ✅ |
 | `Line` | ✅ | ✅ | ✅ |
-| `ParaboloidSurface` | ✅ | ✅ | ✅ |
+| `ParaboloidSurface` | ✅ | ⚠️ | ✅ |
 | `ParametrizedCurve` | ✅ | ✅ | ✅ |
 | `Plane` | ✅ | ✅ | ✅ |
 | `Polyarea` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
 | `Pyramid` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
-| `Quadrangle` | ✅ | ✅ | ✅ |
+| `Quadrangle` | ✅ | ⚠️ | ✅ |
 | `Ray` | ✅ | ✅ | ✅ |
 | `Ring` | ✅ | ✅ | ✅ |
 | `Rope` | ✅ | ✅ | ✅ |
 | `Segment` | ✅ | ✅ | ✅ |
 | `SimpleMesh` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) |
 | `Sphere` in `𝔼{2}` | ✅ | ✅ | ✅ |
-| `Sphere` in `𝔼{3}` | ✅ | ✅ | ✅ |
+| `Sphere` in `𝔼{3}` | ✅ | ⚠️ | ✅ |
 | `Tetrahedron` in `𝔼{3}` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | ✅ | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) |
 | `Triangle` | ✅ | ✅ | ✅ |
-| `Torus` | ✅ | ✅ | ✅ |
+| `Torus` | ✅ | ⚠️ | ✅ |
 | `Wedge` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
diff --git a/src/integral.jl b/src/integral.jl
index 0fee71fb..de69db5d 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -33,7 +33,7 @@ function integral(
 end
 
 ################################################################################
-#                    Generalized (n-Dimensional) Worker Methods
+#                            Integral Workers
 ################################################################################
 
 # GaussKronrod
@@ -41,14 +41,23 @@ function _integral(
         f,
         geometry,
         rule::GaussKronrod;
-        kwargs...
-)
-    # Pass through to dim-specific workers in next section of this file
+        FP::Type{T} = Float64,
+        diff_method::DM = _default_method(geometry)
+) where {DM <: DifferentiationMethod, T <: AbstractFloat}
+    # Implementation depends on number of parametric dimensions over which to integrate
     N = Meshes.paramdim(geometry)
     if N == 1
-        return _integral_gk_1d(f, geometry, rule; kwargs...)
+        integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
+        return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
     elseif N == 2
-        return _integral_gk_2d(f, geometry, rule; kwargs...)
+        # Issue deprecation warning
+        Base.depwarn("Use `HAdaptiveCubature` instead of \
+                     `GaussKronrod` for surfaces.", :integral)
+
+        # Nested integration
+        integrand2D(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method)
+        ∫₁(v) = QuadGK.quadgk(u -> integrand2D(u, v), zero(FP), one(FP); rule.kwargs...)[1]
+        return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
     else
         _error_unsupported_combination("geometry with more than two parametric dimensions",
             "GaussKronrod")
@@ -71,14 +80,14 @@ function _integral(
     nodes = Iterators.product(ntuple(Returns(xs), N)...)
 
     # Domain transformation: x [-1,1] ↦ t [0,1]
-    t(x) = FP(1 // 2) * x + FP(1 // 2)
+    t(x) = (1 // 2) * x + (1 // 2)
 
     function integrand((weights, nodes))
         ts = t.(nodes)
         prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method)
     end
 
-    return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes))
+    return (1 // (2^N)) .* sum(integrand, zip(weights, nodes))
 end
 
 # HAdaptiveCubature
@@ -105,30 +114,3 @@ function _integral(
     # Reapply units
     return value .* integrandunits
 end
-
-################################################################################
-#                    Specialized GaussKronrod Methods
-################################################################################
-
-function _integral_gk_1d(
-        f,
-        geometry,
-        rule::GaussKronrod;
-        FP::Type{T} = Float64,
-        diff_method::DM = _default_method(geometry)
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
-    return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
-end
-
-function _integral_gk_2d(
-        f,
-        geometry2d,
-        rule::GaussKronrod;
-        FP::Type{T} = Float64,
-        diff_method::DM = _default_method(geometry2d)
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method)
-    ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1]
-    return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
-end
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index 7b0f44ab..c3540cbc 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -43,12 +43,12 @@ function integral(
     xs, ws = _gausslegendre(FP, rule.n)
 
     # Change of variables: x [-1,1] ↦ t [0,1]
-    t(x) = FP(1 // 2) * x + FP(1 // 2)
+    t(x) = (1 // 2) * x + (1 // 2)
     point(x) = curve(t(x), alg)
     integrand(x) = f(point(x)) * differential(curve, (t(x),), diff_method)
 
     # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length]
-    return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
+    return (1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
 end
 
 function integral(
diff --git a/test/Project.toml b/test/Project.toml
index 81203128..e38225db 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,7 +4,6 @@ CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
-QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
@@ -17,7 +16,6 @@ CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
 ExplicitImports = "1.6.0"
 LinearAlgebra = "1"
 Meshes = "0.50, 0.51, 0.52"
-QuadGK = "2.1.1"
 SpecialFunctions = "2"
 TestItemRunner = "1"
 TestItems = "1"