Numerical simulations of unsteady groundwater flow in homogeneous and artificially generated heterogeneous geological formations have been presented. The heterogeneous structure of the geological patterns has been generated using the coupled chain Markov model developed by Elfeki and Dekking [2001]. Solution of the governing equations is achieved through the application of a finite difference approach to the partial differential equation of unsteady groundwater flow in horizontal plane with heterogeneous properties. It has been shown that global gradient (regional gradient) magnitude variability coupled with aquifer heterogeneity generates local directional and magnitude gradient variabilities. Increasing of the storage coefficient leads to smoothing of the aquifer response in terms of hydraulic head and Darcy's velocity. However, the smoothing effect is more pronounced in the hydraulic heads when compared with Darcy's velocity. In heterogeneous aquifer presented in this study, the aquifer response in terms of hydraulic head field and the lateral Darcy's velocity are in phase with the input time series, however the longitudinal Darcy's velocity is out of phase.