A numerical algorithm for simulation of two‐phase flow in porous media is presented. The algorithm is based on a modified Picard linearization of the governing equations of flow, coupled with a lumped finite element approximation in space and dynamic time step control. Numerical results indicate that the algorithm produces solutions that are essentially mass conservative and oscillation free, even in the presence of steep infiltrating fronts. When the algorithm is applied to the case of air and water flow in unsaturated soils, numerical results confirm the conditions under which Richards's equation is valid. Numerical results also demonstrate the potential importance of air phase advection when considering contaminant transport in unsaturated soils. Comparison to several other numerical algorithms shows that the modified Picard approach offers robust, mass conservative solutions to the general equations that describe two‐phase flow in porous media.