The extent of remediation of contaminated industrial sites depends on spatial heterogeneity of contaminant concentration and spatially explicit risk characterization. We used sequential Gaussian simulation (SGS) and indicator kriging (IK) to describe the spatial distribution of polycyclic aromatic hydrocarbons (PAHs), pH, electric conductivity, particle aggregate distribution, water holding capacity, and total organic carbon, and quantitative relations among them, in a creosote polluted soil in southern Sweden. The geostatistical analyses were combined with risk analyses, in which the total toxic equivalent concentration of the PAH mixture was calculated from the soil concentrations of individual PAHs and compared with ecotoxicological effect concentrations and regulatory threshold values in block sizes of 1.8 x 1.8 m. Most PAHs were spatially autocorrelated and appeared in several hot spots. The risk calculated by SGS was more confined to specific hot spot areas than the risk calculated by IK, and 40-50% of the site had PAH concentrations exceeding the threshold values with a probability of 80% and higher. The toxic equivalent concentration of the PAH mixture was dependent on the spatial distribution of organic carbon, showing the importance of assessing risk by a combination of measurements of PAH and organic carbon concentrations. Essentially, the same risk distribution pattern was maintained when Monte Carlo simulations were used for implementation of risk in larger (5 x 5 m), economically more feasible remediation blocks, but a smaller area became of great concern for remediation when the simulations included PAH partitioning to two separate sources, creosote and natural, of organic matter, rather than one general.