The Brinkman equation is used to model the isothermal flow of Newtonian fluids in a highly permeable porous medium between two parallel plates. Due to the multi-scale nature of flow regimes represented by the Brinkman equation, the standard Galerkin method requires excessive mesh refinement to achieve stable and accurate results. To avoid this, a residual free method is used to derive appropriate bubble functions for inclusion in the finite element discretisations. It is shown that by using bubble enriched shape functions the standard Galerkin method can yield accurate and stable solutions for multi-scale problems. In this paper the performances of the third and fifth order bubble functions used in conjunction with bilinear Lagrangian elements to solve the Brinkman equation via a penalty finite element scheme, are reported. (c) 2005 Published by Elsevier Ltd.