A structure-preserving discretization of the isoneutral diffusion of Redi and of the eddy induced skew diffusion scheme of Gent and McWilliams on unstructured grids is described. The discretization is based on a variational principle and its essential element is a construction of a discrete Hilbert space by means of admissible vector reconstructions. These reconstructions are applied to discretize the isoneutral slopes as the key component of the parametrizations. Through the discrete scalar product of the Hilbert space for the tracers we are able to represent isoneutral mixing tensor and eddy advection in weak form and to prove that the resulting discrete algorithm preserves essential physical properties of the continuous operator. The discretization follows in principle the established variational and skew-flux approach but implements these concepts via a natural extension of the discretization of the dynamical core of the unstructured grid ocean model ICON-O. The theoretical analysis is supplemented by numerical experiments that illustrate that the discretization indeed implements the functionality of isoneutral diffusion and skew diffusion as complementary parameterizations of unresolved mesoscale eddy effects.