Horizontal flow constructed wetlands are engineered systems capable of eliminating a wide range of pollutants from the aquatic environment. Nevertheless, poor hydrodynamic behavior is commonly found resulting in preferential pathways and variations in both (i) the hydraulic residence time distribution (HRTD) and, consequently, (ii) the wetland's treatment efficiency. The aim of this work was to outline a methodology for wetland design that accounts for the effect of heterogeneous hydraulic properties of the porous substrate on the HRTD and treatment efficiency. Biodegradation of benzene was used to illustrate the influence of hydraulic conductivity heterogeneity on wetland efficiency. Random, spatially correlated hydraulic conductivity fields following a log-normal distribution were generated and then introduced in a subsurface flow numerical model. The results showed that the variance of the distribution and the correlation length in the longitudinal direction are key indicators of the extent of heterogeneity. A reduction of the mean hydraulic residence time was observed as the extent of heterogeneity increased, while the HRTD became broader with increased skewness. At the same time, substrate heterogeneity induced preferential flow paths within the wetland bed resulting in variations of the benzene treatment efficiency. Further to this it was observed that the distribution of biomass within the porous bed became heterogeneous, rising questions on the representativeness of sampling. It was concluded that traditional methods for wetland design based on assumptions such as a homogeneous porous medium and plug flow are not reliable. The alternative design methodology presented here is based on the incorporation of heterogeneity directly during the design phase. The same methodology can also be used to optimize existing systems, where the HRTD has been characterized with tracer experiments.