The self-consistent modeling of runaway electrons (REs) and their interaction with non-linear MHD is a difficult problem owing to the vast separation of time scales between the particle transport and the growth rate of the instabilities, as well as the 3D geometry needed to capture these dynamics. Because of this REs are often treated in a fluid framework to facilitate the crossing of longer time scales. Kinetic effects such as the drift orbits of the particles and changes in the force balance may however affect the growth and saturation of instabilities as well as the rate at which REs are lost during the stochastic phase. To better study these interactions a 3D global non-linear full-f particle-in-cell model for REs has been implemented in JOREK which retains the phase space information of the REs and can resolve the full gyro-orbit of the particles. In this work we introduce the hybrid kinetic-MHD scheme and discuss the steps taken to benchmark the model, both in 2D and 3D configurations. Results from initial studies of a beam termination scenario in JET are then shown, with particular emphasis on differences seen when varying the RE energy. Finally, we present ongoing work aiming to adopt more realistic phase space distributions in the simulations through dedicated studies using the DREAM code.