Higher-order finite volume flux operators for transport algorithms used within Runge-Kutta time integration schemes on irregular Voronoi (hexagonal) meshes are proposed and tested. These operators are generalizations of 3rd- and 4th-order operators currently used in atmospheric models employing regular, orthogonal rectangular meshes. 2D least-squares-fit polynomials are used to evaluate the higher-order spatial derivatives needed to cancel the leading order truncation error terms of the standard 2nd-order centered formulation. Positive definite or monotonic behavior is achieved by applying an appropriate limiter during the final Runge-Kutta stage within a given time step. The 3rd- and 4th-order formulations are evaluated using standard transport tests on the sphere. The new schemes are more accurate and significantly more efficient than the standard second-order scheme and other schemes in the literature we have examined. The 3rd-order formulation is equivalent to the 4th-order formulation plus an additional diffusion term that is proportional to the Courant number. An optimal value for the coefficient scaling this diffusion term is chosen based on qualitative evaluation of the test results. Improvements using the higher-order scheme in place of the traditional second-order centered approach are illustrated within 3D unstable baroclinic wave simulations produced using two global nonhydrostatic models employing spherical Voronoi meshes.