Porous media play an increasing role in industry. Some important areas of application include: chemical conversion, heat transfer, and particle filtration. Having a detailed understanding of such processes (e.g., to enhance process efficiency) often requires the description of transport phenomena on the level of the pore scale. Nowadays, one can incorporate complex pore structures obtained using imaging techniques, e.g., computer tomography, directly into modeling strategies. How to do this efficiently and accurately is the topic of the present research. We present the development and application of a volume-penalizing immersed boundary (IB) method  for the simulation of fluid transport in geometrically complex porous media. In this methodology, the transport in the multiphase (fluid and solid) domain is governed by the set of equations describing the motion of a fluid. These equations, i.e., the incompressible Navier-stokes equations, are extended with a source term (or forcing function) to approximate the net momentum exchange between the fluid and the solid phase across their mutual interfaces. The forcing function is essentially responsible for approximating the no-slip condition (vanishing velocity components) on solid boundaries. By combining the forcing strategy with a Cartesian grid representation of the multiphase domain a simple yet powerful computational method is obtained for modeling fluid transport in geometrically complex domains. Using a Cartesian grid representation allows the avoidance of body-conforming grids, as they are difficult and time-consuming to generate for geometrically complex domains. Also, a Cartesian grid has the benefit that highly efficient solvers can be used. In the present case, spatial discretization is performed using a finite-volume method that preserves the symmetry properties of the convective and diffusive operators in the governing equations . In fig. 1 we show a tomographic representation of a cross-sectional plane through a model porous medium (left) and its Cartesian grid representation (right) using (for illustration purposes) a coarse spatial resolution of (64)2 grid points. A three-dimensional porous medium is reconstructed by placing several cross-sectional planes in the third dimension. In fig. 2 we present simulation results for the out-of-plane velocity component in a representation using (256)2 grid points per cross-section. On the basis of such detailed simulations one can predict macroscopic transport quantities such as the permeability to support efficient macroscopic engineering transport models in porous media. Topics of ongoing research include: improvements in accuracy and order of the IB method, and an extension to incorporate heat and mass transfer. These topics will pave the way for the development of a reliable simulation tool capable of realistic geometries and complex physical processes, and will be addressed in the proposed contribution.  R. Mittal and G. Laccarino, Immersed boundary methods, Ann. Rev. Fluid Mech., 37, 239--261 (2005).  R. W. C. P. Verstappen and A. E. P. Veldman, Symmetry-preserving discretization of turbulent flow, J. Comput. Phys., 187, 343-368 (2003).