We use a new discretization technique to solve the higher-order thermomechanically coupled equations of glacier evolution. We find that under radially symmetric continuum equations, small perturbations in symmetry due to the discretization are sufficient to produce the initiation of nonsymmetric thermomechanical instabilities which we interpret as ice streams, in good agreement with previous studieswhich have indicated a similar instability. We find that the inclusion of membrane stresses regularizes the size of predicted streams, eliminating the ill-posedness evident in previous investigations of ice stream generation through thermomechanical instability. Ice streams exhibit strongly irregular periodicity which is influenced by neighboring ice streams and the synoptic state of the ice stream. Ice streams are not always the same size but instead appear to follow a temperature-dependent distribution of widths that is robust to grid refinement. The morphology of the predicted ice streams corresponds reasonably well to extant ice streams in physically similar environments.