Today, the classical homogenised models for simulating excitable tissue are challenged by new mathematical frameworks that explicitly represent and resolve the geometry of extracellular and intracellular spaces and cellular membranes (EMI models). These mixed-dimensional models crucially enable an abstract representation of heterogeneous distribution of membrane ion channels and realistic cellular morphologies. EMI models typically predict the electrical properties of the tissue, with the underlying assumption that the ion concentrations are constant. Although this assumption is only an approximation, the resulting models still give accurate predictions of neuronal electrodynamics in many scenarios. They do however fail in capturing the numerous phenomena related to shifts in the extracellular ion concentrations. Here, we discuss an alternative approach to detailed modelling of excitable tissue that also accounts for ion concentrations and electrodiffusion. In particular, we consider a mathematical model describing the distribution and evolution of ion concentrations and potentials in a geometrically explicit setting. As a combination of the electroneutral Kirchhoff-Nernst-Planck (KNP) model and the EMI framework, we refer to this model as the KNP-EMI model. We introduce and numerically evaluate two new, finite element-based schemes for the KNP-EMI model, capable of efficiently and flexibly handling geometries of arbitrary dimension and arbitrary polynomial degree. The first scheme is based on a multi-dimensional formulation discretized with mortar finite elements, and the second is based on a single-dimensional formulation discretized with DG finite elements. We discuss stability and convergence properties of the two numerical schemes. Finally, we compare the electrical potentials predicted by the KNP-EMI and EMI models and study ephaptic coupling induced in an unmyelinated axon bundle and demonstrate how the KNP-EMI framework can give new insights in this setting.