In order to understand the interactions between surface processes and multilayer folding systems, we here present fully coupled three-dimensional numerical simulations. The mechanical model represents a sedimentary cover with internal weak layers, detached over a much weaker basal layer representing salt or evaporites. Applying compression in one direction results in a series of three-dimensional buckle folds, of which the topographic expression consists of anticlines and synclines. This topography is modified through time by mass redistribution, which is achieved by a combination of fluvial and hillslope erosion, as well as deposition, and which can in return influence the subsequent deformation. Model results show that surface processes do not have a significant influence on folding patterns and aspect ratio of the folds. Nevertheless, erosion reduces the amount of shortening required to initiate folding and increases the exhumation rates. Increased sedimentation in the synclines contributes to this effect by amplifying the fold growth rate by gravity. The main contribution of surface processes is rather due to their ability to strongly modify the initial topography and hence the initial random noise, prior to deformation. If larger initial random noise is present, folds amplify faster, which is consistent with previous detachment folding theory. Variations in thickness of the sedimentary cover (in one or two directions) also have a significant influence on the folding pattern, resulting in linear, large aspect ratio folds. Our simulation results can be applied to folding-dominated foldand-thrust belt systems, detached over weak basal layers, such as the Zagros Folded Belt.