Bedrock undulations of glaciers and ice sheets cause undulations of the ice surface and of internal layers. Models are constructed, that describe this phenomenon, using perturbation techniques and utilizing the stress function as harmonic solution to the full biharmonic equation. This approach is an extension of a previous treatment of the problem, which was based on the solution of the Laplace equation. This solution is shown to be inadequate. Two types of ice flow, with and without bottom sliding, are considered. Furthermore, in order to cope with the variation of viscosity with depth, semi-analytical multi-layer models have been constructed.