The Navier–Stokes equation, integrated vertically to yield a two-dimensional transport equation, is combined with the three-dimensional continuity equation. By assuming a linear relation between volume transport and bottom shear stress, a system of numerically tractable differential equations is developed that contains a nonlinear surface slope term and a term in which horizontal frictional forces are stated explicitly. These equations, which are not based on perturbation analysis, are integrated directly on a high speed digital computer. The coefficient of each term in the flow equations is composed of physical parameters (gravity, density, viscosity); therefore, since no quantities such as ice velocities enter the equations, it is only necessary to specify an arbitrary bed configuration and a net-balance versus altitude curve to obtain a glacier solution.For two beds, valley and cirque, steady-state glacier solutions are found for a variety of bed and lateral frictional values. By reducing the bed friction coefficient for one year over the entire bed of the glacier, surges are induced in each of these cases. A reduction in bed friction to 5% of its original value yields a surge resembling typical observed surges. Interesting wave forms occur during the recovery. During the initial recovery, after the bed friction has been restored to its original value, height falls continue to occur even in the lower accumulation zone. The original steady-state shape is approached asymptotically, with total recovery taking hundreds of years.