Abstract Viscous folding in mechanically heterogeneous layers is modelled numerically in two dimensions for linear and power-law viscous fluids. Viscosity heterogeneities are expressed as circular-shaped variations of the effective viscosity inside and outside the layers. The layers are initially perfectly flat and are shortened in the layer-parallel direction. The viscosity heterogeneities cause a perturbation of the velocity field from the applied bulk pure shear, which perturb geometrically the initially flat-layer interfaces from the first numerical time step. This geometrical perturbation triggers interfacial instabilities, resulting in high-amplitude folding. We compare simulations with heterogeneities with corresponding simulations in which the heterogeneities are removed after the first time step, and, hence, only the initial small geometrical perturbations control wavelength selection and high-amplitude folding. Results for folding in heterogeneous and homogeneous layers are similar, showing that viscosity heterogeneities have a minor to moderate impact on fold wavelength selection and high-amplitude folding. Our results indicate that the interfacial instability is the controlling process for the generation of buckle folds in heterogeneous rock layers. Therefore, existing analytical and numerical solutions for folding in homogeneous layers, in which folding was triggered by geometrical perturbations, are useful and applicable to study folding in natural, heterogeneous rock layers.