Modified Newtonian Dynamics (MOND) can account for a variety of phenomena on galactic scales without the need for dark matter, but it cannot fully explain the mass contained in galaxy clusters. We explore two possible solutions to this problem: relativistic extensions of MOND, and FEM simulations of the apparent matter distribution in clusters, utilising the non-linear AQUAL formulation to analyse non spherically symmetric systems. We consider Covariant Emergent Gravity, and we show that the theory is inconsistent with the original formulation of Emergent Gravity. Moreover, we show that either the theory is incompatible with observations, or it presents grave theoretical difficulties. We then suggest that a covariant formulation of EG can be obtained through a Generalised Einstein Aether (GEA) theory, which is capable of retrieving the MOND PDE, and is, at the same time, consistent with observational constraints. Regarding the FEM simulations of the apparent matter distribution in galaxy clusters, we construct a sample of 15 clusters from the catalogs of Reiprich and Abell for the baryonic mass distribution present in galaxy clusters. We choose FEM for its ability to treat the combination of continuous and discrete mass distributions without the need for smoothing. This is necessary, as in MOND we cannot apply the principle of superposition or the weak lensing formalism. We then utilise the FeniCS software package to study the properties of these clusters. We simulate each cluster with elements up to degree 3, thanks to a speedup by a factor of ∼ 100 obtained by the use of local mesh refinement for the serial case. In addition, we run the code in parallel on a single-threaded 8-core CPU, achieving near optimal weak scaling for the regime of interest. For the mass distribution in the galaxy clusters, we analyse the distribution of baryons, Phantom Dark Matter (PDM) and apparent mass both close to the core and around each galaxy. We find that the PDM tends to clump around the galaxies, regardless of the gas to galaxy mass ratio. Moreover, we show that both the apparent mass and the PDM can exhibit negative masses as predicted by Milgrom. Our observations on the density of PDM around the galaxies match recent observations of small scale weak lensing. In addition, our results for negative mass distributions provide an opportunity to test a prediction of MOND that can never be replicated in the dark-matter paradigm, and shed light on the properties of non-spherically symmetric mass distributions, that have, up to now, not been studied in the literature for the fully non-linear case.