Abstract We consider the following repulsive-productive chemotaxis model: find u ≥ 0, the cell density, and v ≥ 0, the chemical concentration, satisfying $$ \left\{ \begin{array}{l} \partial_t u - {\Delta} u - \nabla\cdot (u\nabla v)=0 \ \ \text{ in}\ {\Omega},\ t&gt;0,\\ \partial_t v - {\Delta} v + v = u^p \ \ { in}\ {\Omega},\ t&gt;0, \end{array} \right. $$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mfenced> <mml:mrow> <mml:mtable> <mml:mtr> <mml:mtd> <mml:msub> <mml:mrow> <mml:mi>∂</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>t</mml:mi> </mml:mrow> </mml:msub> <mml:mi>u</mml:mi> <mml:mo>−</mml:mo> <mml:mi>Δ</mml:mi> <mml:mi>u</mml:mi> <mml:mo>−</mml:mo> <mml:mo>∇</mml:mo> <mml:mo>⋅</mml:mo> <mml:mo>(</mml:mo> <mml:mi>u</mml:mi> <mml:mo>∇</mml:mo> <mml:mi>v</mml:mi> <mml:mo>)</mml:mo> <mml:mo>=</mml:mo> <mml:mn>0</mml:mn> <mml:mspace /> <mml:mspace /> <mml:mspace /> <mml:mtext>in</mml:mtext> <mml:mspace /> <mml:mi>Ω</mml:mi> <mml:mo>,</mml:mo> <mml:mspace /> <mml:mi>t</mml:mi> <mml:mo>&gt;</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> </mml:mtd> </mml:mtr> <mml:mtr> <mml:mtd> <mml:msub> <mml:mrow> <mml:mi>∂</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>t</mml:mi> </mml:mrow> </mml:msub> <mml:mi>v</mml:mi> <mml:mo>−</mml:mo> <mml:mi>Δ</mml:mi> <mml:mi>v</mml:mi> <mml:mo>+</mml:mo> <mml:mi>v</mml:mi> <mml:mo>=</mml:mo> <mml:msup> <mml:mrow> <mml:mi>u</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>p</mml:mi> </mml:mrow> </mml:msup> <mml:mspace /> <mml:mspace /> <mml:mi>i</mml:mi> <mml:mi>n</mml:mi> <mml:mspace /> <mml:mi>Ω</mml:mi> <mml:mo>,</mml:mo> <mml:mspace /> <mml:mi>t</mml:mi> <mml:mo>&gt;</mml:mo> <mml:mn>0</mml:mn> <mml:mo>,</mml:mo> </mml:mtd> </mml:mtr> </mml:mtable> </mml:mrow> </mml:mfenced> </mml:math> with p ∈ (1, 2), ${\Omega }\subseteq \mathbb {R}^{d}$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mi>Ω</mml:mi> <mml:mo>⊆</mml:mo> <mml:msup> <mml:mrow> <mml:mi>ℝ</mml:mi> </mml:mrow> <mml:mrow> <mml:mi>d</mml:mi> </mml:mrow> </mml:msup> </mml:math> a bounded domain ( d = 1, 2, 3), endowed with non-flux boundary conditions. By using a regularization technique, we prove the existence of global in time weak solutions of (1) which is regular and unique for d = 1, 2. Moreover, we propose two fully discrete Finite Element (FE) nonlinear schemes, the first one defined in the variables ( u , v ) under structured meshes, and the second one by using the auxiliary variable σ = ∇ v and defined in general meshes. We prove some unconditional properties for both schemes, such as mass-conservation, solvability, energy-stability and approximated positivity. Finally, we compare the behavior of these schemes with respect to the classical FE backward Euler scheme throughout several numerical simulations and give some conclusions.