We present a discrete Helmholtz–Hodge decomposition for H(div)-conforming Brezzi–Douglas–Marini (BDM) finite elements on triangulated surfaces of arbitrary topology. The divergence-free BDM subspace is split L2-orthogonally into rotated gradients of a continuous streamfunction space and a finite-dimensional space of discrete harmonic fields whose dimension equals the first Betti number of the surface. Consequently, any incompressible flow discretized on this subspace can be reformulated with a scalar streamfunction and finitely many harmonic coefficients as the only unknowns. This eliminates the pressure and the saddle-point structure while ensuring exact tangentiality, pointwise divergence-freeness, and pressure-robustness. We present a randomized algorithm for constructing the harmonic basis and discuss implementation aspects including hybridization, efficient treatment of the harmonic unknowns, and pressure reconstruction. Numerical experiments for unsteady surface Navier–Stokes equations on a trefoil knot and a multiply-connected sculpture surface demonstrate the method and illustrate the physical role of the harmonic velocity component.