We consider the following repulsive-productive chemotaxis model: find u ≥ 0, the cell density, and v ≥ 0, the chemical concentration, satisfying1∂tu−Δu−∇⋅(u∇v)=0inΩ,t>0,∂tv−Δv+v=upinΩ,t>0,\\documentclass[12pt]{minimal}\t\t\t\t\\usepackage{amsmath}\t\t\t\t\\usepackage{wasysym}\t\t\t\t\\usepackage{amsfonts}\t\t\t\t\\usepackage{amssymb}\t\t\t\t\\usepackage{amsbsy}\t\t\t\t\\usepackage{mathrsfs}\t\t\t\t\\usepackage{upgreek}\t\t\t\t\\setlength{\\oddsidemargin}{-69pt}\t\t\t\t\\begin{document}$$ \\left\\{ \\begin{array}{l} \\partial_t u - {\\Delta} u - \\nabla\\cdot (u\\nabla v)=0 \\ \\ \\text{ in}\\ {\\Omega},\\ t>0,\\\\ \\partial_t v - {\\Delta} v + v = u^p \\ \\ { in}\\ {\\Omega},\\ t>0, \\end{array} \\right. $$\\end{document}with p ∈ (1, 2), {Omega }subseteq mathbb {R}^{d} 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.