The coupled equations for flow in unsaturated soil as proposed by Beliaev and Hassanizadeh  are described. These equations account for mechanism of dynamic capillary pressure via a first order relaxation function. A form of the relaxation function for the dynamic capillary pressure-saturation relation is proposed based on physical reasoning and a semi-analytical solution to the flow equations. A mass conservative and computational efficient numerical solution to the coupled equations in two space dimensions is derived and applied to the simulation of gravity-driven unstable flow. Simulated fingers have all the morphological features of fingers observed in laboratory experiments. The results demonstrate that the dynamic capillary pressure mechanism causes initial destabilization of the flow, while the mechanism of capillary hysteresis leads to finger persistence.