Bipolar resistive switching (BRS) cells based on the valence change mechanism show great potential to enable the design of future non-volatile memory, logic and neuromorphic circuits and architectures. To study these circuits and architectures, accurate compact models are needed, which showcase the most important physical characteristics and lead to their specific experimental behavior. If BRS cells are to be used for computation-in-memory or for neuromorphic computing, their dynamical behavior has to be modeled with special consideration of switching times in SET and RESET. For any realistic assessment, variability has to be considered additionally. This study shows that by extending an existing compact model, which by itself is able to reproduce many different experiments on device behavior critical for the anticipated device purposes, variability found in experimental measurements can be reproduced for important device characteristics such as I-V characteristics, endurance behavior and most significantly the SET and RESET kinetics. Furthermore, this enables the study of spatial and temporal variability and its impact on the circuit and system level.