Stochastic simulation of soil properties serves multiple purposes, including filling in missing values with most probable alternatives, or extrapolating both in space and across scales. Using multifractal characteristics in the stochastic simulation of soil properties has a potential for producing simulations that better represent features of interest, including distribution properties of the high/low values, than the simulations based on geostatistical characteristics. The objective of this study is to compare performance of stochastic simulations that reproduce certain multifractal characteristics with simulations that reproduce variogram model parameters. Data on soil phosphorus concentrations and soil organic matter contents collected at approximately 1,700 sampling locations were used in the simulations. Phosphorus data represented an example of a variable exhibiting a multifractal scaling and organic matter an example of a variable exhibiting a monofractal scaling. In both data sets the data were separated into “model” and “test” data sets. The “model” data set included samples located on a 12x14 regular grid with 100 m distance between the grid points and additional 20 samples randomly selected such that the distance from them to the grid points was equal to 50 m. The remaining 836 samples constituted the “test” data set. The simulations were conducted using a simulated annealing procedure with conditioning the simulations on the “model” data sets. The initial random values of the simulated annealing are drawn from the population distribution of the soil property, which in our example was constructed based on the “model” data. Then, the structure functions of the initial simulated data set are calculated for q values ranging from 1 to 5 in 1.0 increments. The structure function exponents that describe the scaling properties are then obtained for each q value by linear regressions performed on the log-transformed values. The obtained set of exponents is compared with the exponent values desired as characteristics for the final simulated field. As the desired characteristics I used the exponent values obtained based on the whole data set, assuming them to be true structure function exponents of the studied field. The simulation process consists of perturbing the data set followed by calculations of structure function exponent values and assessing the improvement as compared with the desired characteristics. The process continues until the structure function exponent values of the perturbed data set closely match those of the desired spatial structure. To perform the simulations, the program SASIM from the GSLIB package (Deutsch and Journel, 1998) was modified by adding the structure function exponents to the procedure's objective functions.Two simulation scenarios were compared. The first scenario consisted of using the histogram and the structure function for q equal to 2, which is equivalent to simulating a data set with a certain variogram. The second scenario used the histogram and the structure function with exponent values for q ranging from 1 to 5. The first scenario is equivalent to assuming that the studied variable is a monofractal, while the second scenario assumes the variable to be a multifractal. Based on 50 simulations the performances of the two scenarios in representation of the highest/lowest values of the “test” data set were compared. Out of a total of 836 test data points, the comparisons of the performances were done, first, by observing how many of the highest and lowest 46 observations (10th quantile and 95th quantile) were accurately predicted in the simulated field, and, second, by how many of the 80 highest (90th quantile) and 130 lowest (20th quantile) observations were accurately predicted. There was no significant difference in performance of monofractal and multifractal simulations for soil organic matter. Such a result was expected since organic matter data of this study were not exhibiting multifractal scaling. However, performance of monofractal and multifractal simulations were different for soil phosphorus. Multifractal simulations performed better than monofractal in correct predictions of the locations with low values. Multifractal approach produced significantly more accurate predictions of the lowest 46 observations and the lowest 130 values (p<0.1), but significantly less accurate predictions of the 80 highest observations. There was no difference between performance of the monofractal and multifractal simulations in prediction accuracy for the highest 46 phosphorus values. The results indicate that when accurate representations of extreme values of the studied soil variables are of practical importance and the studied soil variables are expected to exhibit multifractal scaling, the stochastic simulations based on the multifractal characteristics should be preferred to those based on variogram parameters.