Accretion onto black holes and other compact objects occurs across a wide range of scales. Despite the diversity in the physics involved, the variability shows a remarkable similarity in its properties. The theory of propagating fluctuations, in which random fluctuations within an accretion disk travel inwards and combine, has long been used to explain this variability. Recent numerical work has expanded on the extensive analytical literature but has been restricted to using the 1D diffusion equation for modelling the disk behaviour. In this talk, I present a novel numerical approach for 2D (vertically integrated), stochastically driven α-disk simulations, generalising existing 1D models. This is achieved through defining α as a function of a stochastic random variable β which is advected as a tracer field and stochastically evolved in time according to an Ornstein-Uhlenbeck process. I show that these simulations reproduce a wide range of observational features of AGN, including the linear rms-flux relationship, a broadband power spectrum which is best fit with a broken power-law and lags between emission in different energy bands. Additionally, there are two key differences between the existing theory of propagating fluctuations and our new 2D simulations. Firstly, we find that the presence of epicyclic motion in 2D (which cannot be captured within the 1D diffusion equation) has an important impact on local disk dynamics. Secondly, the light-curves are only log-normal when the disk is sufficiently thick, with thinner disks instead showing a normal distribution. We also find that these thinner disks are significantly less variable than thicker ones. I further show that the break frequency in the luminosity power spectrum is strongly dependent on the driving timescale of the stochastic evolution of the β parameter within the disk, providing a possible observational signature for probing the magnetorotational instability (MRI) dynamo.