A method for incorporating surface heterogeneity into models of wind erosion and dust emission is presented. Parameters that control threshold shear velocity (bare-soil threshold shear velocity, vegetation cover, and vegetation shape) are parameterized as modified-normal distributions with defined means and coefficients of variation. A parameterized bootstrap (Monte Carlo) approach is used to estimate histograms of threshold shear velocities for vegetated surfaces using the model of Raupach et al. (1993). Horizontal mass flux is calculated from threshold shear velocity histograms and a wind record. Vertical mass flux (dust emission) can be inferred from horizontal mass flux. The results from these simulations indicate that wind erosion and dust emission are dominated by the few instances with unusually high surface erodibility (e.g., low threshold shear velocity). This finding is supported by observations of dust emission hot spots in landscapes undergoing dust emission. As a result, wind erosion and dust emission flux rates are much higher in simulations with variability than they are in nonvariable cases using only parameter means. The magnitude of highly nonlinear, threshold-controlled processes such as wind erosion and dust flux is essentially controlled by variability in the landscape. This result argues for inclusion of surface heterogeneity in global and regional dust emission models.