Compression of the lithosphere, sedimentary sequences or quartz veins may result in a folding instability. We perform numerical simulations of viscous single-layer folding to study this instability in 3D. It is demonstrated that linear theories correctly describe the instability for small amplitudes. At larger amplitudes, however, the theory breaks down. For these stages we present a new nonlinear amplification equation. Numerical simulations of folding of an initially horizontal layer, perturbed with random noise, demonstrate that in most cases fold axes form perpendicular to the main shortening direction. Aspect ratios of folds are finite and the patterns are relatively insensitive to the applied background shortening directions. Furthermore, the 3D folding instability reduces the averaged differential stress within the folded ("strong") layer, in agreement with 2D results. This implies that the Christmas-tree approach to represent the strength of the crust and lithosphere may be invalid if folding occurs during the deformation.