This paper presents a detailed analysis of a numerical solution of three-dimensional transient Darcy flow. The numerical solution of the governing parabolic partial differential equations is obtained by using stabilized mixed Galerkin method and backward Euler method for the discretization of space and time, respectively. The resulting well-posed system of algebraic equations is subsequently solved using conjugate gradient method. The proposed model is validated against Mongan's analytical model for underground water flow using a set of hexahedral and tetrahedral meshes. The model is used to analyze the transient behavior by simulating the Darcy flow through homogeneous and heterogeneous as well as isotropic and anisotropic media. For large meshes, a parallel algorithm of the transient Darcy flow is also developed for shared memory architecture using OpenMP library. For structured meshes, a speedup of over 22 is obtained on dual AMD Opteron processors. The proposed numerical method for transient Darcy flow offers stability, ease of implementation in higher dimensions and parallel solution for large and complex geometry using standard finite element spaces.