The edge function method is considered for the analysis of plane strain problems in glacier mechanics. The essence of the approach is the approximation of the solution by a linear combination of analytical solutions (based on the complex variable formulation of anisotropic elasticity) of the field equations. The unknowns in the linear combination are obtained from a system of equations which follows from the approximation of the boundary conditions by a boundary Galerkin energy method. A range of representative glacier geometries are examined and the presence of areas of tension on the surface is evident. The introduction of small crevasses in the form of widely spaced notches does not seem to have a significant effect in the reduction of that tension, away from the notches. A comparison with a finite element approach indicates that an accurate solution is obtained from the edge function method with fewer degrees of freedom and reduced setup effort. The method is well suited to capturing the singularity at the tip of the crevasse but is limited to wide crevasses and the examination of realistic narrow crevasses would be handled more optimally by a combination of the edge function method with its local accuracy and a conventional method with more robust general geometry capabilities.