Focus is on comparing stochastic, process-based and deterministic methods for modelling heterogeneity in hydraulic properties of fluvial geothermal reservoirs. Models are considered a generalized representation of a fluvial sequence in the upper part of the Gassum Formation in northern Denmark. Two ensemble realizations of process-based and stochastic heterogeneity were simulated using the state-of-the-art numerical modelling software, Delft Advanced Research Terra Simulator (DARTS), to assess differences on a statistically relevant sample size. Simulator settings were optimized to achieve two orders of magnitude improvement in simulation time. Our general findings show that the stochastic and process-based methods produce nearly identical results in terms of predicted breakthrough time and production temperature decline for high net-to-gross ratios (N/G). Simple homogenous and layered models overestimate breakthrough and underestimate temperature decline. More complex representation of facies in process-based models show smaller variance in results and stay within the limits of ensemble runs compared to simpler facies representation. Results indicate that a multivariate Gaussian based stochastic representation of heterogeneity provides comparable thermal response to a process-based model in fluvial systems of similar quality.