Radio-echo sounding (RES) has shown that large-scale folds in the englacial stratigraphy is ubiquitous in Greenland’s ice sheet. However, there is no consensus yet on how these folds form. Here, we use the full-Stokes code Underworld2 to simulate ice movements in three-dimensional convergent flow, mainly investigating the effect of ice anisotropy due to a crystallographic preferred orientation, vertical viscosity and density contrasts in ice layers, and bedrock topography. Our simulated folds show complex patterns and are classified into three types: large-scale folds, small-scale folds and basal-shear folds. The amplitudes of large-scale folds tend to be at their maximum in middle ice layers and decrease towards the surface, in accordance with observations in RES data. We conclude that bedrock topography contributes to perturbations in ice layers, and that ice anisotropy amplifies these into large-scale folds, while vertical viscosity contrasts in ice layers are insufficient for large-scale fold amplification.