Constant-order derivative (COD) time fractional diffusion models have been successfully employed to describe transport processes where the rate of diffusion or the diffusion behavior differs from the classical Brownian motion (i.e. anomalous diffusion). More so, the Variable-order derivative (VOD) time fractional diffusion models have been recently acknowledged as an alternative and more rigorous approach to handle time-dependent anomalous diffusion behavior observed in highly heterogeneous fractured porous media.This study exhibits two finite difference approximations based on control volume finite difference approximations to handle a VOD time fractional diffusion mathematical model describing fluid flow in a heterogeneous porous medium. Subsequently, the numerical models are validated against the exact solution of the simplified COD time fractional diffusion mathematical model with constant rock and fluid properties. The well-established COD time fractional diffusion mathematical model is a special case of the herein presented VOD time fractional diffusion mathematical model. Furthermore, a modified incremental balance check and a modified wellbore flow model are presented to account for the variable order fractional exponent. The results from the numerical simulations are presented to reveal the peculiar features of the VOD time fractional diffusion model. This study would herald more applications of the VOD time fractional diffusion models in numerical modelling of fluid flow in porous media of fractal geometry or heterogeneous media.