diff --git a/src/quadrature.jl b/src/quadrature.jl
index 9f1d9613..4bc5f201 100644
--- a/src/quadrature.jl
+++ b/src/quadrature.jl
@@ -372,3 +372,89 @@ function _curvature(f, Q)
     end
     return curv
 end
+
+"""
+    tangential_gradient_matrix(Q::Quadrature{N,T})
+
+Return a sparse matrix `G` of size `(length(Q), length(Q))` with `SVector{N,T}`
+entries such that `G * u` computes the surface gradient `∇_Γ u` at each
+quadrature node for scalar values `u`.
+
+The surface gradient is computed by locally interpolating `u` in parameter space
+using the Lagrange basis on quadrature nodes, differentiating, and applying the
+chain rule through the element parametrization:
+
+```math
+∇_Γ u = J (Jᵀ J)^{-1} ∇_{\\hat{x}} \\tilde{u}
+```
+
+where `J` is the Jacobian of the element map and `∇_{\\hat{x}} \\tilde{u}` is
+the gradient of the interpolant in reference coordinates.
+"""
+function tangential_gradient_matrix(Q::Quadrature{N, T}) where {N, T}
+    msh = mesh(Q)
+    Is = Int[]
+    Js = Int[]
+    Vs = SVector{N, T}[]
+    ntotal = sum(
+        E -> length(qcoords(quadrature_rule(Q, E)))^2 * size(Q.etype2qtags[E], 2),
+        element_types(msh)
+    )
+    sizehint!(Is, ntotal)
+    sizehint!(Js, ntotal)
+    sizehint!(Vs, ntotal)
+    for (E, qtags_mat) in Q.etype2qtags
+        _tangential_gradient_kernel!(Is, Js, Vs, Q, elements(msh, E), qtags_mat)
+    end
+    return sparse(Is, Js, Vs, length(Q), length(Q))
+end
+
+"""
+    surface_gradient(u::AbstractVector, Q::Quadrature)
+
+Compute the surface gradient `∇_Γ u` at each quadrature node, returning a
+`Vector{SVector{N,T}}`.
+
+See also: [`tangential_gradient_matrix`](@ref)
+"""
+function surface_gradient(u::AbstractVector, Q::Quadrature)
+    return tangential_gradient_matrix(Q) * u
+end
+
+@noinline function _tangential_gradient_kernel!(
+        Is,
+        Js,
+        Vs,
+        Q::Quadrature{N, T},
+        els::AbstractVector{E},
+        qtags_mat::Matrix{Int},
+    ) where {N, T, E}
+    M = geometric_dimension(domain(E))
+    qrule = quadrature_rule(Q, E)
+    x̂_nodes = qcoords(qrule)
+    nq = length(x̂_nodes)
+    L = lagrange_basis(qrule)
+    # precompute derivative of Lagrange basis at each reference node (shared across elements)
+    # dL_rows[q][j] is the SVector{M,T} gradient of Lⱼ at x̂_q
+    dL_rows = map(x̂_nodes) do x̂
+        dL = ForwardDiff.jacobian(L, x̂)
+        ntuple(j -> SVector{M, T}(ntuple(k -> T(dL[j, k]), M)), nq)
+    end
+    for n in 1:size(qtags_mat, 2)
+        el = els[n]
+        qtags = view(qtags_mat, :, n)
+        for q in 1:nq
+            J_q = SMatrix{N, M, T}(jacobian(el, x̂_nodes[q]))
+            invA_q = inv(J_q' * J_q)
+            B_q = J_q * invA_q
+            i_global = qtags[q]
+            for j in 1:nq
+                coeff = B_q * dL_rows[q][j]
+                push!(Is, i_global)
+                push!(Js, qtags[j])
+                push!(Vs, coeff)
+            end
+        end
+    end
+    return nothing
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 61b6439e..4ff22c58 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,6 +21,8 @@ using Aqua
 
 @safetestset "Quadrature" include("quadrature_test.jl")
 
+@safetestset "Surface gradient" include("surface_gradient_test.jl")
+
 @safetestset "Normal orientation" include("normal_orientation_test.jl")
 
 @safetestset "Kernels" include("kernels_test.jl")
diff --git a/test/surface_gradient_test.jl b/test/surface_gradient_test.jl
new file mode 100644
index 00000000..8d701949
--- /dev/null
+++ b/test/surface_gradient_test.jl
@@ -0,0 +1,69 @@
+using Test
+using LinearAlgebra
+using StaticArrays
+using Inti
+
+@testset "2D circle, linear function" begin
+    Inti.clear_entities!()
+    Γ = Inti.parametric_curve(t -> SVector(cos(t), sin(t)), 0, 2π) |> Inti.Domain
+    quad = Inti.Quadrature(Γ; meshsize = 0.1, qorder = 5)
+    # f(x,y) = x on unit circle
+    u = map(q -> q.coords[1], quad)
+    grad = Inti.surface_gradient(u, quad)
+    # analytic: ∇_Γ(x) = (y², -xy) on unit circle
+    grad_exact = map(quad) do q
+        x, y = q.coords
+        SVector(y^2, -y * x)
+    end
+    @test norm(grad .- grad_exact, Inf) < 1.0e-2
+end
+
+@testset "3D ellipsoid, linear function" begin
+    Inti.clear_entities!()
+    Ω = Inti.GeometricEntity("ellipsoid") |> Inti.Domain
+    Γ = Inti.external_boundary(Ω)
+    quad = Inti.Quadrature(Γ; meshsize = 0.2, qorder = 4)
+    # f(x) = a ⋅ x with known a
+    a = SVector(1.0, 0.5, -0.3)
+    u = map(q -> dot(a, q.coords), quad)
+    grad = Inti.surface_gradient(u, quad)
+    # analytic: ∇_Γ f = a - (a⋅n̂)n̂
+    grad_exact = map(quad) do q
+        n̂ = q.normal
+        a - dot(a, n̂) * n̂
+    end
+    err = maximum(norm(gc - ga) for (gc, ga) in zip(grad, grad_exact))
+    @test err < 5.0e-2
+end
+
+@testset "matrix vs direct consistency" begin
+    Inti.clear_entities!()
+    Γ = Inti.parametric_curve(t -> SVector(cos(t), sin(t)), 0, 2π) |> Inti.Domain
+    quad = Inti.Quadrature(Γ; meshsize = 0.1, qorder = 3)
+    u = map(q -> sin(q.coords[1]), quad)
+    G = Inti.tangential_gradient_matrix(quad)
+    grad1 = G * u
+    grad2 = Inti.surface_gradient(u, quad)
+    @test grad1 == grad2
+end
+
+@testset "convergence under refinement (2D circle)" begin
+    # f(x,y) = x² on unit circle
+    # ∇_Γ f = (2x - 2x³, -2x²y)
+    errs = Float64[]
+    for h in [0.4, 0.2, 0.1, 0.05]
+        Inti.clear_entities!()
+        Γ = Inti.parametric_curve(t -> SVector(cos(t), sin(t)), 0, 2π) |> Inti.Domain
+        quad = Inti.Quadrature(Γ; meshsize = h, qorder = 5)
+        u = map(q -> q.coords[1]^2, quad)
+        grad = Inti.surface_gradient(u, quad)
+        grad_exact = map(quad) do q
+            x, y = q.coords
+            SVector(2x - 2x^3, -2x^2 * y)
+        end
+        push!(errs, maximum(norm(gc - ga) for (gc, ga) in zip(grad, grad_exact)))
+    end
+    # check monotonic convergence
+    ratios = errs[1:(end - 1)] ./ errs[2:end]
+    @test all(r > 1.5 for r in ratios)
+end