In this paper, we consider a version of the mathematical model introduced in (Wang et al. in Commun. Nonlinear Sci. Numer. Simul. 42:571–584, 2017) to describe the interaction between vegetation and soil water in arid environments. The model corresponds to a nonlinear parabolic coupled system of partial differential equations, with non-flux boundary conditions, which incorporates, in addition to the natural diffusion of water and plants, a cross-diffusion term given by the hydraulic diffusivity due to the suction of water by the roots. The model also considers a monotonously decreasing vegetation death rate capturing the infiltration feedback between plants and ground water. We first prove the existence and uniqueness of global solutions in a large class of initial data, allowing non-regular ones. These solutions are in a mild setting and under additional regularity assumptions on the initial data and the domain, they are classical. Second, we propose a fully discrete numerical scheme, based on a semi-implicit Euler discretization in time and finite element discretization (with “mass-lumping”) in space, for approximating the solutions of the continuous model. We prove the well-posedness of the numerical scheme and some qualitative properties of the discrete solutions including, positivity, uniform weak and strong estimates, convergence towards strong solutions and optimal error estimates. Finally, we present some numerical experiments in order to showcase the good behavior of the numerical scheme including the formation of Turing patterns, as well as to validate the convergence order in the error estimates obtained in the theoretical analysis.