Sedimentary formations that compose most aquifers are difficult to model as a result of the nature of their deposition. Their formation generally involves multiple processes (alluvial, glacial, lacustrine, etc.) that contribute to the complex organization of these deposits. Representative models can be obtained using process-based or rule-based methods. However, such methods have several drawbacks: complicated parametrization, large computing time, and challenging, if not impossible, conditioning. To address these problems, we propose a new hierarchical surface-based algorithm, named EROSim. First, a predefined number of stochastic surfaces are simulated in a given order (from older to younger). These surfaces are simulated independently but interact with each other through erosion rules. Each surface is either an erosive or a deposition surface. The deposition surfaces represent the boundaries of depositional events, whereas the erosive surfaces can remove parts of the previously simulated deposits. Finally, these surfaces delimit sedimentary regions that are filled with facies. The approach is quite simple, general, flexible, and can be conditioned to borehole data. The applicability of the method is illustrated using data from fluvio-glacial sedimentary deposits observed in the Bümberg quarry in Switzerland.