Computational modelling of transport phenomena in packed beds is of interest in a wide range of engineering applications. We simulate flow in packed beds with similar conditions to heap leaching processes. The governing equations are the conservation of mass and Darcy's law. The primary variables are pressure and saturation. To discretize the governing equation, a control volume finite element method is used. The numerical scheme is implemented in Fluidity, a general purpose flow simulator that is equipped with anisotropic mesh adaptivity techniques to control local solution errors and increase computational accuracy. We demonstrate the application of the developed method to solve for incompressible flow of air and liquid within packed beds of mono-dispersed non-porous particles and compare our results with experimental measurements by Ilankoon and Neethling (2012). To study the effect of heterogeneity on flow regime in packed beds, an experiment is conducted on a packed bed consisting of two regions with different particle sizes. The non-uniform packed bed is also simulated numerically and our comparison shows very good agreement. Mesh adaptivity resolves saturation fronts in two regions of the column with high-resolution meshes and reduces the discretization error. This study serves as a validation for our numerical approach for prediction of hydrodynamics of heap leaching processes at large scales. (C) 2014 Elsevier Ltd. All rights reserved.