Simulation of subsurface heterogeneity is important for modeling subsurface flow and transport processes. Previous studies have indicated that subsurface property variations can often be characterized by fractional Brownian motion (fBm) or (truncated) fractional Levy motion (fLm). Because Levy-stable distributions have many novel and often unfamiliar properties, studies on generating fLm distributions are rare in the literature. In this study, we generalize a relatively simple and computationally efficient successive random additions (SRA) algorithm, originally developed for generating Gaussian fractals, to simulate fLm distributions. We also propose an additional important step in response to continued observations that the traditional SRA algorithm often generates fractal distributions having poor scaling and correlation properties. Finally, the generalized and modified SRA algorithm is validated through numerical tests.