In this paper we analyze a fully practical piecewise linear finite element approximation involving numerical integration, backward Euler time discretization, and possibly regularization of the following degenerate parabolic system arising in a model of reactive solute transport in porous media: find {u(x,t),v(x,t)} such that partial derivative(t)u+partial derivative(t)v - Delta u=f in Omega x (0,T] u=0 on partial derivative Omega x (0,T] partial derivative(t)v=k(phi(u)-v) in Omega x (0,T] u(.,0) = g(1)(.) v(.,0) = g(2)(.) in Omega subset of R(d), 1 less than or equal to d less than or equal to 3 for given data k epsilon R(+), f, g(1), g(2) and a monotonically increasing phi epsilon C-0(R) boolean AND C-1(-infinity, 0] boolean OR (0, infinity) satisfying phi(0) = 0, which is only locally Holder continuous with exponent p epsilon (0, 1) at the origin, e.g., phi(s) drop [s](p)(+). This lack of Lipschitz continuity at the origin limits the regularity of the unique solution {u, v} and leads to difficulties in the finite element error analysis.