Many emerging, and some legacy, pollutants pose risks to humans and ecosystems near the detection limits (DL) of existing analytical systems. As a result, site assessments and management options are often presented with data sets that are sparse, highly skewed, and left-censored. Existing analysis methods are unable to differentiate effects of treatment from covariates, such as space, obscuring influences of site management. As a case study, we computed the mean and variance of censored soil benzene data across four sites over a three year period by gamma distribution with a maximum likelihood. Further, a combined hurdle model to accommodate left-censored concentrations was applied to analyze factors affecting benzene variation. This approach allowed us to assess the success and spatial dependency of a biostimulatory solution in reducing benzene concentrations at very low concentrations. Benzene concentrations decreased due to the addition of biostimulatory solution and spatial effects, but the detection of soil benzene after biostimulation was highly spatially dependent. By combining computed values for censored observations estimated by the hurdle-gamma model and uncensored observations, we can get the pseudocomplete data sets. The combined model is ideally suited to evaluate existing and emerging pollutants, that pose risks to humans and ecosystems but are typically at or near analytical detection limits.