This work is devoted to studying unconditionally energy stable and mass-conservative numerical schemes for the following repulsive-productive chemotaxis model: find $u \geq 0$, the cell density, and $v \geq 0$, the chemical concentration, such that \begin{equation*} \left \{ \begin {array} [c]{lll}\partial _t u - \Delta u - \nabla \cdot (u\nabla v)=0 \ \ \text {in}\ \Omega ,\ t>0,\\ \partial _t v - \Delta v + v = u \ \ \text {in}\ \Omega ,\ t>0, \end{array} \right . \end{equation*} in a bounded domain $\Omega \subseteq \mathbb {R}^d$, $d=2,3$. By using a regularization technique, we propose three fully discrete Finite Element (FE) approximations. The first one is a nonlinear approximation in the variables $(u,v)$; the second one is another nonlinear approximation obtained by introducing ${\boldsymbol \sigma }=\nabla v$ as an auxiliary variable; and the third one is a linear approximation constructed by mixing the regularization procedure with the energy quadratization technique, in which other auxiliary variables are introduced. In addition, we study the well-posedness of the numerical schemes, proving unconditional existence of solution, but conditional uniqueness (for the nonlinear schemes). Finally, we compare the behavior of such schemes throughout several numerical simulations and provide some conclusions.