Water infiltration and unsaturated flow through heterogeneous soil control the distribution of soil moisture in the vadose zone and the dynamics of groundwater recharge, providing the link between climate, biogeochemical soil processes and vegetation dynamics. Infiltration into dry soil is hydrodynamically unstable, leading to preferential flow through narrow wet regions (fingers). In this paper we use numerical simulation to study the interplay between fingering instabilities and soil heterogeneity during water infiltration. We consider soil with heterogeneous intrinsic permeability. Permeabilities are random, with point Gaussian statistics, and vary smoothly in space due to spatial correlation. The key research question is whether the presence of moderate or strong heterogeneity overwhelms the fingering instability, recovering the simple stable displacement patterns predicted by most simplified model of infiltration currently used in hydrological models from the Darcy to the basin scales. We perform detailed simulations of constant-rate infiltration into soils with isotropic and anisotropic intrinsic permeability fields. Our results demonstrate that soil heterogeneity does not suppress fingering instabilities, but it rather enhances its effect of preferential flow and channeling. Fingering patterns strongly depend on soil structure, in particular the correlation length and anisotropy of the permeability field. While the finger size and flow dynamics are only slightly controlled by correlation length in isotropic fields, layering leads to significant finger meandering and bulging, changing arrival times and wetting efficiencies. Fingering and soil heterogeneity need to be considered when upscaling the constitutive relationships of multiphase flow in porous media (relative permeability and water retention curve) from the finger to field and basin scales. While relative permeabilities remain unchanged upon upscaling for stable displacements, the inefficient wetting due to fingering leads to relative permeabilities at the field scale that are significantly different from those at the Darcy scale. These effective relative permeability functions also depend, although less strongly, on heterogeneity and soil structure.