Radio‐echo sounding (RES) shows large‐scale englacial stratigraphic folds are ubiquitous in Greenland's ice sheet. However, there is no consensus yet on how these folds form. Here, we use the full‐Stokescode Underworld2 to simulate ice movements in three‐dimensional convergent flow, mainly considering iceanisotropy due to a crystallographic preferred orientation, vertical viscosity and density gradients in ice layers,and bedrock topography. Our simulated folds show complex patterns and are classified into: large‐scale folds(>100 m amplitude), small‐scale folds (<<100 m) and basal‐shear folds. The amplitudes of large‐scale foldstend to be at their maximum in the middle of the ice column or just below, in accordance with observations inRES data. We conclude that ice anisotropy amplifies the perturbations in ice layers (mainly due to bedrock topography) into large‐scale folds during flow. Density differences between the warm deep ice and cold iceabove may enhance fold amplification.